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.

An accurate location of first order transition lines at a reasonable computational cost is of paramount importance for a wide spectrum of condensed matter systems, ranging from hard to soft materials and biological matter.Basic principles of equilibrium thermodynamics imply that for a given temperature and pressure the system resides in the phase of lowest Gibbs free energy.Phase transitions occur where Gibbs free energy differences between phases vanish, determining phase boundaries in the pressure-temperature plane.From the computational point of view, however, the task of evaluating a phase diagram represents a significant challenge, as phase transitions occur on long time scales [1] such that they cannot be studied using straightforward molecular dynamics or Monte Carlo simulations.
Several numerical approaches have been proposed to cope with this problem [2,3]: (i) in the indirect approach, often based on thermodynamic integration, the Gibbs free energy is computed individually for each of the phases [4][5][6][7][8] and the coexistence line is then calculated by imposing the coexistence condition of equal Gibbs free energy.(ii) Alternatively, in the direct approach, an explicit interface is introduced between the two phases which are then simulated simultaneously in the same simulation box.At fixed pressure and temperature, the system moves towards the phase with the lower Gibbs free energy.Exactly at coexistence the thermodynamic driving force on the interface vanishes and the interface stops moving except for thermal fluctuations.Successful applications of this approach have been reported for a broad spectrum of materials [9][10][11][12][13][14][15][16][17][18][19][20][21].
In this contribution, we present and validate a method to compute the Gibbs free energy difference, ∆G, between two phases.The basic idea of this approach is to compute the average force required to pin the interface of a two-phase system via a harmonic bias potential.This external field couples to a suitably defined order parameter, Q, which distinguishes between the phases of interest.The application of the bias potential effectively transforms the out-of-equilibrium process of the conventional moving interface method into a well-defined equilibrium computation, in which the free energy difference ∆G is determined directly.We refer to this approach as the "interface pinning" method.Coexistence points may subsequently be determined using Newton's root finding method.
To validate our new approach, we have first applied it to the Lennard-Jones (LJ) model [22].Our calculations reproduce with high accuracy the solid-liquid coexistence line identified previously with other approaches [15,23,24] and provide Gibbs free energies that are in excellent agreement with those obtained by thermodynamic integration.We have then used interface pinning in combination with ab initio simulations to compute the melting temperatures of sodium (Na), magnesium (Mg), aluminum (Al) and silicon (Si), demonstrating that this method is efficient, flexible and widely applicable.
Compared to the conventional direct and indirect methods used in the literature so far, interface pinning offers several advantages.In contrast to the direct approaches, interface pinning operates at well-defined equilibrium conditions, thus permitting the explicit calculation of free energy differences and interface properties.The selection of the order parameter Q does not need to be a reaction coordinate capturing the entire transformation mechanism.Finally, interface pinning inherits the general applicability and conceptual simplicity of the direct approaches.The latter makes it easy to implementation into existing programs.
To introduce the method, consider a two phase crystalliquid system [25] in a periodic orthorhombic box (see figure 1) at temperature T and pressure p.The box lengths X and Y are kept constant at values for which the crystal is unstrained, while the box length Z is allowed to change in order to maintain constant pressure.We refer to this ensemble as the N p z T -ensemble.To lower the interface Gibbs free energy G i , the system will have two interfaces in the XY -plane minimizing the interface surface area.We assume that the system is large enough to represent bulk phases at least at the center of the liquid and solid slabs.Particles may then either be labeled crystalline (subscript c), liquid (subscript l) or interfacial (subscript i), so that the total number of particles is N = N c + N l + N i .The contributions to the total Gibbs free energy of particles in the bulk phases is determined by the chemical potentials µ c and µ l of the solid and liquid, respectively, and the total Gibbs free energy is When the relative distance between the interfaces changes, particles are transferred between the bulk phases.Assuming that the interface quantities G i and N i do not change when the interfaces shift due to the growth of one bulk phase at the cost of the other, the number of liquid particles may be written as and the Gibbs free energy is given by where ∆µ ≡ µ c − µ l .Throughout the paper we will use the subscripts c and l for crystal and liquid properties, respectively, and let "∆" denote "[crystal] − [liquid]".
To sample configurations in the two-phase region and to prevent the system from complete transformation into one of the pure phases, we apply a harmonic bias potential that pins the relative position of the interfaces.Let U (R) be the energy of the unbiased system for configuration R = {r 1 , r 2 , . . ., r N }, and the energy of the system plus the bias potential.Here, Q(R) is a global order parameter with a linear dependence on the number of particles in the solid phase: In the biased system, the position of the interfaces relative to each other will fluctuate around an average value and the order parameter Q will fluctuate accordingly.The probability distribution of is the Gibbs free energy along the Q coordinate of the biased system, and Z is the corresponding partition function.The Gibbs free energy of the biased system may be written in terms of the unbiased free energy as By insertion of Equs.( 1) and (3), it follows that P (Q) is Gaussian, where α = N ∆µ/∆Q is the slope of G(Q) in the twophase region, displayed in the lower panel of (2) with κ = 2 and a range of a's.
distribution P (Q) has variance σ 2 Q = k B T /κ and average Q = a − α/κ, and the chemical potential difference between the two phases may be computed as As a guideline, we choose κ such that typical fluctuations in Q correspond to one or a fraction of a crystal plane, and a such that the system is approximately half liquid and half crystal.In practice, we find that a wide range of field parameters give the same precision of the ∆µ estimate [27].
To apply the interface pinning method in practice, we must choose an order parameter Q that grows linearly with the number of crystalline particles N c in the twophase region.Moreover, Q should be computationally inexpensive.Unlike liquids, crystals have long-ranged translational order, allowing us to use the collective density field as order parameter: Here, k = (2πn x /X, 2πn y /Y, 0) for some fixed integers (n x , n y ) that should be chosen such that k correspond to a Bragg peak.This choice will maximize the contrast between the liquid and the crystal.
. Derivatives of Q with respect to the particle coordinates, required to determine the forces resulting from the bias, can be computed with an algorithm scaling as O(N ).We note that this order parameter may be problematic in the supercooled regime, since a crystal can lower |ρ k | by introducing long wave length displacements of particles.The energy penalty of such displacements is low and decreases with increasing system size.We have chosen to use |ρ k | as order parameter for most computations, since it is generally applicable and simple.For some computations we have in addition used the Steinhardt Q = Q 6 order parameter [28], which has the advantages of being robust in the supercooled regime.The two choices of order parameter give the same ∆µ's within statistical error.A more detailed description of the method will be given in a forthcoming publication [27].
To verify the method, we first used it to determine the solid-liquid coexistence line of the LJ model with truncated pair interactions: U (R) = N i>j u(r ij ), where u(r) = 4(r −12 − r −6 ) − 4(6 −12 − 6 −6 ) for r < 6 and zero otherwise (LJ units are used for this model throughout the paper).MD simulations with a time step of t step = 0.004 were performed using the LAMMPS software package [32] modified to include the bias potential.The Parrinello-Rahman barostat was used [33] with a time constant of τ PR = 8 together with a Nosé-Hoover [34,35] thermostat with a time constant of τ NH = 4.
Results presented in Fig. 2 demonstrate that the solidliquid coexistence line for the LJ model computed by interface pinning agrees within high precision with data obtained using other methods [24,30].The coexistence points displayed in Fig. 2 were computed as follows.First, a crystal structure of 8×8×20 face centered cubic (fcc) unit cells (N = 5120) was constructed and simulated at p = 1 and T = 0.8 for t sim = 800.All box lengths were allowed to fluctuate in order to determine the geometry of the unstrained crystal, giving X = Y = 12.85.The unstrained crystal was then simulated for t sim = 800  [24] and [30], respectively.The asterisk indicates the gas-liquid critical point (T CP = 1.31; p CP = 0.15) of the full LJ model [31].
in the N p z T ensemble, and Q c = 56.31(n x = 16, n y = 0) and the average partial volume v c = 1.036 was recorded.Next, a liquid was prepared by melting the crystal in a constant volume simulation at high temperature (T = 5).
The N p z T -ensemble (using X = Y = 12.85) of the liquid was simulated for t sim = 800, and Q l = 0.94 and the average specific volume of the liquid v l = 1.163 was recorded.Then, a two phase configuration was constructed by performing a high temperature constant volume simulation where particles at z < Z/2 were kept at their crystal positions using harmonic springs anchored at crystal sites, with the box volume (length Z) in between that of the crystal and the liquid.The N p z T -ensemble with the bias-field of Equ. ( 2) with parameters a = 26 and κ = 4 was simulated for t sim = 4000 to compute Q = 25.055.Application of Equ. ( 5) yielded a chemical potential difference of ∆µ = 0.040.The coexistence pressure was then determined iteratively using Newton's root finding method along the isotherm: p (i+1) = p (i) − ∆µ (i) /∆v (i) , providing pressures of p (i) = {1.0,1.320, 1.337(1)}.In the last iteration the estimated chemical potential difference is zero within the error bar, ∆µ = −0.0007(10) (throughout the paper numbers in parentheses indicate the statistical errors of the last digits).In Fig. 3 we confirm that ∆µ(p, T ) computed with interface pinning (symbols) is consistent with thermodynamic integration (lines).Due to its efficiency and flexibility, the interface pinning approach can be combined with electronic structure methods and ab initio molecular dynamics to com- p = 1 T = 0.8

