Computing phase diagrams of model systems is an essential part of computational condensed matter physics. In this paper, we discuss in detail the interface pinning (IP) method for calculation of the Gibbs free energy difference between a solid and a liquid. This is done in a single equilibrium simulation by applying a harmonic field that biases the system towards two-phase configurations. The Gibbs free energy difference between the phases is determined from the average force that the applied field exerts on the system. As a test system, we study the Lennard-Jones model. It is shown that the coexistence line can be computed efficiently to a high precision when the IP method is combined with the Newton-Raphson method for finding roots. Statistical and systematic errors are investigated. Advantages and drawbacks of the IP method are discussed. The high pressure part of the temperature-density coexistence region is outlined by isomorphs.

## I. INTRODUCTION

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.^{1} 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,^{15} 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.^{16} 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.

## II. THE INTERFACE PINNING METHOD

To introduce the IP method, imagine a two-phase system^{17} 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*_{z}*T*-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.^{18} (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**_{1}, **r**_{2}, …, **r**_{N}} (for simplicity, we assume that particles are point-like) particles can then be labeled either

*N*=

*N*

_{s}(

**R**) +

*N*

_{l}(

**R**) +

*N*

_{i}(

**R**). Define

*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 “

### A. Harmonic field biasing towards two-phase configurations

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

*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*_{B}*T*). 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.

### B. Computing Δμ from the average force exerted by the applied field on the system

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

### C. Computing coexistence state points with the Newton-Raphson method for finding roots

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.

### D. Algorithm for computing coexistence state points

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*_{z}*T*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*_{z}*T*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*_{z}*T*simulation of the two-phase system with an interface pinning$\frac{\kappa }{2}(Q(\mathbf {R})-a)^2$$\kappa 2(Q(R)\u2212a)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.

## III. TRANSLATIONAL ORDER PARAMETER

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*_{z}*T*-ensemble. The factor

*Q*

_{l}∝1, while the intensity of a (single) crystal will scale as

*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*

_{6}order-parameter.

^{23}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}, 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*

^{2}) 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

*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.

## IV. SOLID-LIQUID COMPUTATIONS OF THE LENNARD-JONES MODEL

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

*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.

^{24}The ρ

_{k}-field was implemented into the package. The Parrinello-Rahman barostat is used

