Direct calculation of the solid-liquid Gibbs free energy difference in a single equilibrium simulation

An important aspect of computational condensed matter physics is to compute phase diagrams of model systems. The naïve approach is to preform a long-time simulation at a selected state point and hope that the system by itself finds its preferred phase, i.e., the phase with the lowest Gibbs free energy. In most cases, this strategy is not viable since first-order transitions are associated with hysteresis. Thus, the system is likely to be stuck in a metastable phase. The origin of this hysteresis effect is the formation of an interface between two phases. The surface tension of the interface gives rise to a free energy barrier that the system has to overcome to transform from one phase to the other.¹ A conceptually appealing and widely used strategy to overcome the hysteresis problem is to preform simulations starting from an initial configuration with two phases in a periodic box.^2–14 When a steady state situation is reached the stable phase will grow at the expense of the other phase. The disadvantages of this approach are that: (i) it relates to a non-equilibrium computation that cannot be done ad infinitum; (ii) a sufficiently large system is needed to reach the steady state growth; (iii) near coexistence, thermal fluctuations will result in some probability that the system will move towards the disfavored phase. In a recent paper,¹⁵ these problems were resolved by introducing the so-called “interface pinning” (IP) method. In short, the idea of this method is to compute the average force needed to keep the system in the two-phase state. This is done by connecting the system to a harmonic field which couples to an order-parameter that discriminates between the two phases of interest. The Gibbs free energy difference between the phases is determined by the average force that the applied field exerts on the system. Thus, a standard ad infinitum equilibrium simulation gives the information needed to compute the Gibbs free energy difference between the two phases of interest.

The purpose of this paper is to give a detailed description of the IP method and show that it is a viable way of computing the solid-liquid Gibbs free energy and construct phase diagrams. As a test case, we investigate the Lennard-Jones (LJ) model.¹⁶ The remainder of the paper is organized as follows. In Sec. II, we describe the IP method in general terms. In Sec. III, we define an order-parameter that distinguishes between solid and liquid by measuring long-range order. In Sec. IV, the IP method is applied to the LJ model. The coexistence line is computed by combining the IP method with the Newton-Raphson method for finding roots, and statistical and systematic errors are investigated. In Sec. V, we compare the IP method to other ways of computing Gibbs free energies and phase diagrams. The paper is completed with a summary.

To introduce the IP method, imagine a two-phase system¹⁷ in a periodic orthorhombic box elongated in the z-direction, i.e., with box lengths X ⩽ Y < Z (Fig. 1). Consider configurations of the Np_zT-ensemble defined as the constant temperature and pressure ensemble where the box lengths X and Y are fixed at values so that the crystal is unstrained, while the box length Z is allowed to fluctuate in order to maintain a constant pressure. Two interfaces will form orthogonally to the long axis. This orientation will minimize the interface area and thereby the interface Gibbs free energy G_i. The system is only barostated in the z-direction since the surface tension will add an a priori unknown pressure to the p_x and p_y components of the pressure tensor.¹⁸ (It is worth noting that an orthorhombic box is not a requirement. The angle between the box vectors X and Y may differ from 90°, but should then be kept constant at an angle that does not strain the crystal.)

Assume that the system is sufficiently large so that the central regions of the pure phase slabs exhibit bulk properties up to some arbitrary threshold values. For a given configuration R = {r₁, r₂, …, r_N} (for simplicity, we assume that particles are point-like) particles can then be labeled either

$s=[ solid ]$⁠,

$l=\textrm {[liquid]}$

$l=[ liquid ]$⁠, or

$i=\textrm {[interface]}$

$i=[ interface ]$ and the total number of particles can be written as N = N_s(R) + N_l(R) + N_i(R). Define

$N_l=\langle N_l(\mathbf {R})\rangle _{N_s}$

$Nl=⟨Nl(R)⟩Ns$ and

$N_i=\langle N_i(\mathbf {R})\rangle _{N_s}$

$Ni=⟨Ni(R)⟩Ns$ to be the average number of liquid and interface particles, respectively, for a given amount of solid particles (i.e., an average over all realizations of capillary waves so that

$N_i=\textrm {const.}$

$Ni= const .$ is independent of N_s)

