The precise estimation of the location of phase transitions is an essential task in the study of many condensed matter systems. A recently developed technique denoted interface pinning (IP) [U. R. Pedersen, F. Hummel, G. Kresse, G. Kahl, and C. Dellago, Phys. Rev. B.88, 094101 (2013); U. R. Pedersen, J. Chem. Phys.139, 104102 (2013)] can accurately estimate the location of fluid-solid transition using the NPzT ensemble for single-component systems by computing the free energy difference between a solid and a fluid. The IP method is extended here to be applicable to different ensembles for both single-component systems and binary mixtures. A more general scheme is also proposed for the extrapolation of properties targeting coexistence conditions. This framework is used to estimate the coexistence pressure for the isotropic-rotator phase transition for three single-component polyhedral systems and to estimate isotropic-crystal coexistence compositions for a binary mixture of hard cubes and spheres. In addition, by exploring various choices for the order parameter used to distinguish between the isotropic and ordered phases, it is found that volume provides a reasonable alternative to translational order parameters which can be either more expensive to calculate or unable to pin a two-phase interfacial state.

At present, a major focus in material science is the self-assembly of anisotropic particles due to its potential to help engineer novel materials from colloidal nanoparticles.1,2 Assemblies of anisotropic particles often undergo order-disorder phase transitions involving changes in both translational and orientational degrees of freedom and can lead to phases with partial structural order or “mesophases.”3–7 The accurate location of these order-disorder phase transitions is of utmost importance to better understand not only the thermodynamic behavior of these systems but also the kinetics of these transitions. Several numerical approaches8 have been developed to pinpoint phase transitions and are mainly classified into 2 groups popularly known as direct and indirect approaches. The indirect approach involves the computation of free energies of pure phases separately via thermodynamic integration,9–12 and locating coexistence lines by imposing the condition of equal free energies of both phases at equilibrium. Though the indirect approach has proven to be quite accurate, the calculation of the coexistence points typically requires numerous simulations, making it at times less appealing than the more straightforward direct approach. In the direct approach, devised by Ladd and Woodock,13 simulations with explicit interfaces between two phases are performed to obtain the conditions at which the system is at coexistence. The direct approach has been found to be successful in application to a broad spectrum of systems including hard spheres,14–17 Lennard Jones spheres,18 hard anisotropic particles,4 mixtures,19,20 and molecular models such as water.21–28 

A recent variant of the direct approach is the Interface Pinning (IP) method proposed by Pedersen et al.29,30 The IP method is easy to implement, conceptually simple, and unlike other direct approach methods, it computes the free energy difference between two phases by performing an ad infinitum equilibrium simulation. The IP method is based on the NPzT simulations of coexisting phases to determine the average force acting on the interface. The force is estimated by applying a harmonic field coupled with the order parameter that biases the system towards a two-phase configuration. The theory and the implementation of the method are described in detail by Pedersen.30 

One key aspect of the IP method is the use of an appropriate order parameter capable of discriminating between the two different phases present in the simulation box. The IP operates by “pinning” or stabilizing the two-phase interfacial state of the system via a biasing harmonic potential that forces the system to approach a global order parameter that is intermediate between those of the two single-phase states (at the same external conditions). Pedersen et al.29,30 consider solid-liquid transitions and hence the types of order parameters adopted were such that they can detect and quantify the extent of translational order in the system. While translational order parameters would appear to be the most informative for such situations, no consideration was given to other “blunter” but still phase discriminating order parameters. In particular, the specific volume (or its reciprocal “density”) was not tested as order parameter, even though some versions of interfacial simulation operate by fixing the total specific volume of the system to a value that is intermediate between those for the two single-phase states.4,6,19,20 In a sense, one could see such an implementation as using the IP method with specific volume as order parameter and pinning it to a specific value with an infinite spring constant.

In this article, we first generalize the IP method to make it readily implementable for different ensembles including mixture systems. We also adopt a second-order scheme for the extrapolation of free energies to target coexistence conditions; this has the potential advantage of reducing the number of iterations required to pinpoint coexistence conditions when compared to the first-order scheme used in Pedersen et al.29 We then apply the IP method to four anisotropic-particle systems, (1) hard cuboctahedra (CO), (2) truncated octahedra (TO), (3) rhombic dodecahedra (RD), and (4) a binary mixture of cubes and spheres. We chose these systems as their liquid-solid phase behavior is already known;4,6,7 the fourth system provides an illustration of a mixture system while the first three are examples for the formation of an interesting plastic solid phase called rotator phase (whose disorder-to-order kinetics have also been recently studied31). For the pure-component systems, we also compare the performance of translational order parameters and specific volume in their ability to estimate coexistence conditions via IP, showing the suitability of the latter. In Secs. II–IV we describe the methods used, the main results obtained, and the conclusions reached.

In Sec. II A, a generalized notation is first introduced to describe typical statistical mechanical ensembles, which is then used in Sec. II B to outline the basis of the IP method for an arbitrary ensemble. In Sec. II C, an improved extrapolation technique is presented to estimate coexistence conditions. In Sec. II D, it is shown how key generalized formulas developed in Secs. II A–II C are specialized for the particular ensembles to be used in the forthcoming simulations. In Sec. II E, several order parameters suitable for studying solid-fluid phase equilibrium are described; in particular, the definitions of both the conventional translational order parameters used in the original IP method and the specific volume as an alternative order parameter. Finally, the choice of the pinning input parameters is explained in Sec. II F and the main steps in the IP method are summarized in Sec. II G.

We follow Ref. 32 to adopt a general notation to allow application of the IP method to different ensembles. Consider a system where the ensemble properties M, f1, and f2 are fixed, where M is an extensive property and the f's are thermodynamic fields. The probability of any microstate R in that system can be written as

\begin{equation}p({\bf R}|f_1,f_2) = \exp( { - f_1 X_1 - f_2 X_2 })/Q(f_1,f_2),\end{equation}
p(R|f1,f2)=exp(f1X1f2X2)/Q(f1,f2),
(1)

where Xi is the extensive property conjugate to field fi, and Q is the partition function which is associated with a dimensionless free energy Φ (with ϕ being the intensive property) via

\begin{equation}\Phi = \phi M = - {\rm ln}\,Q(f_1,f_2).\end{equation}
Φ=ϕM= ln Q(f1,f2).
(2)

The probability of observing the system with order parameter λ (R) = λi can be written as

\begin{equation}p(\lambda _i) \propto \exp( { - \Phi _c (\lambda _i)} ),\end{equation}
p(λi)exp(Φc(λi)),
(3)

where Φc(λi) is the free energy function containing only the contribution of microstates with order parameter λi,

\begin{equation}\Phi _c (\lambda _i) = - {\rm ln}\left[ {\sum\limits_{\bf R} {\delta [\lambda ({\bf R}) - \lambda _i]\exp\left( { - f_1 X_1 - f_2 X_2 } \right)} } \right].\end{equation}
Φc(λi)= ln Rδ[λ(R)λi]expf1X1f2X2.
(4)

The fundamental thermodynamic equation for this Mf1f2 ensemble is

\begin{equation}d\phi = \bar x_1 df_1 + \bar x_2 df_2,\end{equation}
dϕ=x¯1df1+x¯2df2,
(5)

where

\begin{equation}\bar x_i = \bar X_i /M = ( {\partial \phi /\partial f_i } )_{f_j,j \ne i}.\end{equation}
x¯i=X¯i/M=(ϕ/fi)fj,ji.
(6)

