As coupling flows need to split coordinates either across particles or across spatial dimensions or a combination thereof, the conditioner C does not have access to complete three-dimensional interatomic distances, as its input can only consist of a subset of the coordinates40<\/a> (see Supplementary Note\u00a01<\/a> for a detailed discussion on the limitations of standard coupling flows). Consequently, previous architectures often rely on the absolute coordinates of all atoms in the system15<\/a>,19<\/a>, and we will refer to this approach as a \u2018global\u2019 architecture in the following. To overcome this architectural limitation of coupling flows, we leverage the augmented coupling flow framework27<\/a>. Augmented flows introduce auxiliary variables \\({{{\\bf{a}}}}\\in {{\\mathbb{R}}}^{3N}\\) and enable the splitting to be performed between physical and auxiliary variables, thereby updating the physical variables \\({{{\\bf{x}}}}\\in {{\\mathbb{R}}}^{3N}\\) conditioned on a<\/b> and vice versa. Importantly, this scheme retains full three-dimensional coordinates within the physical and auxiliary space and, thus, allows for computing interatomic distances. We model the auxiliary system as a noised copy of the physical system41<\/a> and define the joint base distribution of the augmented system as <\/p>\n $$\\,q({{{\\bf{x}}}},{{{\\bf{a}}}})=q({{{\\bf{x}}}})\\,{{{\\mathcal{N}}}}({{{\\bf{a}}}};{{{\\bf{x}}}},{\\eta }_{q}^{2}{{{\\bf{I}}}})\\,.$$<\/p>\n \n (1)\n <\/p>\n Here, q(x<\/b>) is the base distribution of the physical variables and \\({{{\\mathcal{N}}}}(\\cdot ;{{{\\bf{x}}}},{\\eta }_{q}^{2}{{{\\bf{I}}}})\\) denotes the normal distribution centered at x<\/b> with covariance matrix \\({\\eta }_{q}^{2}{{{\\bf{I}}}}\\), where I<\/b> is the identity matrix and \\({\\eta }_{q}\\in {{\\mathbb{R}}}_{\\, > \\,0}\\). In this construction, the auxiliary system retains structural information of the physical variables that can be exploited by the flow transformation.<\/p>\n Utilizing the auxiliary variables, we define the flow transformation for particle i by <\/p>\n $${{{{\\bf{x}}}}}_{i}^{{\\prime} }=g({{{{\\bf{x}}}}}_{i}| {{{{\\bf{h}}}}}_{i}^{{{{\\bf{a}}}}})\\,,$$<\/p>\n \n (2)\n <\/p>\n where g is parametrized by the embedding \\({{{{\\bf{h}}}}}_{i}^{{{{\\bf{a}}}}}\\), which is computed using a graph neural network (GNN)26<\/a> that aggregates information from each auxiliary particle\u2019s local neighborhood (see Fig.\u00a01<\/a>). Details on the implementation can be found in Supplementary Note\u00a03<\/a>. A crucial feature of the local approach presented here is its linear scaling with system size, achieved by fixing the number of neighbors. This stands in stark contrast to global architectures based on, for example, attention mechanisms, which learn interactions between all particles, resulting in quadratic scaling with system size15<\/a>,42<\/a>.<\/p>\n Fig. 1: Size-transferable augmented flow.<\/b> a<\/b> A GNN learns environment-dependent particle embeddings, h<\/b>i, in a small system and is transferable across system size via each particle\u2019s local neighborhood, \\({{{{\\mathcal{N}}}}}_{i}\\). b<\/b> Physical and auxiliary particles are updated sequentially using the learned embeddings. In each step, either all auxiliary particles or all physical variables are updated.<\/p>\n To ensure that the bijector, which we implement as rational quadratic splines43<\/a>, remains applicable to larger system sizes, it is crucial that the magnitude of its input does not scale with the system size. One possibility that we follow in our proposed architecture is to model the displacements from the ideal crystal lattice rather than absolute positions. This makes the input size-independent and preserves the model\u2019s ability to generalize to arbitrary system sizes. Accordingly, we define the physical base distribution as \\(q({{{\\bf{x}}}})={{{\\mathcal{N}}}}({{{\\bf{x}}}};0,{\\eta }_{q}^{2}{{{\\bf{I}}}})\\).<\/p>\n For a flow acting on the augmented space, the change-of-variable for the augmented system writes <\/p>\n $${q}_{\\theta }({{{{\\bf{x}}}}}^{{\\prime} },{{{{\\bf{a}}}}}^{{\\prime} })=q({{{\\bf{x}}}},{{{\\bf{a}}}})\\left| \\det {J}_{{{{{\\mathcal{F}}}}}_{\\theta }}({{{\\bf{x}}}},{{{\\bf{a}}}})\\right|{\\scriptstyle{{-1}}\\atop} \\!\\,,$$<\/p>\n \n (3)\n <\/p>\n where \\(\\det {J}_{{{{{\\mathcal{F}}}}}_{\\theta }}\\) denotes the Jacobian determinant of the flow transformation \\({{{{\\mathcal{F}}}}}_{\\theta }\\). The target distribution of the auxiliary variables is \\(\\pi ({{{\\bf{a}}}}| {{{\\bf{x}}}})={{{\\mathcal{N}}}}({{{\\bf{a}}}};{{{\\bf{x}}}},{\\eta }_{p}^{2}{{{\\bf{I}}}})\\), such that the flow is trained to optimize the joint target distribution <\/p>\n $$p({{{\\bf{x}}}},{{{\\bf{a}}}})=p({{{\\bf{x}}}})\\,{{{\\mathcal{N}}}}({{{\\bf{a}}}};{{{\\bf{x}}}},{\\eta }_{p}^{2}{{{\\bf{I}}}})\\,,$$<\/p>\n \n (4)\n <\/p>\n where p(x<\/b>) is the equilibrium distribution of the ensemble of interest. Here, we focus on the canonical (or NVT) ensemble, whose equilibrium distribution is the Boltzmann distribution. The flow is trained by minimizing the Kullback-Leibler (KL) divergence between the generated and target distributions in the joint space. More details on the canonical ensemble and the explicit form of the loss function can be found in \u201cMethods\u201d.<\/p>\n Within the augmented setting, only the joint generated density \\({q}_{\\theta }({{{{\\bf{x}}}}}^{{\\prime} },{{{{\\bf{a}}}}}^{{\\prime} })\\) can be evaluated exactly. However, the marginal of the physical system can be obtained by integrating out the auxiliary degrees of freedom, which can be approximated as41<\/a><\/p>\n $${q}_{\\theta }({{{{\\bf{x}}}}}^{{\\prime} })\\approx \\frac{1}{M}{\\sum }_{m=1}^{M}\\frac{{q}_{\\theta }({{{{\\bf{x}}}}}^{{\\prime} },{{{{\\bf{a}}}}}_{m})}{\\pi ({{{{\\bf{a}}}}}_{m}| {{{{\\bf{x}}}}}^{{\\prime} })}\\,,$$<\/p>\n \n (5)\n <\/p>\n with \\({{{{\\bf{a}}}}}_{m} \\sim \\pi (\\cdot | {{{{\\bf{x}}}}}^{{\\prime} })\\) and M denoting the number of auxiliary samples drawn per generated physical sample. Unless stated otherwise, all results obtained from the marginal densities presented in this work were obtained using M\u00a0=\u00a0200, for which the free energy estimates from the marginal density were found to be converged (see Supplementary Note\u00a01<\/a>). Importantly, while the evaluation of the marginal distribution does require M additional inverse passes through the flow per generated sample (since \\({q}_{\\theta }({{{{\\bf{x}}}}}^{{\\prime} },{{{{\\bf{a}}}}}^{{\\prime} })\\,=\\,q({{{{\\mathcal{F}}}}}_{\\theta }^{-1}({{{{\\bf{x}}}}}^{{\\prime} },{{{{\\bf{a}}}}}^{{\\prime} }))| \\det {J}_{{{{{\\mathcal{F}}}}}_{\\theta }^{-1}}({{{{\\bf{x}}}}}^{{\\prime} },{{{{\\bf{a}}}}}^{{\\prime} })|\\)), it does not require additional evaluations of the target potential.<\/p>\n As detailed in \u201cMethods\u201d, the calculation of the (reduced) Helmholtz free energy associated to a distribution \u03bc, \\({f}_{\\mu }={F}_{\\mu }\/{k}_{\\mathrm {B}}{T}_{\\mu }=-\\log {Z}_{\\mu }\\), requires an estimate of the partition function Z\u03bc. While this quantity cannot be computed directly, the trained flow can be used within targeted free energy perturbation (TFEP)44<\/a> to estimate the free energy difference between the base and target states, \u0394fqp\u00a0=\u00a0fp\u00a0\u2212\u00a0fq. Within the augmented flow framework, the partition function of the generated joint distribution depends on both physical and auxiliary variables. Since the target and base distributions of the auxiliary system are Gaussians, the partition function factorizes (see Supplementary Note\u00a01<\/a>) and the free energy of an augmented system is given by <\/p>\n $${f_{\\mu }^{{{{\\rm{aug}}}}}}={f_{\\mu }}+{f_{\\mu }^{{{{\\rm{aux}}}}}}\\,,$$<\/p>\n \n (6)\n <\/p>\n where f\u03bc and \\({f}_{\\mu }^{{{{\\rm{aux}}}}}\\) denote the free energy of the physical and auxiliary system, respectively. The corresponding TFEP estimator for the augmented system is <\/p>\n $$\\Delta {f_{qp}^{{{{\\rm{aug}}}}}}\t={f_{p}^{{{{\\rm{aug}}}}}}-{f_{q}^{{{{\\rm{aug}}}}}}={f_{p}}+{f_{p}^{{{{\\rm{aux}}}}}}-({f_{q}}+{f_{q}^{{{{\\rm{aux}}}}}})\\\\ \t=\\Delta {f_{qp}}+\\Delta {f_{qp}^{{{{\\rm{aux}}}}}}\\,.$$<\/p>\n \n (7)\n <\/p>\n If \u03b7q\u00a0=\u00a0\u03b7p, as chosen in this work, then \\(\\Delta {f_{qp}^{{{{\\rm{aux}}}}}}=0\\) and \\(\\Delta {f_{qp}^{{{{\\rm{aug}}}}}}=\\Delta {f_{qp}}\\).<\/p>\n Flow-based approaches also allow to compute Gibbs free energies corresponding to the NPT ensemble, in which the shape of the simulation box, \\({{{\\bf{h}}}}\\in {{\\mathbb{R}}}^{3\\times 3}\\) with \\(V=\\det {{{\\bf{h}}}}\\), changes during the simulation at constant pressure P1<\/a>. While prior work17<\/a>,18<\/a> studied flows modeling shape and configurational distributions simultaneously, in our current approach, the Gibbs free energy is computed via a Legendre transformation as \\(G={\\min }_{{{{\\bf{h}}}}}[F({{{\\bf{h}}}})+\\det ({{{\\bf{h}}}})P]\\)45<\/a>. To obtain a reliable estimate for F(h<\/b>), we train the flow in a shape-conditional fashion18<\/a>,46<\/a>, optimizing the flow to generate samples for varying box dimensions. The Gibbs free energies computed in this manner also enable the calculation of derived quantities, such as isothermal compressibilities (see Supplementary Notes\u00a03<\/a> and 4<\/a> for details on the conditional training and derived quantities, respectively).<\/p>\n The principle advantages of our proposed BG architecture based on local augmented flows are reduced training times, improved sampling efficiencies, and, most importantly, its ability to generalize to larger systems. In the following, we demonstrate the applicability of our method across a range of crystal structures for various materials systems, exemplified by the Lennard-Jones (LJ) potential and multiple parameterizations of the Stillinger-Weber potential28<\/a>, including ice phases of monatomic water (mW) and silicon. Computational details are provided in\u00a0\u201cMethods\u201d.<\/p>\n Size transferable training and evaluation<\/p>\n We illustrate the effectiveness of the size transferability by applying local BGs, trained on cubic mW ice with N\u00a0=\u00a0216 particles and face-centred cubic (FCC) LJ with N\u00a0=\u00a0256, to systems containing 512, 1000, and 1728 particles for cubic mW ice and 500, 864, 1000, and 1372 particles for FCC LJ.<\/p>\n Figure\u00a02<\/a> shows the radial distribution functions (RDFs) and potential energy histograms for the largest investigated system sizes, as computed from samples generated by the BGs, MD, and drawn from the base distribution. For both systems, the RDFs produced by the BGs are nearly indistinguishable from those obtained via MD over length scales going far beyond the size of the training systems. While by design the RDF of the base distribution already captures the main features of the target, the corresponding configurations do not reflect the correct correlations between particle positions in the target potential, which is evident from the pronounced differences in the energy distributions of the base and MD ensembles. In contrast, the local BGs accurately reproduce both structural and energetic features in both systems, yielding excellent agreement in the histograms already without reweighting. These results clearly demonstrate that coordinate transformations based on local environments are entirely sufficient for generating accurate equilibrium structures and are reliably transferable to larger system sizes.<\/p>\n Fig. 2: Radial distribution functions and energy histograms.<\/b> a<\/b>, c<\/b> show the radial distribution functions (RDFs) in units of the interaction length \u03c3 for FCC LJ with N\u00a0=\u00a01372 and cubic mW ice with N\u00a0=\u00a01728, respectively, as obtained from MD, the base distributions, and the local BGs. b<\/b>, d<\/b> show the corresponding reduced energy histograms. The local BG results were computed using models trained with N\u00a0=\u00a0216 (mW) and N\u00a0=\u00a0256 (LJ). The red line indicates half the box length of the training system, which is the largest distance for which the RDF can be evaluated for this system. For the transferred systems, the RDF is similarly evaluated up to their corresponding half-box length. No reweighting was applied. Source data are provided as a Source Data file.<\/p>\n While structural and energetic accuracy are necessary conditions for accurate free energy estimates, they alone do not guarantee high sampling efficiency or precise free energy calculations. A more rigorous measure is the effective sample size (ESS) (see Methods) evaluated across the various system sizes within the transferable framework and summarized in the lower part of Table\u00a01<\/a>. As the system size increases, the ESS naturally decreases, as even constant per-particle errors in the extensive potential energy lead to exponentially amplified errors in the Boltzmann distribution. But even for large systems such as cubic mW ice with N\u00a0=\u00a01728 and FCC LJ with N\u00a0=\u00a01372, the ESS of the local BGs is large enough to ensure reliable and highly accurate statistics of ensemble averages.<\/p>\n Table 1 Effective sample size (ESS) in percent for local and global BGs for cubic mW ice and FCC LJ with different number of particles<\/b><\/p>\n Importantly, the ESS for large systems is significantly higher than that obtained with global approaches, while the computational cost is significantly reduced. For mW ice with N\u00a0=\u00a0512, global BGs reported in refs. 15<\/a>,\u200917<\/a> required more than 330 GPU days of training until convergence, yet achieved ESS values of only around 0.2%. In contrast, for the mW system with N\u00a0=\u00a0216, our local BG reaches convergence in approximately four GPU days on the same hardware (see Supplementary Note\u00a06<\/a> for details on training times), while yielding ESS values for the 512-particle system up to an order of magnitude higher than those obtained with the global BG. Moreover, scaling global architectures to larger systems becomes computationally prohibitive, making particle counts beyond 500 effectively impossible due to excessive training costs and deteriorating sampling efficiency.<\/p>\n The improved training performance of the local architecture compared to a global one is further illustrated in Fig.\u00a03<\/a>, evaluating the ESS of both types as a function of training time for three different systems (cubic mW ice with N\u00a0=\u00a064 and N\u00a0=\u00a0216, and FCC LJ with N\u00a0=\u00a0216). For the global BGs, the architecture proposed in Ref. 15<\/a> was used, which employs an attention-based coupling flow. Across all systems, the local BGs consistently achieve higher ESS than the global ones at any given training time. This is particularly evident for the larger system sizes, demonstrating the superior scaling of the local BGs with number of particles. The local BGs also exhibit substantially lower variance in ESS, indicating that reliable estimates can be obtained even with small sample sizes and remain stable as the number of samples increases.<\/p>\n Fig. 3: Effective sample size.<\/b> Effective sample size (ESS) plotted against the number of training steps for global and local BGs of a<\/b> cubic mW ice with N\u00a0=\u00a064 and b<\/b> N\u00a0=\u00a0216, and c<\/b> FCC LJ with N\u00a0=\u00a0256. The ESS was evaluated every 1k steps using 1k samples. Solid lines represent running averages of ESS across five independent training runs, with thin, light-colored lines indicating the raw ESS values from individual runs. For the local BGs, the joint efficiency is shown. Learning rate reductions were applied after 250k and 500k steps. Source data are provided as a Source Data file.<\/p>\n These trends are corroborated by the ESS obtained from the converged models, which are reported in the upper part of Table\u00a01<\/a>. Remarkably, the marginal ESS of the physical system, which is the metric of primary interest, exceeds 80% for mW ice with N\u00a0=\u00a064 and over 40% for N\u00a0=\u00a0216, and reaches about 7% for FCC LJ with N\u00a0=\u00a0256, highlighting the significantly improved efficiency of our local flow architecture in capturing the relevant short-range environmental information. Where available, we have also provided the ESS of the global BGs as reported in Ref. 17<\/a>, which were obtained using models comprising twice as many layers and trained significantly longer. While the final ESS reached by these models is slightly higher than the ones we obtained for the global BGs within this work, they remain far lower than those of the local BGs.<\/p>\n We note that the BGs do not perform equally well across all systems. For example, both local and global BGs achieve significantly higher sampling efficiencies for the cubic mW ice compared to FCC LJ, despite only a modest increase in particle number. However, we also observed that certain crystal structures within the same potential can exhibit rather different sampling efficiencies, which asks for further investigation beyond the scope of this work.<\/p>\n Free energy estimation<\/p>\n A key advantage of a trained BG is that it enables a direct evaluation of free energy differences within the framework of TFEP44<\/a>, as discussed in \u201cMethods\u201d. Utilizing our local architecture, BGs trained on relatively small systems allow for an accurate evaluation of absolute Helmholtz and Gibbs free energies for far larger systems, which reduces the computational effort dramatically compared to traditional free energy estimators such as MBAR47<\/a>.<\/p>\n The left panel of Fig.\u00a04<\/a> shows the absolute free energy of the cubic ice phase as a function of the number of particles computed with the BG trained on 216 particles. In addition to the joint free energy estimates, we also present values computed from the marginal generated density and further compare to high-accuracy MBAR estimates obtained from MD simulations interpolating between the Einstein crystal and the physical system of interest48<\/a>. Joint and marginal free energy estimates from the BG show excellent agreement with the reference values for system sizes up to 1000 particles, deviating by less than 10\u22123kBT per particle. For the largest system size (N\u00a0=\u00a01728), the joint estimates slightly overestimate the true values, while the marginal evaluation partially corrects for this, yielding more accurate results in agreement with the ESS reported in Tab.\u00a01<\/a>. In this context, it is worth mentioning that it is difficult to verify the reliability of a given free energy estimate at low ESS, since the variance of the estimates can remain very small even when ESS is low (see Supplementary Note\u00a05<\/a>). Based on the examples considered, a mean ESS of 0.1% or higher appears to be a sufficient criterion to ensure deviations in the free energy estimates remain below 10\u22123kBT per particle when using 5 \u22c5 104 samples. A similar behavior was observed for the FCC LJ crystal, with the corresponding free energy estimates reported in Supplementary Note\u00a05<\/a>.<\/p>\n Fig. 4: Helmholtz Free energy estimates.<\/b> a<\/b> Absolute reduced Helmholtz free energy estimates per particle for the cubic ice systems against the particle number as obtained from the local BG. Blue and orange lines indicate the evaluation of joint and marginal densities, respectively. Reference MD+MBAR values are shown as a gray shaded area corresponding to their mean \u00a0\u00b1\u00a010\u22123kBT. The red shaded area corresponds to the MD+MBAR results evaluated at N\u00a0=\u00a04096. BG results were obtained using a model trained at N\u00a0=\u00a0216. Uncertainties are smaller than the marker size (see Supplementary Note\u00a05<\/a>). b<\/b> Reduced free energy difference between cubic and hexagonal ice, \u0394f\u00a0=\u00a0fhex\u00a0\u2212\u00a0fcubic, against the particle number. Shaded areas have the same meaning as for the left plot. Source data are provided as a Source Data file.<\/p>\n Obtaining accurate free energy difference between different crystalline phases is particularly challenging but crucial to explore phase diagrams. To reach a converged value of the free energy difference between cubic and hexagonal mW ice, system sizes up to 1728 particles are necessary, as shown in the right panel of Fig.\u00a04<\/a>. Both joint and marginal estimates of \u0394f are in very close agreement with the reference values, where the good performance of the joint estimates is likely due to error cancellation between the two phases. The large number of particles necessary to fully converge the free energy difference between the cubic and hexagonal phase underscores the importance of being able to evaluate extended system size which is infeasible with global architectures.<\/p>\n To quantify the computational savings associated with the size-transferable free energy estimation, we compare the number of energy evaluations required to the one of the standard MD+MBAR approach using typical parameters to interpolate between the Einstein crystal and the physical system21<\/a>,49<\/a>,50<\/a>. The MBAR setup requires approximately 50 intermediate states with 103 uncorrelated samples per state, obtained by storing one sample every 103 MD steps, resulting in a total of around 5 \u22c5 107 energy evaluations. In addition, MBAR evaluates all 5 \u22c5 104 stored samples across all 50 intermediate potentials, resulting in additional 2.5 \u22c5 106 energy evaluations. Training the flow model until convergence involves a comparable number of evaluations, approximately 6 \u22c5 107 for 5 \u22c5 105 training steps with batches of 128 samples (see Fig.\u00a03<\/a>). However, the training is performed for small system sizes only. Once trained, evaluating the free energy with the BG requires only a few tens of thousands of samples and the corresponding energy evaluations for the large system. In contrast, MD+MBAR necessitates several tens of millions of energy evaluations for each system size, with the computational cost scaling quadratically for pairwise potentials and even more steeply for quantum mechanical methods such as density functional theory and higher level electronic structure methods51<\/a>. The computational savings of the flow-based approach are further amplified when performing convergence checks, which can be accomplished with a single model. Even further efficiency gains are possible by conditioning the model on external parameters, enabling reuse across multiple thermodynamic conditions18<\/a>, as exemplified below.<\/p>\n While BG models are highly efficient compared to MD+MBAR in terms of the number of required energy evaluations, this advantage does not directly translate to wall clock time for a single free energy estimate. For the cubic ice system studied here with N\u00a0=\u00a01728, a full MD+MBAR free energy estimate only takes a few hours, including simulations at all intermediate states and the MBAR evaluation, whereas training the flow model to convergence required approximately 8 GPU days (see Supplementary Note\u00a06<\/a> for details). However, this comparison is specific to inexpensive interaction models such as mW, which have a cost per MD step on the order of 10\u22127 s per atom. Machine learning interatomic potentials, today the workhorse of computational materials science, are typically more expensive than empirical potentials by a factor of 103 to 10452<\/a>. In these regimes, energy evaluations dominate the cost, while the overhead from optimization and neural network evaluation for the flow model remains largely independent of the potential. Consequently, the reduced number of energy evaluations required by the BG approach is expected to make the relative wall clock cost considerably more favorable, or even cheaper, than conventional MD+MBAR calculations.<\/p>\n Including volume fluctuations<\/p>\n To illustrate the efficiency and accuracy of our approach in the isothermal-isobaric ensemble, we determine the Gibbs free energy difference between FCC and the hexagonal close-packed (HCP) LJ crystals at constant pressure. Similar to the hexagonal and cubic ice phases, these structures are known to have extremely small free energy differences, which require large system sizes to converge. Crucially, the relative stability of FCC and HCP depends on the cutoff radius of the LJ potential and small cutoffs can even lead to qualitatively wrong predictions of the stable phase22<\/a>,53<\/a>. The local BGs for each phase are trained in a volume-conditional way (see Supplementary Note\u00a03<\/a>) using simulation cells comprising 180 particles. The trained BGs are subsequently applied to systems with N\u00a0=\u00a01080 to evaluate the Gibbs free energy over a range of cutoff radii.<\/p>\n Figure\u00a05<\/a> summarizes the BG-based Gibbs free energy estimates for the two crystalline phases. The left panel shows the difference in reduced Gibbs free energy per particle between HCP and FCC as a function of the cutoff radius. The predictions of the BGs (based on joint density estimates) show excellent agreement with reference values obtained from MD combined with MBAR. Consistent with previous studies22<\/a>,53<\/a>, a strong dependence of the free energy difference on the cutoff radius is revealed, leading to a qualitative change in the predicted stable phase. As the cutoff increases, this sensitivity decreases, emphasizing the need for large simulation cells to accommodate sufficiently large cutoff radii. Additionally, the particle densities obtained via the Legendre transformation using the BGs closely match the mean densities observed in NPT MD simulations, showing nearly identical dependence with respect to the cutoff radius (right panel of Fig.\u00a05<\/a>). Remarkably, the BG-based estimates yield highly accurate results over the entire range of cutoff radii. This strongly suggests that the local flow models can be trained using a single, short cutoff in a relatively small system, yet still be transferable to much larger systems while continuing to generate accurate configurations.<\/p>\n Fig. 5: Gibbs free energy estimates.<\/b> a<\/b> Reduced Gibbs free energy difference \u0394g\u00a0=\u00a0gHCP\u00a0\u2212\u00a0gFCC per particle between HCP and FCC crystal structures with 1080 particles in the LJ potential as a function of the cutoff radius. Shown are BG-based predictions (joint density estimates) alongside reference results from MD+MBAR. The red vertical line denotes the cutoff radius applied during training with 180 particles (quantities marked with an asterisk are expressed in LJ units, see Supplementary Note\u00a02<\/a>). b<\/b> Densities \u03c1* as obtained from the BGs (circles) and mean densities from MD NPT simulations (crosses) for FCC (upper panel) and HCP (lower panel). Error bars for both BG and MD+MBAR, estimated from three independent runs each, are smaller than the plotted marker size. Source data are provided as a Source Data file.<\/p>\n Conditioning on atom type<\/p>\n As another important application, we show that the flow can be trained conditioned on parameters of the interaction potential, which allows a single model to generalize over a whole class of systems. As an example, we employ the Stillinger-Weber (SW) potential28<\/a>, combining a pair potential and a three-body term. The strength of the three-body interaction is controlled by a factor \u03bb3, yielding a total potential energy function <\/p>\n $${U}_{{{{\\rm{SW}}}}}({{{\\bf{x}}}})={\\Phi }_{2}({{{\\bf{x}}}})+{\\uplambda }_{3}{\\Phi }_{3}({{{\\bf{x}}}})\\,.$$<\/p>\n \n (8)\n <\/p>\n For each system of interest, the SW potential is defined in terms of characteristic length and energy scales, which are tuned to reproduce specific properties of the system. However, when expressed in reduced units, the only effective parameter distinguishing different parameterizations is the strength of the three-body interaction, \u03bb354<\/a>. Importantly, different values of \u03bb3 stabilize or destabilize different crystal structures.<\/p>\n To study the phase diagram under varying three-body strengths and system sizes, we train local BGs for diamond cubic and \u03b2-tin structures with N\u00a0=\u00a0216, and for the body-centered cubic (BCC) structure with N\u00a0=\u00a0128. The BGs are trained conditioned on \u03bb3 \u2208 [15, 24] and the shape of the simulation box, allowing the model to generalize across the full class of SW potentials and cover an entire range of different materials, as well as to evaluate Gibbs free energies. Figure 6<\/a> shows the Gibbs free energies at zero pressure of the three different crystal structures as a function of \u03bb3 for the system sizes used during training (left panel) and transferred to N\u00a0=\u00a01000 for diamond cubic and \u03b2-tin and 686 for BCC (right panel). Across the entire range of \u03bb3, the free energy estimates of the BGs align very closely with reference values from MD+MBAR. While the transition point between the cubic and \u03b2-tin structures is largely independent of system size, a slight shift is observed in the transition point between the BCC and \u03b2-tin structures, moving towards lower values of the three-body interaction strength for larger system sizes. To capture this subtle effect, it is imperative to be able to determine highly accurate free energies for rather large system sizes. An additional advantage of the flow-based approach is that the models can be trained across all structures over the full \u03bb3 spectrum, even though some structures are only stable in certain parts of this region in MD simulations.<\/p>\n Fig. 6: Phase diagram SW potential.<\/b> a<\/b>, b<\/b> show the Gibbs free energies per particle for small and large particle numbers, respectively, as obtained from the local BGs at zero pressure for different crystal structures as a function of the three-body interaction strength \u03bb3. Vertical dashed lines indicate the value of this parameter in the parametrizations for germanium, silicon, and monatomic water. The color of the shaded areas corresponds to the most stable structure. MD+MBAR values are shown as black dashed lines. Error bars for both BG and MD+MBAR, estimated from three independent runs each, are smaller than the plotted line width. Source data are provided as a Source Data file.<\/p>\n Silicon phase diagram<\/p>\n Focusing on the SW parametrization for silicon (Si), we train local BGs conditioned on both temperature and the shape of the simulation box to map out the (T,\u00a0P)-phase diagram. Specifically, we are aiming to determine the coexistence line between the diamond cubic and \u03b2-tin phases. Both training and evaluation of the BGs were performed at N\u00a0=\u00a0216. While size-transferability was again excellent for the cubic structure, the \u03b2-tin structure appeared to be more challenging already at the training size, limiting the accuracy of free energy estimates in much larger systems. To improve the accuracy of the predicted phase diagram, the Gibbs free energies of the \u03b2-tin phase were evaluated based on the marginal density using M\u00a0=\u00a0200.<\/p>\n Even without transferring to larger systems, the conditional training provides substantial cost amortization across thermodynamic states, since reference MD+MBAR calculations require fully converged simulations on at least an 7\u00a0\u00d7\u00a07 grid to cover the full (T,\u00a0P)-range, whereas the local BGs enable efficient evaluation of Gibbs free energies across all thermodynamic states. The resulting phase diagram of Si is presented in Fig.\u00a07<\/a>. It captures the competition between the cubic and \u03b2-tin phases across a wide range of conditions and accurately predicts the corresponding coexistence line. Validation against free energy estimates obtained from MD+MBAR shows excellent agreement, as reflected in the practical indistinguishability of the coexistence lines.<\/p>\n Fig. 7: Silicon phase diagram.<\/b> Si phase diagram for N\u00a0=\u00a0216 particles, computed using local BGs and MD+MBAR. The color scale represents the reduced Gibbs free energy difference per particle between \u03b2-tin and diamond cubic phases, \u0394g\u00a0=\u00a0gcubic\u00a0\u2212\u00a0g\u03b2\u2212tin, as obtained from the local BGs. The red line marks the coexistence line between the two phases, while the black dashed line indicates the corresponding result from MD+MBAR estimates. Error bars for both BG and MD+MBAR, estimated from three independent runs each, are smaller than the plotted line width. Source data are provided as a Source Data file.<\/p>\n","protected":false},"excerpt":{"rendered":"Scaling Boltzmann generators to large systems Central to scaling any architecture in the context of many-body systems is…\n","protected":false},"author":2,"featured_media":521758,"comment_status":"","ping_status":"","sticky":false,"template":"","format":"standard","meta":{"footnotes":"","_share_on_mastodon":"0"},"categories":[271],"tags":[62101,18,1099,19,17,60792,1100,452,133,7016],"class_list":["post-521757","post","type-post","status-publish","format-standard","has-post-thumbnail","category-physics","tag-computational-methods","tag-eire","tag-humanities-and-social-sciences","tag-ie","tag-ireland","tag-method-development","tag-multidisciplinary","tag-physics","tag-science","tag-thermodynamics"],"share_on_mastodon":{"url":"https:\/\/pubeurope.com\/@ie\/116703407889669170","error":""},"_links":{"self":[{"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/posts\/521757","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/comments?post=521757"}],"version-history":[{"count":0,"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/posts\/521757\/revisions"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/media\/521758"}],"wp:attachment":[{"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/media?parent=521757"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/categories?post=521757"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.europesays.com\/ie\/wp-json\/wp\/v2\/tags?post=521757"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}![\"Fig.](\"https:\/\/www.europesays.com\/ie\/wp-content\/uploads\/2026\/06\/41467_2026_73900_Fig1_HTML.png\")
The alternative text for this image may have been generated using AI.<\/p>\n![\"Fig.](\"https:\/\/www.europesays.com\/ie\/wp-content\/uploads\/2026\/06\/41467_2026_73900_Fig2_HTML.png\")
The alternative text for this image may have been generated using AI.<\/p>\n![\"Fig.](\"https:\/\/www.europesays.com\/ie\/wp-content\/uploads\/2026\/06\/41467_2026_73900_Fig3_HTML.png\")
The alternative text for this image may have been generated using AI.<\/p>\n![\"Fig.](\"https:\/\/www.europesays.com\/ie\/wp-content\/uploads\/2026\/06\/41467_2026_73900_Fig4_HTML.png\")
The alternative text for this image may have been generated using AI.<\/p>\n![\"Fig.](\"https:\/\/www.europesays.com\/ie\/wp-content\/uploads\/2026\/06\/41467_2026_73900_Fig5_HTML.png\")
The alternative text for this image may have been generated using AI.<\/p>\n![\"Fig.](\"https:\/\/www.europesays.com\/ie\/wp-content\/uploads\/2026\/06\/41467_2026_73900_Fig6_HTML.png\")
The alternative text for this image may have been generated using AI.<\/p>\n![\"Fig.](\"https:\/\/www.europesays.com\/ie\/wp-content\/uploads\/2026\/06\/41467_2026_73900_Fig7_HTML.png\")
The alternative text for this image may have been generated using AI.<\/p>\n