Let μ_s and μ_l be the chemical potential of the solid and the liquid, respectively. The total Gibbs free energy of the two-phase system along N_s is then

as sketched in Fig. 2. When the relative position of the interfaces change within the two-phase regime, G_i and N_i can be replaced by an constant if “wiggles” can be neglected (wiggles^19,20 are discussed later in the paper). Thus, combining the last two equations gives

where Δμ ≡ μ_s − μ_l and (N − N_i)μ_l + G_i is assumed to be constant (this is discussed and verified later in the paper). Throughout the paper, we let “Δ” denote “

$[ solid ]−[ liquid ]$” and let “

$\textrm {const.}$

$const .$” refer to an constant.

To maintain the system in configurations having two phases, i.e., “pinning the interfaces,” we apply a harmonic field that couples to an order-parameter which relates to the amount of crystal particles in the simulation box. To this aim, we introduce a global order parameter Q(R). Let Q_s and Q_l be the average values of Q(R) when the system is completely solid or liquid, respectively (at a given pressure p and temperature T). We define Q so that it has a linear dependency on the amount of solid particles in a two-phase state when additional degrees of freedom are integrated out (such as phonon vibrations in the slabs of the pure phases)

where

$NiNQi$ is a constant contribution from interface particles. Let U(R) be the energy of the unperturbed system, and

be the energy when a “spring-like” harmonic field is applied. We refer to the field parameters a and κ as the anchor point and the spring constant, respectively. The Gibbs free energy along the Q coordinate when the biasing field applied is (dashed line in Fig. 2)

To give an expression for the probability distribution of Q, we use that P′(Q) ∝ exp (−G′(Q)/k_BT). By insertion of Eq. (3) and elimination of N_s with Eq. (4), we get that

in the two-phase regime. This distribution is shown for the LJ model in Fig. 3.

The chemical potential difference Δμ can be computed from the average force

that the field exerts on the system (along the Q coordinate). When equilibrium is established, the relative position does not change up to thermal fluctuations and F^field = −F^system where

$F system =−∂G∂Q$⁠. By applying the chain rule

$\Delta \mu =\frac{\partial G}{\partial N_s}\break=\frac{\partial G}{\partial Q}\frac{\partial Q}{\partial N_s}$

$Δμ=∂G∂Ns=∂G∂Q∂Q∂Ns$⁠, then

where

$∂Q∂Ns=ΔQN$ is obtained from Eq. (4). Alternatively, a statistical mechanical deviation of this is possible by isolating Δμ from the average of the P′(Q) distribution (Eq. (7)).

Coexistence state points are defined as Δμ(p, T) = 0 and may be computed efficiently using the Newton-Raphson algorithm for finding roots. The required derivatives of Δμ along isobars and isotherms are given by the standard thermodynamic expressions

and

where

In these relations, Δv, Δs, and Δu are changes in specific volume, entropy, and energy, respectively. Thus, coexistence points can be computed by iteration of

or

Iterations are continued until the computed Δμ is zero within the statistical error.

To conclude this section, we give an algorithm for computing coexistence state points in the phase diagram: First, select a pressure p and temperature T for the initial set of simulations. Then,

• construct a crystal configuration in an elongated orthorhombic box;

• determine the lattice constants of the unstrained crystal by performing an NpT simulation in which the box lengths X, Y, and Z are allowed to fluctuate independently to maintain constant pressure;

• compute Q_s and v_s in an Np_zT simulation of the unstrained crystal;

• construct a liquid configuration in an elongated orthorhombic box having the same box lengths X and Y as the unstrained crystal;

• compute Q_l and v_l in an Np_zT simulation of the liquid;

• construct a two-phase configuration having the same box lengths X and Y as the unstrained crystal;

• compute ⟨Q⟩′ in a Np_zT simulation of the two-phase system with an interface pinning

  $\frac{\kappa }{2}(Q(\mathbf {R})-a)^2$
  $κ2(Q(R)−a)2$ field applied;