This formulation can be easily generalized to cases where more variables enter into the ensemble formulation. For instance, cases with more than 2 field properties can be treated identically as long as the value of only one field property is unknown at the coexistence conditions (the others being specified by the user). An example with 3 field properties is given in Sec. II D.

The IP method by Pedersen et al. was originally applied in the NPzT ensemble to estimate the difference in the chemical potential between a solid and a liquid present in an interfacial setup. In this section, we recapitulate the basis of the IP method as pertaining to the generalized ensemble described in Sec. II A. For an Mf1f2 ensemble, if two phases I and II (used therefrom as superscripts) are at coexistence, then it must be true that not only |$f_1^I = f_1^{II}$|f1I=f1II and |$f_2^I = f_2^{II}$|f2I=f2II but also that ϕI = ϕII. In an interfacial simulation at fixed M (total), f1 and f2, two phases with bulk-like quality are present in quantities MI and MII, but if they are not at equilibrium then ϕI ≠ ϕII; the IP is essentially a method to estimate the difference ϕI − ϕII to hence determine if phase coexistence occurs, and if not, to estimate what change in one of the f fields (say f2) would be needed to attain coexistence. If the two bulk-like phases are separated by an interface present in quantity Mi, then the total free energy, Φ, for this two-phase system is given by

\begin{equation}\Phi = M^I \phi ^I + M^{II} \phi ^{II} + M_i \phi _i,\end{equation}
Φ=MIϕI+MIIϕII+Miϕi,
(7)

where ϕi is the interface contribution to Φ. The IP method assumes the values of ϕi and Mi to be constant. Using this assumption and the fact that M = MI + MII + Mi, Eq. (7) reduces to

\begin{equation}\Phi \approx M^I (\phi ^I - \phi ^{II}) + {\rm const.}\end{equation}
ΦMI(ϕIϕII)+ const .
(8)

The IP method estimates ϕIϕII while maintaining the system as a two-phase system, which is enacted by pinning the interfaces by application of a harmonic field coupled to an order parameter λ that is capable of discriminating between two phases I and II and varies linearly with the quantity MI. If phases I and II in the bulk have order parameter values λI and λII, respectively, then the order parameter of the two phase system can be simply written as

\begin{equation}\lambda = \frac{{M^I }}{M}\lambda ^I + \frac{{M^{II} }}{M}\lambda ^{II} + \frac{{M_i }}{M}\lambda _i = \frac{{M^I }}{M}(\lambda ^I - \lambda ^{II}) + {\rm const.},\end{equation}
λ=MIMλI+MIIMλII+MiMλi=MIM(λIλII)+ const .,
(9)

where the interfacial contribution to λ is assumed to be constant. The IP method then combines Eqs. (8) and (9) to yield

\begin{eqnarray}\Phi \approx \Phi _c \left( {\left\langle \lambda \right\rangle } \right) + {\rm const.} &=& \frac{{M(\phi ^I - \phi ^{II})}}{{\lambda ^I - \lambda ^{II} }}\lambda + {\rm const}.\nonumber\\& =& \frac{{\Phi ^I - \Phi ^{II} }}{{\lambda ^I - \lambda ^{II} }}\lambda + {\rm const.},\end{eqnarray}
ΦΦcλ+ const .=M(ϕIϕII)λIλIIλ+ const .=ΦIΦIIλIλIIλ+ const .,
(10)

which highlights the idea that one essentially assumes that Φc(λ) can be assumed to depend linearly on λ. The potential energy of the system is now modified by adding a pinning potential:

\begin{equation}\Delta U' = \frac{\kappa }{2}(\lambda - a)^2 + {\rm const.},\end{equation}
ΔU=κ2(λa)2+ const .,
(11)

where κ is a spring constant and a is the anchor value. This modifies the probability of microstates of Eq. (3) of the augmented Hamiltonian to

\begin{equation}p'(\lambda) \propto \exp( { - \Phi _c (\lambda) - \beta \Delta U'}),\end{equation}
p(λ)exp(Φc(λ)βΔU),
(12)

which on account of Eqs. (10) and (11) can be shown to be a near-Gaussian distribution in λ around the average value ⟨λ that allows one to estimate ϕIϕII as

\begin{equation}\phi ^I - \phi ^{II} = - \frac{{\beta \kappa (\lambda ^I - \lambda ^{II})}}{M}[ {\langle \lambda \rangle ^\prime - a} ],\end{equation}
ϕIϕII=βκ(λIλII)M[λa],
(13)