FIG. 3. (color online)
. Upper panels (a and c): ∆µ computed with interface pinning method along an isobar and an isotherm, respectively.Lower panels (b and d): specific entropy ∆s vs. T and specific volume ∆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 ∆µ-data.puted free energy differences from first principles.For the present simulations, the method was implemented in the Vienna ab inito Simulation Package [38].As an example, we used interface pinning to compute the melting temperatures T m of the period three elements Na, Mg, Al and Si at ambient pressure.Computed T m 's are shown in Table I along with simulation details.Melting temperatures were computed first for crystalline Si in the fourfold coordinated cubic diamond (cd) structure (see Fig. 1).To be compatible to previous calculations [8,39], density functional theory (DFT) in the local density approximation (LDA) within the framework of the projector augmented wave method was used [40].N pT and N p z T simulations were performed using a time step of t step = 3 fs with a Parrinello-Rahman barostat [33] and a Langevin thermostat [41].To compute the Si coexis-tence temperature at ambient pressure, we use a similar strategy as outlined for the LJ model: bulk properties of the crystal and the liquid (Q c , Q l , v c , v l , X and Y ) were evaluated in simulations for 216 Si atoms (3×3×3 conventional cells; t sim = 60 ps) at T = 1200 K. Next, solid-liquid simulations with a bias field (t sim > 30 ps) were performed for four system sizes: {2×2×4, 2×2×7, 3×3×6, 3×4×7} conventional cubic cells corresponding to N = {128, 224, 432, 672} atoms.Coexistence temperatures were estimated to be {1189,1218,1225,1241} K using T m T + ∆µ ∆s .Finally, finite size effects were extrapolated assuming a 1/N decay of the finite size error yielding T m = 1250(10) K.The finite size effects are particularly large for liquid silicon, since the metallic liquid is embedded in a semiconducting host, resulting in a discretization of the electronic states in the metal (electron in a box).The present value is fully consistent with previous LDA calculations [8,39], and the discrepancy to the experimental value of T m = 1635 K originates from an underestimation of the energy difference between four fold coordinated semiconducting Si in the cd structure and six fold coordinated metallic Si in liquid Si resembling the β-tin structure [8].The entropy of fusion ∆s(T m ) = 3.5(1) k B /atom and the slope of the melting curve d T m /d p = −51(7) K/GPa (computed using the Clausius-Clapeyron relation) are also in agreement with previous theoretical results [8,39].For the other elements, Na, Mg and Al, finite size effects are less critical, and we only considered system sizes comparable to the largest Si system.For these three elements, the calculations were performed using PBEsol (Perdew, Burke, Ernzerhof functional for solids) [42], which yields more accurate lattice constants than the LDA.The computed T m 's of Na and Mg are in excellent agreement with experimental values, while for Al the computed T m is about 6% too large (see Table I).
In summary, we have introduced a computational method that allows a direct evaluation of the Gibbs free energy differences between two phases.In contrast to previous approaches, simulations are carried out at equilibrium conditions by pinning the interface between the phases of interest via a harmonic bias potential that couples to a suitably defined order parameter.Application of interface pinning to the LJ model demonstrates the accuracy and efficiency of this new approach: the solidliquid coexistence line agrees to high accuracy with data obtained by other methods and the computed Gibbs free energy is consistent with data obtained via thermodynamic integration.The practical applicability and flexibility of the method was demonstrated by computing the melting points of Na, Mg, Al and Si at ambient pressure using first principles simulations.The results demonstrate that present density functionals yield very accurate melting temperatures for crystalline metal to liquid metal transitions, but errors are sizable for semiconductor (cd-Si) to metal transitions (liquid Si).Furthermore, the present approach allows to compute directly and straightforwardly structural and thermodynamic properties of the interface, such as surface tension from the pressure tensor [43] or the crystal growth rate from Q(t) fluctuations [44,45].
This work was financially supported by the Austrian Science Fund FWF within the SFB ViCoM (F41).Supercomputing time on the Vienna Scientific cluster (VSC) is gratefully acknowledged.

Fig. 1 FIG. 1 .
FIG. 1. (color online).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=[fourfold coordinated] and blue otherwise).Lower panel: schematic sketch of the Gibbs free energy G(Q) (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 Gi (red) and the bulk contribution N ∆µ (blue), respectively.The dashed green curve indicates the Gibbs free energy G (Q) with bias potential applied.The inset shows that the computed G(Q)/(kBT ) in the two-phase region is indeed linear; here, G(Q) was computed for the LJ model (N = 5120) via umbrella sampling[2,26] using Equ.(2) with κ = 2 and a range of a's.

TABLE I .
Ab inito and experimental (Refs.36 and 37) melting temperatures Tm (in K) of period 3 elements using either |ρ k | or Q6 as order parameter."Super cell" indicates the applied super cell built from the conventional cubic cell (including liquid and solid part).N is the total particle number.