• calculate Δμ using Eq. (9);

• if Δμ is non-zero within the statistical error, repeat steps i–ix at the pressure given by Eq. (13) or the temperature given by Eq. (14).

We note that an algorithm only involving two-phase simulations can be designed, since a two-phase simulation contains information about the bulk properties of the liquid and the crystal. In this paper, we choose to use the above algorithm for practical reasons.

To utilize the method, we need to define an order parameter Q(R) that distinguishes between the phases of interest. Unlike liquids, crystals have long-ranged translational order and the collective density field may be used to define Q(R)

where

and k = (2πn_x/X, 2πn_y/Y, 0). The integers n_x and n_y are chosen so that the wave vector k corresponds to a Bragg peak. This will maximize the contrast between the liquid and the crystal. The k vector is in the xy-plane (n_z = 0) since Z fluctuates in the Np_zT-ensemble. The factor

$N−12$ ensures scale invariance for the average liquid value, Q_l∝1, while the intensity of a (single) crystal will scale as

$Q_s\propto N^{\frac{1}{2}}$

$Qs∝N12$⁠. For two-phase configurations, Q will scale linearly with the amount of crystal particles (fulfilling Eq. (4)) since the ρ_k-argument of the crystal slab and liquid slab are independent of each other. We note that this order-parameter may be problematic in the supercooled regime, since a crystal can lower |ρ_k| by introducing a long-wavelength displacement of particles. This can be avoided by using a different order-parameter, e.g., the Steinhard Q₆ order-parameter.²³ This was done in Ref. 15. We choose to use |ρ_k| as order-parameter since it is conceptually appealing, simple, and generally applicable.

Equilibrium trajectories can be constructed using standard Monte Carlo sampling or standard molecular dynamics simulations.^21,22 For the latter, forces exerted on particles by the external field have to be evaluated: The force acting on particle j is

where F_j is the force without external field, and

where

$ℜ[ρk]=∑j=1Ncos(k·rj)/N$ and

$\Im [\rho _{\bf k}]\break=-\sum _{j=1}^N\sin ({\bf k}\cdot {\bf r}_j)/\sqrt{N}$

$ℑ[ρk]=−∑j=1Nsin(k·rj)/N$ are the real and imaginary parts of ρ_k, respectively. Forces can be computed with an efficient O(N) scaling algorithm although the force on particle j depends on the position of all of the particles (this typically results in an O(N²) scaling algorithm). This is done in two O(N) steps: (i) compute ρ_k using Eq. (16) and (ii) compute particle forces using Eqs. (17) and (18). Monte Carlo sampling involves evaluation of the energy change δU′ when a particle is moved or the box length Z is changed

where

$δ|ρk|2=|ρk try |2−|ρk current |2$ and

$\delta |\rho _{\bf k}| = |\rho _{\bf k}^{\rm try}|\break-|\rho _{\bf k}^{\rm current}|$

$δ|ρk|=|ρk try |−|ρk current |$⁠. These changes may be computed by evaluating

$\delta \rho _{\bf k}=\rho _{\bf k}^{\rm try}-\rho _{\bf k}^{\rm current}$

$δρk=ρk try −ρk current$ if the current value of

$\rho _{\bf k}=\rho _{\bf k}^{\rm current}$

$ρk=ρk current$ is stored. Moving particle j yields

Thus, computing δρ_k only involves information about particle j allowing for efficient computations. When the box length Z is changed, then δρ_k = 0 and δU′ = δU since k is perpendicular to the z-direction.

As a test case, we apply the IP method to compute solid-liquid Gibbs free energy differences of the LJ model.¹⁶ LJ interactions are truncated and shifted:

$U=∑i>jNuij$ where

$u_{ij}=4\varepsilon ([\frac{\sigma }{r_{ij}}]^{12}-[\frac{\sigma }{r_{ij}}]^{6})-4\varepsilon ([\frac{\sigma }{r_{c}}]^{12}-[\frac{\sigma }{r_{c}}]^{6})$