^{25}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*_{z}*T* 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*_{z}*T*-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*_{z}*T*-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)} = 1.5. Equation (13) provides pressures of *p*^{(i)} = {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,

T _{ m } | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 | 1.2 | 1.4 | 1.6 | 1.8 | 2.0 | 2.2 | 2.4 | 2.6 |

p _{ m } | −0.212 | 0.928 | 2.185 | 3.514 | 4.939 | 7.921 | 11.181 | 14.632 | 18.180 | 22.007 | 26.029 | 30.050 | 34.314 |

p_{m} + p^{tail} | −1.046 | 0.049 | 1.264 | 2.555 | 3.943 | 6.859 | 10.056 | 13.448 | 16.943 | 20.717 | 24.688 | 28.661 | 32.878 |

v _{ s } | 1.0614 | 1.0452 | 1.0277 | 1.0110 | 0.9951 | 0.9672 | 0.9421 | 0.9202 | 0.9014 | 0.8835 | 0.8671 | 0.8530 | 0.8394 |

v _{ l } | 1.2194 | 1.1714 | 1.1360 | 1.1080 | 1.0830 | 1.0446 | 1.0117 | 0.9838 | 0.9612 | 0.9399 | 0.9211 | 0.9043 | 0.8888 |

u _{ s } | −5.358 | −5.156 | −4.953 | −4.742 | −4.513 | −4.020 | −3.483 | −2.907 | −2.301 | −1.663 | −0.997 | −0.315 | 0.394 |

u _{ l } | −4.294 | −4.218 | −4.075 | −3.888 | −3.683 | −3.183 | −2.627 | −2.041 | −1.400 | −0.727 | −0.009 | 0.696 | 1.446 |

Δs | −1.718 | −1.507 | −1.392 | −1.327 | −1.263 | −1.207 | −1.168 | −1.123 | −1.107 | −1.091 | −1.087 | −1.063 | −1.055 |

$\frac{dp_m}{dT_m}$ $ d p m d T m $ | … | 12.0^{a} | 12.9 | 13.8 | 14.7 | 15.6 | 16.8 | 17.5 | 18.4 | 19.6 | 20.1 | 20.7 | … |

$\frac{\Delta s}{\Delta v}$ $ \Delta s \Delta v $ | 10.9 | 11.9 | 12.9 | 13.7 | 14.4 | 15.6 | 16.8 | 17.6 | 18.5 | 19.3 | 20.1 | 20.7 | 21.4 |

T _{ m } | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 | 1.2 | 1.4 | 1.6 | 1.8 | 2.0 | 2.2 | 2.4 | 2.6 |

p _{ m } | −0.212 | 0.928 | 2.185 | 3.514 | 4.939 | 7.921 | 11.181 | 14.632 | 18.180 | 22.007 | 26.029 | 30.050 | 34.314 |

p_{m} + p^{tail} | −1.046 | 0.049 | 1.264 | 2.555 | 3.943 | 6.859 | 10.056 | 13.448 | 16.943 | 20.717 | 24.688 | 28.661 | 32.878 |

v _{ s } | 1.0614 | 1.0452 | 1.0277 | 1.0110 | 0.9951 | 0.9672 | 0.9421 | 0.9202 | 0.9014 | 0.8835 | 0.8671 | 0.8530 | 0.8394 |

v _{ l } | 1.2194 | 1.1714 | 1.1360 | 1.1080 | 1.0830 | 1.0446 | 1.0117 | 0.9838 | 0.9612 | 0.9399 | 0.9211 | 0.9043 | 0.8888 |

u _{ s } | −5.358 | −5.156 | −4.953 | −4.742 | −4.513 | −4.020 | −3.483 | −2.907 | −2.301 | −1.663 | −0.997 | −0.315 | 0.394 |

u _{ l } | −4.294 | −4.218 | −4.075 | −3.888 | −3.683 | −3.183 | −2.627 | −2.041 | −1.400 | −0.727 | −0.009 | 0.696 | 1.446 |

Δs | −1.718 | −1.507 | −1.392 | −1.327 | −1.263 | −1.207 | −1.168 | −1.123 | −1.107 | −1.091 | −1.087 | −1.063 | −1.055 |

$\frac{dp_m}{dT_m}$ $ d p m d T m $ | … | 12.0^{a} | 12.9 | 13.8 | 14.7 | 15.6 | 16.8 | 17.5 | 18.4 | 19.6 | 20.1 | 20.7 | … |

$\frac{\Delta s}{\Delta v}$ $ \Delta s \Delta v $ | 10.9 | 11.9 | 12.9 | 13.7 | 14.4 | 15.6 | 16.8 | 17.6 | 18.5 | 19.3 | 20.1 | 20.7 | 21.4 |

^{a}

|$\frac{dp_m}{dT_m}$|$dpmdTm$ computed by central difference of values in the first two rows.

T _{ m } | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 | 1.2 | 1.4 | 1.6 | 1.8 | 2.0 | 2.2 | 2.4 | 2.6 |

p _{ m } | −0.970 | 0.132 | 1.337 | 2.629 | 4.012 | 6.930 | 10.145 | 13.549 | 17.104 | 20.857 | 24.850 | 28.916 | 33.041 |

p_{m} + p^{tail} | −1.030 | 0.068 | 1.270 | 2.560 | 3.940 | 6.853 | 10.063 | 13.463 | 17.014 | 20.763 | 24.753 | 28.815 | 32.937 |

v _{ s } | 1.0553 | 1.0400 | 1.0242 | 1.0086 | 0.9933 | 0.9661 | 0.9412 | 0.9194 | 0.9002 | 0.8827 | 0.8663 | 0.8518 | 0.8388 |

v _{ l } | 1.2358 | 1.1804 | 1.1425 | 1.1120 | 1.0864 | 1.0470 | 1.0127 | 0.9852 | 0.9616 | 0.9403 | 0.9208 | 0.9038 | 0.8885 |

u _{ s } | −6.314 | −6.125 | −5.929 | −5.722 | −5.506 | −5.037 | −4.526 | −3.972 | −3.385 | −2.768 | −2.120 | −1.451 | −0.764 |

u _{ l } | −5.024 | −5.008 | −4.902 | −4.755 | −4.570 | −4.106 | −3.603 | −3.021 | −2.409 | −1.762 | −1.085 | −0.379 | 0.325 |

Δs | −1.858 | −1.622 | −1.480 | −1.377 | −1.311 | −1.245 | −1.177 | −1.153 | −1.125 | −1.103 | −1.085 | −1.073 | −1.050 |

$\frac{dp_m}{dT_m}$ $ d p m d T m $ | … | 11.5^{a} | 12.5 | 13.4 | 14.3 | 15.3 | 16.5 | 17.4 | 18.3 | 19.4 | 20.1 | 20.5 | … |

$\frac{\Delta s}{\Delta v}$ $ \Delta s \Delta v $ | 10.3 | 11.6 | 12.5 | 13.3 | 14.1 | 15.4 | 16.5 | 17.5 | 18.3 | 19.1 | 20.0 | 20.6 | 21.1 |

T _{ m } | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 | 1.2 | 1.4 | 1.6 | 1.8 | 2.0 | 2.2 | 2.4 | 2.6 |

p _{ m } | −0.970 | 0.132 | 1.337 | 2.629 | 4.012 | 6.930 | 10.145 | 13.549 | 17.104 | 20.857 | 24.850 | 28.916 | 33.041 |

p_{m} + p^{tail} | −1.030 | 0.068 | 1.270 | 2.560 | 3.940 | 6.853 | 10.063 | 13.463 | 17.014 | 20.763 | 24.753 | 28.815 | 32.937 |

v _{ s } | 1.0553 | 1.0400 | 1.0242 | 1.0086 | 0.9933 | 0.9661 | 0.9412 | 0.9194 | 0.9002 | 0.8827 | 0.8663 | 0.8518 | 0.8388 |

v _{ l } | 1.2358 | 1.1804 | 1.1425 | 1.1120 | 1.0864 | 1.0470 | 1.0127 | 0.9852 | 0.9616 | 0.9403 | 0.9208 | 0.9038 | 0.8885 |

u _{ s } | −6.314 | −6.125 | −5.929 | −5.722 | −5.506 | −5.037 | −4.526 | −3.972 | −3.385 | −2.768 | −2.120 | −1.451 | −0.764 |

u _{ l } | −5.024 | −5.008 | −4.902 | −4.755 | −4.570 | −4.106 | −3.603 | −3.021 | −2.409 | −1.762 | −1.085 | −0.379 | 0.325 |

Δs | −1.858 | −1.622 | −1.480 | −1.377 | −1.311 | −1.245 | −1.177 | −1.153 | −1.125 | −1.103 | −1.085 | −1.073 | −1.050 |

$\frac{dp_m}{dT_m}$ $ d p m d T m $ | … | 11.5^{a} | 12.5 | 13.4 | 14.3 | 15.3 | 16.5 | 17.4 | 18.3 | 19.4 | 20.1 | 20.5 | … |

$\frac{\Delta s}{\Delta v}$ $ \Delta s \Delta v $ | 10.3 | 11.6 | 12.5 | 13.3 | 14.1 | 15.4 | 16.5 | 17.5 | 18.3 | 19.1 | 20.0 | 20.6 | 21.1 |

^{a}

|$\frac{dp_m}{dT_m}$|$dpmdTm$ computed by central difference of values in the first two rows.

### A. Isomorph prediction of the ρ*T* coexistence region

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,”^{29} 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

^{31}is a function of ρ (not to be confused with the specific enthalpy

*h*). This class of “simple”

^{32}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

For systems with inverse power-law pair interactions, *u*_{ij} ∝ *r*^{−n}, isomorph invariance is trivial^{36} with

_{*}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.

^{37}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

*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*ρ^{4} + *B*ρ^{2} (like Eq. (22)) has previously been validated by Khrapak and co-workers.^{39,40} They used Rosenfeld's rule of additivity of melting curves^{41} to motivate the scaling law.

### B. Correcting for missing long-range attractions

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}.^{21} 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^{42} 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.

### C. Statistical error

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^{45} 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,

*t*

_{block}= 100 > τ.

^{44}For the shown κ's, spanning two orders of magnitude, the statistical error is independent of κ (inset on Fig. 7).

### D. Systematic errors at small system sizes

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^{−3}. This corresponds to an error of the computed melting temperature of about 10^{−3} (Eq. (14)).

### E. Gibbs free energy dependency of the interface positions relative to each other

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^{46}). 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.

### F. Guidelines for choosing *a* and κ

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.

## V. ADVANTAGES AND DRAWBACKS OF THE INTERFACE PINNING METHOD

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,^{49} a harmonic solid,^{36} or an Einstein solid.^{50} 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.^{55} 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.^{58} 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,^{59} 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.^{1} 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.

## VI. SUMMARY

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.

## ACKNOWLEDGMENTS

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.

## REFERENCES

*Np*

_{z}

*T*ensemble must be replaced by another ensemble (since it is only defined for single component systems). As an example, for a two component system one may consider the

*N*

_{A}

*N*

_{B}

*p*

_{z}

*T*ensemble or the μ

_{A}μ

_{B}

*VT*ensemble. We leave such investigations for a future study.