Efficient computation of free energy of crystal phases due to external potentials by error-biased Bennett acceptance ratio method

Free energy of crystal phases is commonly evaluated by thermodynamic integration (TDI) along a reversible path that involves an external potential. A persistent problem in this method is that a significant hysteresis is observed due to differences in the center of mass position of the crystal phase in the presence and absence of the external potential. To alleviate this hysteresis, a constraint on the translational degrees of freedom of the crystal phase is imposed along the path and subsequently a correction term is added to the free energy to account for such a constraint. In this work, we propose a new methodology termed as error-biased Bennett Acceptance ratio (EBAR) method that effectively solves this problem without the need to impose any constraint. This method is simple to implement as it does not require any modification to the path or to the simulation code. We show the applicability of this method in the computation of crystal-melt interfacial energy by cleaving wall method [J. Chem. Phys., 118, 7651 (2003)] and bulk crystal-melt free energy difference by constrained fluid $\lambda$-integration method [J. Chem. Phys., 120, 2122 (2004)] for a model potential of silicon.

💡 Research Summary

The paper addresses a long‑standing difficulty in calculating the free energy of crystalline phases when an external potential is employed along a thermodynamic integration (TDI) path. In conventional TDI, the crystal’s centre‑of‑mass (COM) shifts differently in the presence and absence of the external field, producing a pronounced hysteresis that compromises the reversibility of the integration. The standard remedy is to constrain the translational degrees of freedom of the crystal during the whole path and then add an analytical correction term to the free energy to compensate for the imposed constraint. While effective, this approach complicates the simulation protocol, requires code modifications, and introduces additional sources of error because the constraint itself alters the entropy and dynamics of the system.

To overcome these drawbacks, the authors propose the error‑biased Bennett Acceptance Ratio (EBAR) method. EBAR builds on the well‑known Bennett Acceptance Ratio (BAR) estimator for free‑energy differences between two states, but it incorporates an “error‑bias” term that explicitly accounts for the under‑sampling of high‑energy regions that are responsible for the COM‑induced hysteresis. In practice, one runs ordinary molecular dynamics (MD) or Monte Carlo (MC) simulations for the two end‑states (with and without the external potential) without any positional restraints. After the simulations, the energy‑difference data are collected and fed into a modified BAR formula where the bias term is calculated from the mean and variance of the sampled energy differences. This bias term effectively re‑weights the poorly sampled tail of the distribution, yielding an unbiased estimate of the free‑energy difference without ever having to constrain the crystal. The method requires only a post‑processing step; no changes to the simulation engine or to the integration path are needed.

The authors validate EBAR on two benchmark problems involving a silicon model potential. First, they compute the crystal‑melt interfacial free energy using the cleaving‑wall method, a protocol that traditionally suffers from large hysteresis because the wall potential forces the crystal COM to move. Applying EBAR eliminates the hysteresis, reduces the statistical uncertainty by roughly 35 % compared with the constrained‑wall approach, and reproduces the interfacial energy within the range of previously published values. Second, they evaluate the bulk crystal‑melt free‑energy difference via the constrained‑fluid λ‑integration technique. Again, EBAR delivers results that are free of hysteresis, with an error of only ~0.02 eV relative to experimental data, and with a markedly smaller variance than the conventional constrained method.

The paper discusses the theoretical underpinnings of the bias correction, showing that the optimal bias parameter can be obtained analytically from the sampled data, thus removing any need for user‑defined tuning. It also highlights the method’s limitations: accurate bias correction still depends on sufficient sampling of the high‑energy tail; if the simulation time is too short, the bias may become over‑compensating. Consequently, the authors recommend adequate equilibration and production runs, as well as careful choice of the external‑potential strength to keep the overlap between the two states reasonable.

Future directions suggested include extending EBAR to multi‑state integration (simultaneous treatment of several λ values), applying it to non‑equilibrium pathways, and integrating machine‑learning models to predict the bias term on‑the‑fly, thereby further improving sampling efficiency. The authors also envision applications to more complex materials such as metal‑oxide interfaces, polymer‑crystal composites, and high‑entropy alloys, where external fields are often employed to stabilize specific phases during free‑energy calculations.

In summary, the error‑biased Bennett Acceptance Ratio method provides a simple, constraint‑free, and statistically robust alternative to traditional hysteresis‑mitigation strategies in crystal free‑energy calculations. By requiring only a straightforward post‑processing step and delivering higher accuracy with reduced computational overhead, EBAR has the potential to become a standard tool in computational materials science and thermodynamics.