$uij=4ɛ([σrij]12−[σrij]6)−4ɛ([σrc]12−[σrc]6)$ for r_ij < r_c and zero otherwise. LJ units are used throughout the paper: ε = σ = m = k_B = 1. Two truncation distances are considered: r_c = 2.5 and r_c = 6. The first choice of 2.5 is the standard truncation of the LJ model. The latter choice of 6 is a better approximation of the full LJ model (r_c → ∞). Molecular dynamics simulations are performed using the LAMMPS software package.²⁴ The ρ_k-field was implemented into the package. The Parrinello-Rahman barostat is used²⁵ with a time constant of 8 Lennard-Jones time units together with a Nosé-Hoover^26,27 thermostat with a time constant of τ_NH = 4. Trajectories are evaluated using a time step of 0.004.

As an example, we compute Δμ at p = 1.5 and T = 0.8 (r_c = 2.5) as described in the following. First, a crystal structure of 8×8×20 face centered cubic unit cells (N = 5120) is constructed and simulated for t_sim = 4000. All box lengths are allowed to fluctuate in order to determine the lattice constants of the unstrained crystal giving box lengths of X = Y = 12.92. The unstrained crystal is then simulated for t_sim = 4000 in the Np_zT ensemble, and Q_s = 55.04 (n_x = 16, n_y = 0) and the average partial volume v_s = 1.052 is recorded. Next, a liquid configuration is constructed by melting the crystal in a constant volume simulation by simulating at a high temperature (T = 5). The Np_zT-ensemble (using X = Y = 12.92) of the liquid is then simulated for t_sim = 4000. Q_l = 0.93 and an average specific volume of v_l = 1.177 is recorded. Then, a two-phase configuration is constructed by performing a high temperature constant volume simulation where particles at z < Z/2 are kept at their crystal positions using harmonic springs anchored at crystal sites. The box volume is set in between that of the crystal and the liquid by scaling the box length Z. The Np_zT-ensemble with a harmonic bias-field (a = 27; κ = 10) is then simulated for t_sim = 40 000 to compute ⟨Q⟩′ = 25.246. Equation (9) yields a chemical potential difference of Δμ = 0.080. A configuration from this last simulation is shown in Fig. 1. The two upper panels in Fig. 4 show Δμ along the p = 1.5 isobar and the T = 0.8 isotherm, respectively (r_c = 2.5). The solid lines are computed with thermodynamic integration of Δs and Δv, respectively (shown in the lower panels). The agreement is excellent.

Next, we use the Newton-Raphson method along isotherms to compute coexistence state points. As an example, we computed the T = 0.8 coexistence pressure from the state point at p⁽¹⁾ = 1.5. Equation (13) provides pressures of p⁽ⁱ⁾ = {2.141, 2.189, 2.185(2)}. In the last iteration, the estimated chemical potential difference is zero within the statistical error, Δμ = 0.0001(2) (numbers in parentheses indicate the statistical errors on the last digit). Tables I and II list computed coexistence points using r_c = 6 and r_c = 2.5, respectively. As a consistency check, we note that the computed melting line obeys the Clausius-Clapeyron relation,

$dpmdTm=ΔsΔv$ (two last rows in Tables I and II). The left-hand side of the relation is computed by central differences of the computed melting line.

As an aside, we test a recent theoretical prediction^29,30 related to the melting region in the ρT-plane (Fig. 5): a large class of systems have curves in the phase diagram, referred to as “isomorphs,”²⁹ along which structure, dynamics, and some thermodynamic properties are nearly constant. These are defined as

where T_* is the temperature at a reference state point and

$h(ρ)$³¹ is a function of ρ (not to be confused with the specific enthalpy h). This class of “simple”³² systems is characterized by the property that fluctuations of the virial W (the potential part of the pressure) and the potential energy U are strongly correlated in the NVT ensemble:^33–35 if δW = W − ⟨W⟩ and δU = U − ⟨U⟩, then the correlation coefficient

$R=\langle \delta W\delta U\rangle /\sqrt{\langle (\delta W)^2\rangle \langle (\delta U)^2\rangle }$

$R=⟨δWδU⟩/⟨(δW)2⟩⟨(δU)2⟩$ is close to unity.