where β = 1/kBT (kB is the Boltzmann's constant). All terms on the right hand side of Eq. (13) are knowable from simulations of: (i) each of the single-phases to get λI and λII and (ii) the IP interfacial system to get ⟨λ, noting that κ, a, and M values are preset by the user.

Since a departure from phase coexistence is associated with a non-zero difference, ϕIϕII as obtained by applying Eq. (13), a strategy is now needed to leverage the simulation data available to find the conditions for which ϕI = ϕII. This can be achieved by a model for the free energy ϕ that locally fits the available simulation data and can be extrapolated to nearby conditions. A good approximation to the local behavior of the free-energy function ϕ can be provided by a 2nd order Taylor expansion around a point of interest B:

\begin{eqnarray}\phi &=& \phi _B + ( {\partial \phi /\partial f_1 } )_B \Delta f_1 + ( {\partial \phi /\partial f_2 } )_B \Delta f_2 \nonumber\\&&+\, \tfrac{1}{2}\big( {\partial ^2 \phi /\partial f_1^2 } \big)_B (\Delta f_1)^2 \nonumber\\&&+\, \tfrac{1}{2}\big( {\partial ^2 \phi /\partial f_2^2 } \big)(\Delta f_2)^2 + ( {\partial ^2 \phi /\partial f_1 \partial f_2 } )\Delta f_1 \Delta f_2,\qquad\end{eqnarray}
ϕ=ϕB+(ϕ/f1)BΔf1+(ϕ/f2)BΔf2+122ϕ/f12B(Δf1)2+122ϕ/f22(Δf2)2+(2ϕ/f1f2)Δf1Δf2,
(14)

which by virtue of Eq. (6) reduces to

\begin{eqnarray}\phi &=& \phi _B + \bar x_{1,B} \Delta f_1 + \bar x_{2,B} \Delta f_2 + \tfrac{1}{2}( {\partial \bar x_1 /\partial f_1 } )_B (\Delta f_1)^2 \nonumber\\&&+\, \tfrac{1}{2}( {\partial \bar x_2 /\partial f_2 } )_B (\Delta f_2)^2 + ( {\partial \bar x_1 /\partial f_2 } )_B \Delta f_1 \Delta f_2,\end{eqnarray}
ϕ=ϕB+x¯1,BΔf1+x¯2,BΔf2+12(x¯1/f1)B(Δf1)2+12(x¯2/f2)B(Δf2)2+(x¯1/f2)BΔf1Δf2,
(15)

where |$\bar x_i = \langle {x_i }\rangle$|x¯i=xi is the average of xi, a quantity readily obtained in the simulation at point B (recall that Xi = M xi denotes the extensive properties allowed to fluctuate in the ensemble). If it is assumed that one can also collect second moment data for these fluctuating properties so that we can compute the covariances of variables xi and Xj : cov11, cov22, and cov12.

\begin{equation}{\rm cov}_{ij} = \langle {x_i X_j } \rangle - \langle {x_i } \rangle \langle {X_j } \rangle = \langle {X_i x_j } \rangle - \langle {X_i } \rangle \langle {x_j } \rangle.\end{equation}
cov ij=xiXjxiXj=XixjXixj.
(16)

Note that covij = covji and from |$\bar x_i = \int x_i \exp(- f_1 X_1 \break - f_2 X_2 )dX_1 dX_2 /Q(f_1,f_2)$|x¯i=xiexp(f1X1f2X2)dX1dX2/Q(f1,f2) it is straightforward to show that

\begin{equation}( {\partial \bar x_i /\partial f_j } )_{f_k,k \ne j} = - {\rm cov}_{ij}.\end{equation}
(x¯i/fj)fk,kj= cov ij.
(17)

It follows then that Eq. (15) reduces to32 

\begin{eqnarray}\phi &=& \phi _B + \bar x_{1,B} \Delta f_1 + \bar x_{2,B} \Delta f_2 - \tfrac{1}{2}{\rm cov}_{11,B} (\Delta f_1)^2 \nonumber\\&&-\, \tfrac{1}{2}{\rm cov}_{22,B} (\Delta f_2)^2 - {\rm cov}_{12,B} \Delta f_1 \Delta f_2.\end{eqnarray}
ϕ=ϕB+x¯1,BΔf1+x¯2,BΔf212 cov 11,B(Δf1)212 cov 22,B(Δf2)2 cov 12,BΔf1Δf2.
(18)

Equation (18) is applied to each of two phases I and II so that properties from each phase at point B (⁠|$\bar x_{1,B},\bar x_{2,B}$|x¯1,B,x¯2,B, cov11, B, cov22, B, cov12, B) are known from bulk phase simulations in the ensemble M f1,Bf2,B. Then if f1 and f2 are to correspond to phase coexistence so that Δf2 = f2f2,B and Δf1 = f1f1,B, by setting ϕI = ϕII we can write

\begin{eqnarray}&&\phi _B^I + \bar x_{1,B}^I \Delta f_1 + \bar x_{2,B}^I \Delta f_2 - \tfrac{1}{2}{\rm cov}_{11,B}^I (\Delta f_1)^2\nonumber\\&&\qquad - \tfrac{1}{2}{\rm cov}_{22,B}^I (\Delta f_2)^2 - {\rm cov}_{12,B}^I \Delta f_1 \Delta f_2 \nonumber\\&&\quad = \phi _B^{II} + \bar x_{1,B}^{II} \Delta f_1 + \bar x_{2,B}^{II} \Delta f_2 - \tfrac{1}{2}{\rm cov}_{11,B}^{II} (\Delta f_1)^2 \nonumber\\&&\qquad - \tfrac{1}{2}{\rm cov}_{22,B}^{II} (\Delta f_2)^2 - {\rm cov}_{12,B}^{II} \Delta f_1 \Delta f_2.\end{eqnarray}
ϕBI+x¯1,BIΔf1+x¯2,BIΔf212 cov 11,BI(Δf1)212 cov 22,BI(Δf2)2 cov 12,BIΔf1Δf2=ϕBII+x¯1,BIIΔf1+x¯2,BIIΔf212 cov 11,BII(Δf1)212 cov 22,BII(Δf2)2 cov 12,BIIΔf1Δf2.
(19)

Note that if |$\phi _B^I - \phi _B^{II}$|ϕBIϕBIIis also known, then Eq. (19) would allow us to solve for Δf2 if Δf1 is specified; i.e., to estimate the value of f2 corresponding to phase coexistence for any new value of f1 (away from f1,B). For the interfacial simulation, since the state is user-specified by the value of one of the thermodynamic fields, say property f1, then f1 = f1,B always, and the sought-after coexistence state will be determined by a (yet unknown) value of f2. Equation (19) then simplifies and can be rearranged as

\begin{equation}\Delta f_2 = \frac{{\phi _B^I - \phi _B^{II} }}{{ - \big(\bar x_{2,B}^I - \bar x_{2,B}^{II}\big) + \tfrac{1}{2}\big({\rm cov}_{22,B}^I - {\rm cov}_{22,B}^{II}\big)\Delta f_2 }},\end{equation}
Δf2=ϕBIϕBIIx¯2,BIx¯2,BII+12 cov 22,BI cov 22,BIIΔf2,
(20)

which is second order in Δf2 and hence solvable by the quadratic formula. The difference |$\phi _B^I - \phi _B^{II}$|ϕBIϕBII in Eq. (20) is obtained from the IP simulation via Eq. (13), while the |$\bar x_{2,B}$|x¯2,B and cov22, B properties come from the separate single-phase runs. If one were to neglect the contribution of second order terms (covariances), Eq. (20) reduces then to

\begin{equation}\Delta f_2 = \frac{{\phi _B^{II} - \phi _B^I }}{{ - \big(\bar x_{2,B}^{II} - \bar x_{2,B}^I\big)}}.\end{equation}
Δf2=ϕBIIϕBIx¯2,BIIx¯2,BI.
(21)

Equation (21) is a generalized version of the Newton-Raphson method for finding roots that Pedersen et al. suggested as iteration scheme to estimate the f2 value that satisfies phase equilibrium. Equation (20) has the potential advantage of reducing the number of iterations that may be needed for convergence. The 2nd order extrapolation framework also allows predicting the coexisting properties |$\bar x_1$|x¯1 and |$\bar x_2$|x¯2 for each phases if the correction Δf2 [as found in Eq. (20) or (21)] is to be applied to the properties measured at conditions “B;” namely [adapting Eqs. (14) and (15) from Ref. 32],

\begin{equation}\bar x_k^J = \bar x_{k,B}^J - {\rm cov}^{J}_{k2,B} \Delta f_2,\end{equation}
x¯kJ=x¯k,BJ cov k2,BJΔf2,
(22)

where k = 1 or 2 and J = I or II.

Of course, one can use more advanced methods for free-energy extrapolation which can use either polynomial functions32–35 or reweighted histograms;36–42 however, the large system sizes needed for the simulations would render histogram reweighting inefficient. Hence, the relatively small corrections in Δf2 typically associated with the needed single-phase simulations that remained (meta) stable, makes Eq. (20) or (21) more suitable and practical.

In the isothermal-isobaric NPT ensemble for a single component, one would have that M = N (number of particles), if we choose f1 = β, then f2 = βP (P = pressure), and ϕ = βμ (where μ is chemical potential), so that |$\bar x_1 = u$|x¯1=u(molar energy), |$\bar x_2 = v$|x¯2=v (molar volume), and cov11, cov22, and cov12 would relate to the constant-pressure heat capacity, isothermal compressibility, and thermal expansion coefficient, respectively. Equations (13) and (20) then become

\begin{equation}\beta \big( {\mu _B^I - \mu _B^{II} } \big) = - \frac{{\beta \kappa (\lambda ^I - \lambda ^{II})}}{N}[ {\langle \lambda \rangle ^\prime - a} ],\end{equation}
βμBIμBII=βκ(λIλII)N[λa],
(23)
\begin{equation}\beta \Delta P = \frac{{\beta \big(\mu _B^I - \mu _B^{II}\big)}}{{ - \big(v_B^I - v_B^{II}\big) + \tfrac{1}{2}\big({\rm cov}_{vV,B}^I - {\rm cov}_{vV,B}^{II}\big)\beta \Delta P}}.\end{equation}
βΔP=βμBIμBIIvBIvBII+12 cov vV,BI cov vV,BIIβΔP.
(24)

For the semigrand isothermal-isobaric ensemble or SGNPT ensemble for a binary mixture, the total number of particles N, the temperature, pressure, and the chemical potential difference between component 1 and component 2 of the mixture, ξ = β (μ2μ1) are fixed. While in this case there are 3 specifiable f fields, we let both the temperature and pressure to be set (fixed) so that one would effectively considers f1 to be a vector f1 = (β, βP), a distinction that has no effect on Eqs. (13) and (20). We have then M = N, ϕ = βμ2, |$\bar x_1 = v$|x¯1=v, f2 = ξ, |$\bar x_2 = y_1$|x¯2=y1, where yi is the mole fraction of component i. Equations (13) and (20) then become

\begin{equation}\beta \big( {\mu _{2,B}^I - \mu _{2,B}^{II} } \big) = - \frac{{\beta \kappa (\lambda ^I - \lambda ^{II})}}{N}[ {\langle \lambda \rangle ^\prime - a} ],\end{equation}
βμ2,BIμ2,BII=βκ(λIλII)N[λa],
(25)
\begin{equation}\Delta \xi = \frac{{\beta \big(\mu _{2,B}^I - \mu _{2,B}^{II}\big)}}{{ - \big(y_{1,B}^I - y_{1,B}^{II}\big) + \tfrac{1}{2}\big({\rm cov}_{y_1 N_1,B}^I - {\rm cov}_{y_1 N_1,B}^{II}\big)\Delta \xi }}.\end{equation}
Δξ=βμ2,BIμ2,BIIy1,BIy1,BII+12 cov y1N1,BI cov y1N1,BIIΔξ.
(26)

It should be noted that we could have chosen the sought-after (f2) and specified (f1) properties differently, e.g., to seek out the coexistence temperature for a given pressure in the NPT ensemble or the coexistence pressure of a given chemical potential difference in the SGNPT ensemble. However, for hard-core systems the specification of pressure is the more suitable or convenient choice.

To distinguish between isotropic and solid phases, several translational order parameters can be used such as the bond-orientational order parameters Q643 and a collective order parameter |ρk| derived from the structure factor, S(k), which were used in the original implementation of the IP method by Pedersen et al.29 The bond order parameter, Q6 is defined as

\begin{equation}Q_6 = \frac{{4\pi }}{{13}}\left[ {\sum\limits_{m = - 6}^6 {\left| {\bar Q_{6m} (\vec r)} \right|^2 } } \right]^{\frac{1}{2}},\end{equation}
Q6=4π13m=66Q¯6m(r)212,
(27)

where |$\bar Q_{6m} (\vec r)$|Q¯6m(r) is given by

\begin{equation}\bar Q_{6m} (\vec r) = \frac{1}{{N_b }}\sum\limits_{bonds} {Y_{6m} (\vec r)},\end{equation}
Q¯6m(r)=1NbbondsY6m(r),
(28)

where |$Y_{6m} (\vec r)$|Y6m(r) are spherical harmonics for the position vector |$\vec r$|r and the summation is over all the bonds (total number Nb) present in the system. Two particles are considered to be bonded if they are separated by a distance less than 1.2 times the position of the first peak of radial distribution function. The value of Q6 approaches zero for structures lacking translational order while it attains a high value for different types of lattice structures;43 it is hence a good metric to detect translational order for different types of lattice structures. |ρk| has similar characteristics to Q6 in terms of generality and it is defined as

\begin{equation}\left| {\rho _k } \right| = \left| {N^{ - \frac{1}{2}} \sum\limits_{j = 1}^N {\exp( - i{\bf k} \cdot \vec r_j } } \right|,\end{equation}
ρk=N12j=1Nexp(ik·rj,
(29)

where k = (2πnx/LX, 2πny/LY, 0). The values LX and LY are the box lengths along the X and Y axes, respectively. The integers nx and ny are chosen so that the wave vector, k corresponds to a Bragg peak and they can be readily found by trial and error using a sample solid phase. The value of |ρk| is of the order of 1 and N1/2 for the liquid and solid phase, respectively.

Here, we also propose the specific volume as alternative order parameter,

\begin{equation}\upsilon = \frac{V}{{NV_r }},\end{equation}
υ=VNVr,
(30)

where V is the volume of the box, N is the number of particles, and Vr is a reference volume which is taken to be the volume of a single particle in the system (for the mixture case it is the volume of a single cube). Note that for our systems the order parameter υ satisfies the conditions that: (1) It distinguishes the isotropic and solid phase and (2) it varies approximately linearly with the amount of solid particles in the system. In terms of computational cost, υ has the appeal that its value is known at all times for “free;” this is in contrast with Q6 values, for instance, whose update requires a special calculation after any particle translates in the system, adding to the expense of the IP method. The performance of these 3 order parameters when applied with IP method is explained in detail in Sec. III.

In selecting the pinning potential parameters, the value of a should be such that the pure phase slabs have bulk properties and the value of κ should be large enough to keep the system in a two-phase configuration. The formulas suggested in Ref. 30 are

\begin{equation}a \approx \tfrac{1}{2}( {\lambda ^I + \lambda ^{II} } ),\end{equation}
a12(λI+λII),
(31)
\begin{equation}\kappa \approx \frac{{N_z^2 }}{{\beta ( {\lambda ^I - \lambda ^{II} } )^2 }},\end{equation}
κNz2β(λIλII)2,
(32)

where Nz is the number of ordered planes in z-direction. We used Eq. (31) estimate a and Eq. (32) to provide the order of magnitude of κ for λ = Q6, |ρk|, and υ in the pure polyhedra systems. An alternative criterion for the selection of κ can also be devised as follows: The variance in the p(λ) distribution for the IP run (=kBT/κ) should be significantly smaller than that in the unconstrained p (λ) for the stable phase (say, |$s_\lambda ^2$|sλ2) or else λ fluctuations would not be sufficiently suppressed to prevent the disappearance of the metastable phase; hence

\begin{equation}\kappa > \big( {\beta s_\lambda ^2 } \big)^{ - 1}.\end{equation}
κ>βsλ21.
(33)

We used Eq. (33) when using υ as order parameter for the mixture system choosing κ to be an order of magnitude larger than |$1/\beta s_\upsilon ^2$|1/βsυ2.

For completeness and clarity we briefly summarize the main steps in the IP method as applied to solid-fluid equilibrium (further details can be found in Refs. 29 and 30):

  • Build single-phase configurations for phases I and II in boxes elongated along the z-axis.

  • Perform simulations in the Mf1f2 ensemble on each single-phase system using the same values of the properties M, f1, and f2, where f2 is an estimate of the sought-after coexistence value. These conditions define the state “B” in Eq. (20) and simulations should yield the converged values for the order parameter λ and the properties |$\bar x_1,\bar x_2,{\rm cov}_{11},\,{\rm cov}_{22},\, {\rm cov}_{12}$|x¯1,x¯2, cov 11, cov 22, cov 12 for each of the phases. For the solid phase, the unstrained cross section of the box (e.g., along the X and Y axes) must have been identified.

  • Create a two-phase interfacial system with box lengths in X and Y axes that match those of the unstrained solid.

  • Simulate this interfacial system in the Mf1f2 ensemble with the pinning potential of Eq. (11) to evaluate ⟨λ, and allow estimation of ϕIϕII via Eq. (13).

  • Estimate the coexistence value of f2 via extrapolations from the current state conditions B via Eq. (20) or (21).

Iterate the process above until a desired level of convergence in the coexistence f2 value is attained.

In this section, we apply the IP method to three polyhedral shapes; COs, TOs, and RDs. They are representative of a class of polyhedral particles that form crystals with well-defined translational and orientational order at close packing but at the order-disorder transition pressure, the isotropic phase coexists with a rotator or plastic mesophase wherein particles are translationally ordered but essentially orientationally disordered.4,6,7 The IP method is applied in the NPzT ensemble to estimate the coexistence pressure, Pcoex, between isotropic and rotator phases. We report our results in terms of reduced coexistence pressure given by |$P_{coex}^*$|Pcoex* = βPcoexσ3/8, where σ is the diameter of the circumscribing sphere (specific to each particle geometry). In our simulations, the number of particles N is chosen as 4356, 4536, and 4312 for COs, TOs, and RDs, respectively, and correspond to 11 × 11 × 36, 9 × 9 × 28, and 7 × 7 × 22, lattice unit cells in their solid phases. At any given pressure P, we first perform the bulk phase simulations to obtain λiso and λrot for the isotropic and rotator phases, respectively. We then prepare the two-phase configuration using the procedure mentioned in Espinosa et al.44 The box lengths X and Y in the corresponding x and y directions of the two-phase system are the same as those of the unrestrained rotator phase at a pressure P. NPzT Monte Carlo (MC) simulations are performed to obtain the average value of ⟨λ which is used to obtain Δμ(P) through Eq. (23). Each MC simulation consists of 2 × 106 MC cycles and each MC cycle consists of N translational, N rotational, N/10 flip, and 1 volume move attempts. Flip moves attempt to rotate a chosen particle to a random orientation in the plane perpendicular to its present orientation. All trial moves are accepted according to the Metropolis criterion45 which requires ruling out overlaps between any two particles (via the separating axes theorem46). For all three shapes, we have applied the method using three different order parameters Q6,43 |ρk|, and υ. All three order parameters can distinguish the rotator phase from the isotropic phase and vary linearly with the amount of translationally ordered particles. For the wave vector k needed for |ρk|, the integers nx and ny are chosen to maximize the contrast between liquid and rotator phase; in particular we used (nx, ny) = (22, 0), (9, 9), and (14, 0) for the COs, TOs, and RDs, respectively. The harmonic field parameters a and κ used for Q6, |ρk|, and υ for all three shapes are listed in Table I.

Table I.

Pinning potential parameters, a and κ used for different order parameters for COs, TOs, and RDs.

ShapeOrder parameteraκ
CO Q6 0.2 80 000 
  υ 2.06 400 000 
  k15.0 10 
TO Q6 0.2 80 000 
  υ 2.05 400 000 
  |ρk25.5 10 
RD Q6 0.2 80 000 
  υ 2.07 400 000 
  |ρk25.5 10 
ShapeOrder parameteraκ
CO Q6 0.2 80 000 
  υ 2.06 400 000 
  k15.0 10 
TO Q6 0.2 80 000 
  υ 2.05 400 000 
  |ρk25.5 10 
RD Q6 0.2 80 000 
  υ 2.07 400 000 
  |ρk25.5 10 

To test the performance of each order parameter, we apply the IP method for a TO system at a pressure value of 2.5. This pressure is chosen as it is the approximate coexistence value obtained by Umang et al.4 using an interfacial NVT ensemble. We obtain equilibrated configurations for each order parameter by performing NPzT MC simulations of a two-phase system. Figure 1(a) shows the profile of translational order over the distance z for an equilibrated configuration obtained at applying the IP method for each order parameter. We also show the near Gaussian probability distributions obtained for each order parameter in Fig. 1(b). To plot all three distributions in the same figure, we use reduced order parameter, λred = (λa)/|Δλ|, where |Δλ| is the absolute value of the difference between the bulk phase λ values of the isotropic and rotator phases. The profile in Fig. 1(a) shows that Q6 and υ are successful in keeping the system in a two-phase with an interface and the distributions in Fig. 1(b) shows that at P* = 2.5 the rotator phase is the stable phase. The latter follows from the positions of the distributions that show the system pulling towards Q6 values above the anchor point and υ values below the anchor point. The near Gaussian distribution centered close to a obtained with |ρk| is spurious as the system fails to maintain a two-phase state and completely transforms into a rotator phase. It appears that the rotator phase lowers its value of |ρk| by a collective rearrangement of particles while maintaining a high translational order. While using |ρk|, we have observed a similar issue for the other two shapes whenever the pressure employed was an overestimation of the coexistence pressure (i.e., when the liquid phase is metastable). Hence, |ρk| is not a good choice for the systems under study.

FIG. 1.

Application of the IP method to the TO system at P* = 2.5 using separately Q6, |ρk|, and υ as order parameters. (a) Profile of Q6 (z) versus z/σ for three equilibrated configurations and (b) the probability distributions versus λred.

FIG. 1.

Application of the IP method to the TO system at P* = 2.5 using separately Q6, |ρk|, and υ as order parameters. (a) Profile of Q6 (z) versus z/σ for three equilibrated configurations and (b) the probability distributions versus λred.

Close modal

By using both Q6 and υ, Fig. 2 shows the results obtained by using Eqs. (23) and (24) to iteratively obtain the converged values of |$P_{coex}^*$|Pcoex* starting from different initial pressures|$P_0^*$|P0*. To validate the IP method, two values of initial pressures are used for each shape, one below and one above |$P_{coex}^*$|Pcoex*. We stop our simulations when the value of Δμ(P) is close to zero (Δμ(P) = ±0.005). In all cases, we obtain the coexistence pressure in fewer than 4 iterations as shown in Fig. 2. Table II compares the values of |$P_{coex}^*$|Pcoex*obtained using different |$P_0^*$|P0* for both Q6 and υ and shows that the converged values of |$P_{coex}^*$|Pcoex* obtained by using the proposed order parameter υ matches well with that based in Q6. However, a MC cycle in the IP run is approximately 3–4 times longer when using Q6 than when using υ as order parameter, due to the calculation of spherical harmonics after each MC translation move. Such a disparity is dependent on the type of machine used (in our case a Xeon® 5500, 3.06 GHz processor node) and implementation details (in our updates of Q6 the cost was minimized by constraining any new calculation to only those particles whose neighbor's positions had been affected by a translation move).

FIG. 2.

P* values obtained at different iterations of the IP method for COs, TOs, and RDs. Results are shown for two different initial pressure values used for each of the two order parameters Q6 and υ.

FIG. 2.

P* values obtained at different iterations of the IP method for COs, TOs, and RDs. Results are shown for two different initial pressure values used for each of the two order parameters Q6 and υ.

Close modal
Table II.

The values of |$P_0^*$|P0* and |$P_{coex}^*$|Pcoex* using Q6 and υ as order parameters in the IP method for COs, TOs, and RDs. Literature |$P_{coex}^*$|Pcoex* values obtained using different methods are also shown.

    Direct 
 Order  interfacialThermodynamic
Shapeparameter|$P_0^*$|P0*|$P_{coex}^*$|Pcoex*method4 integration6 
CO Q6 2.7 2.94 N/A 2.9 
    3.2 2.94     
  υ 2.7 2.94     
    3.2 2.95     
TO Q6 2.0 2.18 2.5 2.18 
    2.5 2.19     
  υ 2.0 2.19     
    2.5 2.18     
RD Q6 3.2 3.37 3.6 N/A 
    3.6 3.37     
  υ 3.2 3.38     
    3.6 3.38     
    Direct 
 Order  interfacialThermodynamic
Shapeparameter|$P_0^*$|P0*|$P_{coex}^*$|Pcoex*method4 integration6 
CO Q6 2.7 2.94 N/A 2.9 
    3.2 2.94     
  υ 2.7 2.94     
    3.2 2.95     
TO Q6 2.0 2.18 2.5 2.18 
    2.5 2.19     
  υ 2.0 2.19     
    2.5 2.18     
RD Q6 3.2 3.37 3.6 N/A 
    3.6 3.37     
  υ 3.2 3.38     
    3.6 3.38     

We compared our results with available literature values for |$P_{coex}^*$|Pcoex* obtained from a thermodynamic integration scheme6 and a direct coexistence method applied in the NVT ensemble.4 We find that our results match well with those considered most accurate, i.e., from thermodynamic integration, thus validating our implementation. The interfacial NVT simulation overestimate the |$P_{coex}^*$|Pcoex* values we now find, even though the former can be considered a limiting case of the IP method that uses υ as order parameter with an infinitely strong pinning potential. The reasons for this discrepancy are twofold. First, interfacial simulations, like most other simulation methods, are sensitive to finite size effects and in particular to the number of lattice unit cells along the crystal cross section; it turns out that the cross sections used in Ref. 4 were smaller than those adopted here. Second, the NVT simulation was not iterated to better match the cross section of the unstrained rotator phase for a given pressure. Hence, the system was not able to fully relax in the x y plane and the rotator phase might have been under some stress.

For the systems studied in this section, we find that the |$P_{coex}^*$|Pcoex* value obtained using our 2nd order extrapolation technique [Eq. (20) or (24)] are similar to the ones obtained using the 1st order scheme [Eq. (21)] embodied by the Newton Raphson's method used by Pedersen et al.29 This occurs because the initial pressures chosen for three shapes were close enough to the coexistence value that the variation of chemical potential difference with pressure was approximately linear.

In implementing the SGNPT ensemble for the mixture of cubes and spheres, the total number of particles is fixed at N = 3920 so that the solid phase would contain 7 × 7 × 20 fcc unit cells. The temperature is fixed at a reduced value of T* = (βɛ)−1 = 1 (where ɛ is an arbitrary parameter). The composition of the system is not fixed but allowed to fluctuate in response to the specified value of the chemical potential difference ξβΔμ = β(μ2μ1), letting species 1 and 2 correspond to cubes and spheres, respectively. The dimensionless osmotic pressure is given in reduced units as P* = βPσ3, where σ is the cube edge length. Runs for any given system and state involved 106 MC cycles for equilibration and 2 × 106 MC cycles for production. Each MC cycle consisted on average of N translational, N rotational (applied to cubes only), N/20 swap, N/10 mutation, and 2 volume move attempts. Swap moves involved picking randomly two particles of different species and attempting to swap their positions (while keeping the orientation for the cubes). Mutation moves involved picking a particle at random and attempting to mutate its original identity into another species type, keeping its old position (and for a sphere → cube mutation, choosing a random orientation for the cube). In the interfacial runs only the box length along the z axis is allowed to fluctuate via the volume moves (i.e., effectively a SGNPzT ensemble). All trial moves are accepted according to the Metropolis criterion45 which required ruling out overlaps between any two spheres, any cube-sphere pair (via the Arvo's algorithm47) and any two cubes (via the separating axes theorem46).

Both hard cubes and hard spheres exhibit a first-order phase transition from an isotropic phase at low pressures to a solid at high pressures. The pressure-composition phase diagram for sphere + cube mixtures has been reported for sphere diameter (d) to cube edge length ratios of d/σ = 119 and for d/σ = 1.2332 and in both cases an eutectic behavior was observed. For this study we only selected one point in the phase diagram for d/σ = 1.23 where an fcc sphere-rich solid phase coexists with an isotropic phase at P* = 8.0 (lying in between the order-disorder pressure for the pure spheres and the mixture eutectic pressure).

For these simulations we first implemented |ρk| as order parameter but found that, like in the case of the pure polyhedra reported before, it was not able to stabilize two-phase interfacial structures at conditions when the liquid phase was metastable, thus rendering the IP method unworkable. Further, since in Sec. III A we already demonstrated that υ is an inexpensive order parameter that is similarly workable as Q6 for hard-core systems, we only implemented calculations with λ = υ.

Table III shows the results from the IP method for different input values of ξ (whose value at coexistence is the sought-after unknown). The initial guesses of ξ = −2.0 and −2.8 correspond to states where either the liquid phase or solid phase is strongly metastable, respectively, while the ξ = −2.33 and −2.36 can be seen as near-coexistence, improved estimates from successive iterations of the IP extrapolation. It appears that convergence to the coexistence value of ξ can only be achieved with at most 3 significant figures as the extrapolations from ξ = −2.33 and −2.36 tend to overshoot in opposite directions, reflecting the difficulty to precisely estimate differences (⟨υ⟩ − a). We hence estimate the coexistence value as ξ = −2.35 ± 0.01, which is the average between the last trial value (−2.36) and the last extrapolated value (−2.34).

Table III.

Results for cube-sphere mixture at P* = 8.0. Spring constant used was κ = 5.9 × 105and the constant a was chosen to be the average of υ for the liquid and solid phases found in the single-phase simulations.

Input ξ−2.00−2.80−2.33−2.36
|$y_1^{solid}$|y1solid, cube fraction in solid 0.072 0.240 0.122 0.124 
|$y_1^{liquid}$|y1liquid, cube fraction in liquid 0.255 0.388 0.308 0.313 
|$\beta ( {\mu _2^{solid} - \mu _2^{liquid} } )$|β(μ2solidμ2liquid), Eq. (25) −0.061 0.069 −0.0056 0.0038 
Constant a 1.8250 1.8602 1.8427 1.8417 
υ⟩ 1.8229 1.8637 1.8425 1.8419 
New ξ via Eq. (21) −2.333 −2.330 −2.360 −2.340 
New ξ via Eq. (26) −2.337 −2.365 −2.360 −2.340 
Input ξ−2.00−2.80−2.33−2.36
|$y_1^{solid}$|y1solid, cube fraction in solid 0.072 0.240 0.122 0.124 
|$y_1^{liquid}$|y1liquid, cube fraction in liquid 0.255 0.388 0.308 0.313 
|$\beta ( {\mu _2^{solid} - \mu _2^{liquid} } )$|β(μ2solidμ2liquid), Eq. (25) −0.061 0.069 −0.0056 0.0038 
Constant a 1.8250 1.8602 1.8427 1.8417 
υ⟩ 1.8229 1.8637 1.8425 1.8419 
New ξ via Eq. (21) −2.333 −2.330 −2.360 −2.340 
New ξ via Eq. (26) −2.337 −2.365 −2.360 −2.340 

Note that the 2nd order extrapolation scheme for ξ [i.e., via Eq. (26)] only provides a slight improvement over the 1st order scheme (that neglects covariances) when the initial guess for ξ is far enough from the correct value that nonlinear effects are non-negligible, as shown by the extrapolated values obtained for ξ = −2.8. For the coexistence ξ = −2.35 value we can use Eq. (22) to estimate the cube mole fractions of |$y_1^{solid}$|y1solid = 0.123 and |$y_1^{liquid}$|y1liquid = 0.311, which agree well with the results found from the Gibbs-Duhem integration scheme reported in Ref. 32 (that started at the solid-liquid transition point for pure spheres), namely: ξ = −2.354, |$y_1^{solid}$|y1solid = 0.123, and |$y_1^{liquid}$|y1liquid = 0.310. In fact, Eq. (22) would give good estimates even if the base point were quite far away; e.g., taking the ξ = −2.8 data as B point would give |$y_1^{solid}$|y1solid = 0.128 and |$y_1^{liquid}$|y1liquid = 0.306 which constitute much improved estimates compared to the original results for those conditions (3rd column in Table III).

Figure 3 shows the histograms of the specific volume υ in the IP simulations for selected values of the input ξ value. The curves show that the distributions of υ are very narrow, almost delta-function like, indicating that the pinning spring potential value selected as per formula (Eq. (33)) (κ = 5.9 × 105 in our units) was very high, only allowing fluctuations in υ of ±1% relative to (υliquid − υsolid). We experimented with different values of κ and found (results not shown) that the convergence rate (and thus the efficiency) of the IP method as implemented here (with λ = υ) was rather sensitive to the choice of κ and we conjecture that using a translational-order parameter would have provided a more robust convergence. Note also that in choosing κ it is safer to err in the side of excess than in the alternative, especially at far-from-equilibrium conditions, since a too-soft pinning potential may in such cases allow large volume fluctuations that move the interfacial system toward its stable one-phase state. The pinning potential could be made softer as the coexistence state is approached and the driving force toward a one-phase state lessens. Figure 4 shows a sample configuration for the two-phase state of the system during the IP simulation with ξ = −2.36 which clearly shows a significant concentration of cubes (∼12.4%) which are positioned into the fcc lattice sites of the sphere-rich solid phase while exhibiting diverse orientations.

FIG. 3.

Probability density distributions of υ obtained in the IP simulations of the cube-sphere mixture in the semigrand NPzT ensemble at P* = 8 for 3 different choices of the inter-species chemical potential difference ξ (values shown in the legend). The vertical lines correspond to the pinning value of parameter a. Note that for ξ = −2.0 the system pushes left of a to lower υ's (the stable solid phase) while for ξ = −2.8 the system pushes right of a to higher υ's (the stable liquid phase).

FIG. 3.

Probability density distributions of υ obtained in the IP simulations of the cube-sphere mixture in the semigrand NPzT ensemble at P* = 8 for 3 different choices of the inter-species chemical potential difference ξ (values shown in the legend). The vertical lines correspond to the pinning value of parameter a. Note that for ξ = −2.0 the system pushes left of a to lower υ's (the stable solid phase) while for ξ = −2.8 the system pushes right of a to higher υ's (the stable liquid phase).

Close modal
FIG. 4.

Snapshot of the semigrand NPzT interfacial simulation for the cube-sphere mixture for P* = 8 and ξ = −2.33 showing the solid phase region in the leftmost half and the liquid region in the rightmost half of the box.

FIG. 4.

Snapshot of the semigrand NPzT interfacial simulation for the cube-sphere mixture for P* = 8 and ξ = −2.33 showing the solid phase region in the leftmost half and the liquid region in the rightmost half of the box.

Close modal

Of course, the implementation of IP with λ = υ relies on the presence of a significant difference in υ between the two phases, a scenario that may not always occur between the solid and liquid phases in a mixture. In such cases, other phase distinguishing properties (besides translational order) could prove useful such as composition. In fact, we also conduced additional tests with alternative order parameters that could better resolve the two phases (results not shown); one such example was to add a second harmonic term to the pinning potential so that: ΔU = κ/2(υ − a)2 + κ/2(y1a)2 where a is chosen to be the average composition (y1) of those for the two pure phases. For the cube-sphere mixture, however, the simpler υ-based pinning potential (i.e., κ = 0) proved to be similarly effective.

We implemented a generalized version of the IP method and applied it to locate solid-liquid coexistence conditions for hard-core single-component polyhedral particles and a mixture of cubes and spheres. For pure systems of COs, TOs, and RDs, we applied the IP method in the NPzT ensemble to estimate the coexistence pressure, |$P_{coex}^*$|Pcoex*, between the isotropic and rotator phases. For this purpose we tested 3 different order parameters |ρk|, Q6, and υ (specific volume). We find that when using |ρk| the system is unable to stay in a two-phase coexistence state whereas using Q6 (a more robust translational order parameter) does allow for accurate estimates of |$P_{coex}^*$|Pcoex*. The use of υ, which unlike Q6 involves no computational cost to update, also allows for reliable estimates of |$P_{coex}^*$|Pcoex* that match closely those obtained using Q6 and those reported in the literature obtained using thermodynamic integration.

For the mixture of cubes and spheres, we applied the IP method in the semigrand NPT ensemble to estimate isotropic-fcc crystal phase coexistence conditions at P* = 8.0. Similar to single-component systems, the use of |ρk| is found to be ineffectual to consistently restrict the system to an interfacial coexistence state. The use of υ as order parameter was found to be workable and the coexistence properties converged to values that agree with results formerly obtained using a Gibbs-Duhem integration method. In this case, the convergence of the IP iterations was found to be somewhat sensitive to the choice of the pinning spring constant.

This work expands the IP methodology in the following aspects: (1) A generalized framework is provided that allows its implementation for different types of ensembles and for systems with more than one component (the original implementation was restricted to a single ensemble and pure components), (2) a second order extrapolation scheme is presented to estimate coexistence conditions which may converge faster to the solution than the first order scheme originally presented, (3) an alternative criterion is given for choosing the pinning spring constant [Eq. (33)], and (4) the workability of inexpensive order parameters like specific volume or composition is demonstrated, which can be used in the cases where their values change non-negligibly between the two coexisting phases. In those cases, these simple global properties may provide a convenient alternative to translational order parameters when the latter are not easy to implement, renders simulations very expensive, or if translation order is not captured by conventional parameters. On the application side, this work illustrates the use of the IP method to entropy-driven order-disorder transitions for non-trivial colloidal model systems, a departure from the temperature-driven phase transitions of atomic systems studied in Refs. 29 and 30.

The advantages and disadvantages of interfacial simulation approaches (like the IP method) in comparison with other methods to simulate phase coexistence have already been reviewed in previous papers (see, for example, Refs. 8, 16, and 30). Indirect approaches like thermodynamic integration often require more computational effort to simulate multiple state points along a chosen integration path (using reference points and coupling parameters that may be non-trivial to implement as for the pure polyhedral systems studied here6 and mixtures therefrom). Methods that estimate free energies or enact free-energy equality between phases via insertion-deletion of molecules would be inapplicable to solid phases. Compared to other interfacial methods, the advantage of the IP approach centers on its ability to obtain the free-energy difference between the phases present in the system, information that can be used to iteratively target coexistence conditions. But as with all interfacial methods where finite bulk-like phase regions coexist in the same box, the IP method can exhibit larger finite size effects (e.g., associated with tail-corrections for long-range interactions) than in methods where the bulk phases are in separate boxes. Also, simulating solid-solid equilibria is troublesome via any interfacial method, as mechanical straining would necessarily occur for phases with incompatible cross-sections.

While the IP method could be advantageously used in conjunction with a number of other methods that target phase coexistence, it appears to be most complementary to the Gibbs-Duhem integration method (GDI) when dealing with solid-liquid transitions.35,48 A single GDI step only involves two single-phase simulations and it is hence significantly cheaper than an IP calculation (that requires a few iterations of both single-phase and interfacial simulations and for larger system sizes). However, GDI does not have a built-in method to jump start it and after tracing many points along a coexistence line, errors may accrue that can be difficult to diagnose or correct. In this context, the IP method may provide GDI with both a starting point and a checkpoint to correct any systematic deviation. Note also that both IP (see Sec. II C) and GDI32 can be seen to rely on the same extrapolation scheme of free energies.

The authors are grateful to Dr. U. Agarwal and M. Khadilkar for useful exchanges. This work was supported by the U.S. National Science Foundation Grant No. CBET-1402117.

1.
G. M.
Whitesides
and
M.
Boncheva
,
Proc. Natl. Acad. Sci. U.S.A.
99
,
4769
(
2002
).
2.
S. C.
Glotzer
and
M. J.
Solomon
,
Nat. Mater.
6
,
557
(
2007
).
3.
C.
Burda
,
X.
Chen
,
R.
Narayanan
, and
M. A.
El-Sayed
,
Chem. Rev.
105
,
1025
(
2005
).
4.
U.
Agarwal
and
F. A.
Escobedo
,
Nat. Mater.
10
,
230
(
2011
).
5.
P. F.
Damasceno
,
M.
Engel
, and
S. C.
Glotzer
,
Science
337
,
453
(
2012
).
6.
A. P.
Gantapara
,
J.
de Graaf
,
R.
van Roij
, and
M.
Dijkstra
,
Phys. Rev. Lett.
111
,
015501
(
2013
).
7.
U.
Agarwal
and
F. A.
Escobedo
,
J. Chem. Phys.
137
,
024905
(
2012
).
8.
C.
Vega
,
E.
Sanz
,
J. L.
Abascal
, and
E. G.
Noya
,
J. Phys.: Condens. Matter
20
,
153101
(
2008
).
9.
B.
Widom
,
J. Chem. Phys.
39
,
2808
(
1963
).
10.
W. G.
Hoover
,
J. Chem. Phys.
49
,
3609
(
1968
).
11.
W. G.
Hoover
,
J. Chem. Phys.
55
,
1128
(
1971
).
12.
D.
Alfè
and
M.
Gillan
,
Phys. Rev. B
68
,
205212
(
2003
).
13.
A. J. C.
Ladd
and
L. V.
Woodcock
,
Chem. Phys. Lett.
51
,
155
(
1977
).
14.
R. L.
Davidchack
and
B. B.
Laird
,
J. Chem. Phys.
108
,
9452
(
1998
).
15.
R.
Sibug-Aga
and
B. B.
Laird
,
J. Chem. Phys.
116
,
3410
(
2002
).
16.
E. G.
Noya
,
C.
Vega
, and
E.
de Miguel
,
J. Chem. Phys.
128
,
154507
(
2008
).
17.
T.
Zykova-Timan
,
J.
Horbach
, and
K.
Binder
,
J. Chem. Phys.
133
,
014705
(
2010
).
18.
J. R.
Morris
and
X.
Song
,
J. Chem. Phys.
116
,
9352
(
2002
).
19.
M. R.
Khadilkar
,
U.
Agarwal
, and
F. A.
Escobedo
,
Soft Matter
9
,
11557
(
2013
).
20.
M. R.
Khadilkar
and
F. A.
Escobedo
,
J. Chem. Phys.
137
,
194907
(
2012
).
21.
O. A.
Karim
,
P. A.
Kay
, and
A. D. J.
Haymet
,
J. Chem. Phys.
92
,
4634
(
1990
).
22.
J. B.
Sturgeon
and
B. B.
Laird
,
Phys. Rev. B
62
,
14720
(
2000
).
23.
T.
Bryk
and
A. D. J.
Haymet
,
J. Chem. Phys.
117
,
10258
(
2002
).
24.
H.
Nada
,
J.
van der Eerden
, and
Y.
Furukawa
,
J. Cryst. Growth
266
,
297
(
2004
).
25.
R. G.
Fernandez
,
J. L. F.
Abascal
, and
C.
Vega
,
J. Chem. Phys.
124
,
144506
(
2006
).
26.
R.
Handel
,
R.
Davidchack
,
J.
Anwar
, and
A.
Brukhno
,
Phys. Rev. Lett.
100
,
036104
(
2008
).
27.
M. M.
Conde
,
C.
Vega
, and
A.
Patrykiejew
,
J. Chem. Phys.
129
,
014702
(
2008
).
28.
M. S. G.
Razul
and
P. G.
Kusalik
,
J. Chem. Phys.
134
,
014710
(
2011
).
29.
U. R.
Pedersen
,
F.
Hummel
,
G.
Kresse
,
G.
Kahl
, and
C.
Dellago
,
Phys. Rev. B
88
,
094101
(
2013
).
30.
U. R.
Pedersen
,
J. Chem. Phys.
139
,
104102
(
2013
).
31.
V.
Thapar
and
F. A.
Escobedo
,
Phys. Rev. Lett.
112
,
048301
(
2014
).
32.
F. A.
Escobedo
,
J. Chem. Phys.
140
,
094102
(
2014
).
33.
D. A.
Kofke
,
J. Chem. Phys.
98
,
4149
(
1993
).
34.
M.
Mehta
and
D. A.
Kofke
,
Chem. Eng. Sci.
49
,
2633
(
1994
).
35.
R.
Agrawal
and
D. A.
Kofke
,
Mol. Phys.
85
,
43
(
1995
).
36.
K. S.
Rane
,
S.
Murali
, and
J. R.
Errington
,
J. Chem. Theory Comput.
9
,
2552
(
2013
).
37.
N. B.
Wilding
,
Phys. Rev. E
52
,
602
(
1995
).
38.
N. B.
Wilding
,
Am. J. Phys.
69
,
1147
(
2001
).
39.
F. A.
Escobedo
,
J. Chem. Phys.
115
,
5642
(
2001
).
40.
F. A.
Escobedo
,
Phys. Rev. E
73
,
056701
(
2006
).
41.
A. M.
Ferrenberg
and
R. H.
Swendsen
,
Phys. Rev. Lett.
61
,
2635
(
1988
).
42.
A. M.
Ferrenberg
and
R. H.
Swendsen
,
Phys. Rev. Lett.
63
,
1195
(
1989
).
43.
P. J.
Steinhardt
,
D. R.
Nelson
, and
M.
Ronchetti
,
Phys. Rev. B
28
,
784
(
1983
).
44.
J. R.
Espinosa
,
E.
Sanz
,
C.
Valeriani
, and
C.
Vega
,
J. Chem. Phys.
139
,
144502
(
2013
).
45.
N.
Metropolis
,
A. W.
Rosenbluth
,
M. N.
Rosenbluth
,
A. H.
Teller
, and
E.
Teller
,
J. Chem. Phys.
21
,
1087
(
1953
).
46.
S.
Gottschalk
,
M. C.
Lin
, and
D.
Manocha
, “
OBBTree: A hierarchical structure for rapid interference detection
,” in
Proceedings of the 23rd Annual Conference on Computer Graphics and Interactive Techniques, SIG-GRAPH ‘96
(
Association for Computing Machinery
,
New York
,
1996
), p.
171
.
47.
J.
Arvo
,
Graphics Gems
(
Academic
,
San Diego
,
1990
), pp.
335
339
.
48.
M.
Lamm
and
C. K.
Hall
,
Fluid Phase Equilib.
182
,
37
(
2001
).