Computing Gibbs free energy differences by interface pinning

We propose an approach for computing the Gibbs free energy difference between phases of a material. The method is based on the determination of the average force acting on interfaces that separate the two phases of interest. This force, which depends on the Gibbs free energy difference between the phases, is computed by applying an external harmonic field that couples to a parameter which specifies the two phases. Validated first for the Lennard-Jones model, we demonstrate the flexibility, efficiency, and practical applicability of this approach by computing the melting temperatures of sodium, magnesium, aluminum, and silicon at ambient pressure using density functional theory. Excellent agreement with experiment is found for all four elements, except for silicon, for which the melting temperature is, in agreement with previous simulations, seriously underestimated.

Figure 1

Upper panel: Crystal-liquid configuration from an ab initio simulation of 432 Si atoms in a periodic box. Atoms are colored according to the coordination number (red indicates fourfold coordinated and blue otherwise). Lower panel: Schematic sketch of the Gibbs free energy $G\left(Q\right)$ (black solid line) as a function of the crystallinity order parameter $Q$ at a state point where the liquid is thermodynamically stable and the crystal is metastable. The double arrows indicate the interface contribution $G_i$ (red) and the bulk contribution $N\Delta \mu$ (blue), respectively. The dashed green curve indicates the Gibbs free energy $G^{\prime }\left(Q\right)$ with bias potential applied. The inset shows that the computed $G\left(Q\right)/\left(k_{B}T\right)$ in the two-phase region is indeed linear; here, $G\left(Q\right)$ was computed for the LJ model ($N=5120$) via umbrella sampling (Refs. 2, 26) using Eq. (2) with $\kappa =2$ and a range of $a$'s.

Figure 2

Crystal-liquid coexistence of the LJ model (filled black dots) in the $(p,T)$ plane computed with interface pinning for the LJ model. The solid line is a cubic fit: $-0.5223T^3+5.017T^2+5.502T-5.989$. The computed coexistence line agrees well with results of other methods (Ref. 33): $+$'s and $\times$'s are from Refs. 24 and 34, respectively. The asterisk indicates the gas-liquid critical point ($T_{\mathrm{CP}}=1.31$; $p_{\mathrm{CP}}=0.15$) of the full LJ model (Ref. 35).

Figure 3

(a), (c) $\Delta \mu$ computed with interface-pinning method along an isobar and an isotherm, respectively. (b), (d) Specific entropy $\Delta s$ vs $T$ and specific volume $\Delta v$ vs $p$, respectively; the solid lines in the lower panels are quadratic polynomial fits to these data. The solid lines in the upper panels were computed by integration of these fits. The integration constants were chosen to provide the best overall agreement with the $\Delta \mu$ data.