For systems with inverse power-law pair interactions, u_ij ∝ r⁻ⁿ, isomorph invariance is trivial³⁶ with

$h(ρ)=(ρ/ρ*)n/3$ where ρ_* is the density at the reference state point. For systems with LJ pair interactions (a sum of two inverse power-laws), the scaling is approximate and reflects an effective inverse power-law of repulsive interactions.³⁷ The apparent exponent is state point and phase dependent^31,38

where γ_* is a constant that may be determined from virial-energy fluctuations in the NVT ensemble at the reference state point by using the identity

$γ=∂log(T)∂log(ρ)Sex=⟨δWδU⟩⟨(δU)2⟩$⁠. Here, it is used that the excess entropy S_ex is isomorph invariant (γ_* = γ evaluated at the reference state point).

The dashed line in Fig. 5 is a liquid isomorph (γ_* = 4.816; T_* = 2; ρ_* = 1/0.9403 = 1.064; R = 0.991) and the solid line is a crystal isomorph (γ_* = 5.517; T_* = 2; ρ_* = 1/0.8827 = 1.133; R = 0.998). At high temperatures and densities, the coexistence region is outlined by these isomorphs. Deviations from the prediction are, however, significant at low ρT. This is properly due to long-range attractive interactions not being accounted for by an effective inverse power-law. Consistent with this interpretation, the melting region of the r_c = 2.5 truncation of pair interactions deviates from the r_c = 6 in this part of phase space. A scaling form of Aρ⁴ + Bρ² (like Eq. (22)) has previously been validated by Khrapak and co-workers.^39,40 They used Rosenfeld's rule of additivity of melting curves⁴¹ to motivate the scaling law.

To estimate the melting line of the full LJ model (r_c → ∞), we apply an approximate pressure correction p^tail that rectifies missing long-range attractions of the truncated model. To this aim, we first consider the pressure correction in a simulation of solid or liquid in bulk: it is convenient to assume that the radial distribution function is constant at distance larger than the truncation. Then the correction is analytic and only depends on ρ and r_c.²¹ Since the densities of the solid and the liquid are different, so are the corrections for the two phases. For the pressure correction of a two-phase simulation, we use the average of the bulk pressure corrections

The third rows in Tables I and II list the corrected melting pressures (p_m + p^tail). Deviations between the corrected melting pressures when truncating at r_c = 2.5 or r_c = 6 are comparable to statistical error.

Computed melting points are shown in Fig. 6 as filled symbols. In the same figure, +'s and ×'s are coexistence points computed with other methods.^28,42 The agreement is excellent. Differences are highlighted in the inset by showing deviations from a cubic fit to the computed melting points. Deviations from the results of Ref. 28 are within statistical error, while the melting pressure of Ref. 42 is systematically too low by about Δp ≃ 0.05 (except for one data point). Systematic errors in computed melting lines are common⁴² and are typically related to approximate tail corrections (like Eq. (23)), finite size effects or method specific systematic errors. In Subsections IV C–IV F, we will discuss systematic and statistical errors related to the IP method.

How does the statistical error of the Δμ estimate depend on the choice of spring constant κ? To answer this question, we compute the Q(t) autocorrelation function (using the Wiener-Khinchin theorem⁴⁵ with fast Fourier transforms)

where δQ(t) = Q(t) − ⟨Q⟩. Fig. 7 shows C(t) for four choices of κ (p = 1.5; T = 0.8; r_c = 2.5). C(t) reveals two relaxation processes that are occurring on different time-scales. We assign them as follows: (i) a fast process related to phonon vibrations and rearrangements of particles, and (ii) a slower over-damped process related to particles moving between phases,

$Ns(t)=−Nl(t)+ const$⁠. For the investigated κ's, the characteristic time τ for the slow process is nearly constant. τ scales as 1/κ for smaller values of κ on this timescale (data not shown). The relative statistical error on the computed Δμ is estimated by dividing runs into 20 independent blocks of length t_block = 100 > τ.⁴⁴ For the shown κ's, spanning two orders of magnitude, the statistical error is independent of κ (inset on Fig. 7).

To investigate finite size effects, we computed Δμ at p = 1.5 and T = 0.8 with r_c = 2.5 using three systems sizes with the same geometry: N = {640, 2160, 5120}, respectively (Fig. 8). For the smallest system sizes, the computed Δμ depends measurably on the relative interface positions (varied by changing a). This effect is not seen for the two larger system sizes. The inset shows the system size dependency of the computed Δμ. Error bars indicate the statistical error plus the variation related to the positions of the interfaces relative to each other. The dashed 1/N line is a guide to the eye. For the largest system sizes, the error of the estimated Δμ is on the order of 10⁻³. This corresponds to an error of the computed melting temperature of about 10⁻³ (Eq. (14)).

We have assumed that the Gibbs free energy in the two-phase region is independent of interface positions relative to each other. There are, however, two effects that may spoil this assumption: (i) if the distance between the interfaces is sufficiently small, particles in one (or both) phases will not have bulk properties, and (ii) “wiggles” on G(N_s).^19,20 To exemplify the latter effect, think of a square lattice gas with attractions between neighboring particles. Fig. 9 shows two two-phase configurations of this model. The interface Gibbs free energy G_i is different for the two configurations resulting in wiggles of G(N_s). The wiggle period on G(Q) is ΔQ/N_z where N_z is the number of crystal planes in the z-direction. To investigate the Gibbs free energy dependency of interface positions of the LJ model, we perform simulations over a range of a's with overlapping P′(Q) distributions. From this, G(Q) is constructed with the umbrella method^21,22 (histogram reweighing is done with the MBAR algorithm⁴⁶). We find no wiggles (Fig. 10), but conjecture that they are hidden in the statistical noise. We emphasize that wiggles may be accounted for and do not constitute a fundamental limitation of the IP method.

How should the anchor point a and the spring constant κ of the harmonic field be chosen to yield the optimal computation? To answer this question, we note that the average distance between the two interfaces should be as large as possible to ensure that pure phase slabs have bulk properties. This distance can be controlled by the anchor point a. To give a guideline for an optimal a, we consider the limit where κ → ∞ or Δμ = 0. From Eq. (9), we find that a = ⟨Q⟩′. The optimal value of a is

The “≃” indicates that the interface contribution Q_i has been ignored. Next, we consider the spring constant κ. There are two arguments for choosing a stiff spring, i.e., a large value of κ: (i) a small κ would not keep the system in two-phase configurations; (ii) the relaxation time of interface dynamics τ scales as 1/κ for small κ's, thus giving bad statistics. There is, however, also an argument for choosing a small κ: Interface fluctuations should span at least one crystal plane to account for wiggles on G(Q). Thus, κ should optimally be chosen so that interface fluctuations span one crystal plane. Setting the standard deviation of the P′(Q) distribution (Eq. (7)) equal to ΔQ/N_z, we get

where N_z is the number of crystal planes in the z-direction. We emphasize that one of the conclusions of this paper is that the IP method is forgiving towards the choice of field parameters.

Let us briefly review other methods for computing Gibbs free energies and phase diagrams^11,21,22 before discussing the advantages and drawbacks of the IP method. In the moving interface approaches, discussed in the Introduction, a simulation is performed of a two-phase system.^2–14 The thermodynamically favored phase will grow in a constant NpT or μVT simulation, allowing to locate coexistence points by changing intensive variables such as p, T, or μ. An alternative is to use an indirect method where the Gibbs free energy of the pure phases is computed in separate simulations. The Gibbs free energy can be computed by Widom insertion^47,48 or by thermodynamic integration to a state of known Gibbs free energy, e.g., an ideal gas,⁴⁹ a harmonic solid,³⁶ or an Einstein solid.⁵⁰ Umbrella sampling or metadynamics along a good reaction coordinate can be used to compute the Gibbs free-energy of transforming the system from one phase to the other.^51–54 The reader is encouraged to explore Refs. 11 and 21, and 22 and references therein for more about methods for computing phase diagrams.

What are the advances and drawbacks of the IP method? First, the IP method inherits the conceptual simplicity, general applicability, and ease to implement the moving interface method. Since the solid-liquid interface is represented explicitly, simulations give information about interface properties. As an example, the interface stiffness may be computed with the capillary fluctuation method.⁵⁵ Crystal growth rates may be computed from Q(t) fluctuations (Fig. 7) similar to the method suggested by Briels and Tepper^56,57 (in the BT method, two-phase configurations are stabilized by keeping the volume constant). We leave such investigations to future publications. A disadvantage of having an explicit interface is that it is in direct contact with the solid and the liquid phases.⁵⁸ In effect, properties of the liquid and the solid slabs may not have bulk values. This is evident in the computed Δμ's of the system size with 640 particles (Fig. 8). Here, Δμ depends systematically on a. This may lead to larger finite size effects compared to methods where free energies are computed in separate simulations having periodic boundaries. Another disadvantage is that a low interface mobility can result in slow dynamics (Fig. 7). This will result in large statistical errors or, in the worst case, it may even be difficult to reach equilibrium.

The IP method constitutes an alternative when other methods are difficult or impossible to use. Widom insertion^47,48 works best for dilute systems whereas it is not a viable option for dense liquids or solids. Thermodynamic integration to a state of known Gibbs free energy^21,22 is problematic when the path includes additional first order transitions. This happens when a phase is surrounded by other phases in the phase diagram. A reference state point can also be nontrivial to identify. Examples are quasi-crystals, liquid-crystals, plastic-crystals, or other phases having a mixture of order and disorder. The IP method is versatile, and may be used to study these phases. Moreover, it can be generalized to be used with two-phase simulations of multi-components systems,⁵⁹ between two fluid phases (gas-liquid or liquid-liquid) or two solid phases. For the latter, one of the crystals will be strained when simulating the two-phase system (assuming that the lattice constants of the crystals are different). This can be rectified by adding the free energy of straining to the computed Δμ. Low mobility of a solid-solid interface could make it unfeasible to use the IP method. Integration along a reaction path of transformation is not trivial, since it is often difficult to identify a suitable coordinate capturing the entire phase transition. As an example, neither the number of crystalline particles N_s nor Q = |ρ_k| are good reaction coordinates, since the cluster of crystalline particles undergoes geometrical transitions^19,60 (in the periodic elongated orthorhombic simulation cell). The shape of the preferred cluster is determined by the surface area, since a small area gives the lowest G_i due to surface tension. For small N_s, a spherical cluster will have the smallest surface area.¹ At some point, when N_s is increased, a cylinder will be the preferred geometry. When N_s ≃ N/2, then a slab has the smallest surface area (Fig. 1). Thus, there are several geometrical free energy barriers orthogonal to the N_s coordinate (with barrier heights that scales with the area of the transition state clusters). With the IP method, the selection of the order parameter Q only has to discriminate between the two phases of interest.

In summary, we have given a detailed description of the IP method and shown that it can be used for efficient calculations of the Gibbs free energy difference between a solid and a liquid. The melting line can be computed efficiently to a high precision when the method is combined with the Newton-Raphson algorithm for finding roots. As an example, the solid-liquid coexistence line of the truncated LJ model line was computed. As an aside, it was shown that the high pressure part of the temperature-density coexistence region is outlined by isomorphs. An approximate pressure correction for rectifying the truncation of pair interactions was given. Statistical errors and systematic variations were investigated.

An important advantage of the IP method is that the solid-liquid Gibbs free energy difference is computed directly in an ad infinitum simulation at a single state point. This makes it versatile and a viable alternative when it is difficult or impossible to perform Widom insertion, thermodynamic integration, or integration along a reaction path of transformation.

This work was financially supported by the Austrian Science Fund FWF within the SFB ViCoM (F41). The author is grateful for valuable comments and suggestions from Christoph Dellago, Gerhard Kahl, Georg Kresse, Felix Hummel, Andreas Tröster, Emanuela Bianchi, Michael Grünwald, David Chandler, Patrick Varilly, David Limmer, Thomas B. Schrøder, Jeppe C. Dyre, and Peter Harrowell.