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

## I. INTRODUCTION

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 approaches^{8} 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 *NP*_{z}*T* 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 studied^{31}). 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.

## II. METHOD

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.

### A. A generalized ensemble formulation

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*, *f*_{1}, and *f*_{2} 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

where *X*_{i} is the extensive property conjugate to field *f*_{i}, and *Q* is the partition function which is associated with a dimensionless free energy Φ (with *ϕ* being the intensive property) via

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

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

The fundamental thermodynamic equation for this *Mf*_{1}*f*_{2} ensemble is

where

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.

### B. Interfacial pinning

The IP method by Pedersen *et al.* was originally applied in the *NP*_{z}*T* 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 *Mf*_{1}*f*_{2} 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), *f*_{1} and *f*_{2}, two phases with bulk-like quality are present in quantities *M*^{I} and *M*^{II}, 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 *f*_{2}) would be needed to attain coexistence. If the two bulk-like phases are separated by an interface present in quantity *M*_{i}, then the total free energy, Φ, for this two-phase system is given by

where *ϕ*_{i} is the interface contribution to Φ. The IP method assumes the values of *ϕ*_{i} and *M*_{i} to be constant. Using this assumption and the fact that *M* = *M*^{I} + *M*^{II} + *M*_{i}, Eq. (7) reduces to

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 *M*^{I}. 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

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

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:

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

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

where *β* = 1/*k*_{B}*T* (*k*_{B} 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.

### C. Extrapolation of free energies to target coexistence conditions

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:

which by virtue of Eq. (6) reduces to

where |$\bar x_i = \langle {x_i }\rangle$|$x\xafi=\u27e8xi\u27e9$ is the average of *x*_{i}, a quantity readily obtained in the simulation at point B (recall that *X*_{i} = *M x*_{i} 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 *x*_{i} and *X*_{j} : cov_{11}, cov_{22}, and cov_{12}.

Note that cov_{ij} = cov_{ji} 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\xafi=\u222bxiexp(\u2212f1X1\u2212f2X2)dX1dX2/Q(f1,f2)$ it is straightforward to show that

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\xaf1,B,x\xaf2,B$, cov_{11, B}, cov_{22, B}, cov_{12, B}) are known from bulk phase simulations in the ensemble *M f*_{1,B}*f*_{2,B}. Then if *f*_{1} and *f*_{2} are to correspond to phase coexistence so that Δ*f*_{2} = *f*_{2} – *f*_{2,B} and Δ*f*_{1} = *f*_{1} − *f*_{1,B}, by setting ϕ^{I} = ϕ^{II} we can write

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

which is second order in Δ*f*_{2} and hence solvable by the quadratic formula. The difference |$\phi _B^I - \phi _B^{II}$|$\varphi BI\u2212\varphi BII$ in Eq. (20) is obtained from the IP simulation via Eq. (13), while the |$\bar x_{2,B}$|$x\xaf2,B$ and cov_{22, 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

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 *f*_{2} 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\xaf1$ and |$\bar x_2$|$x\xaf2$ for each phases if the correction Δ*f*_{2} [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],

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 functions^{32–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 Δ*f*_{2} typically associated with the needed single-phase simulations that remained (meta) stable, makes Eq. (20) or (21) more suitable and practical.

### D. Examples of ensemble specialization

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

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 *f*_{1} to be a vector **f**_{1} = (β, β*P*), a distinction that has no effect on Eqs. (13) and (20). We have then *M* = *N*, ϕ = βμ_{2}, |$\bar x_1 = v$|$x\xaf1=v$, *f*_{2} = ξ, |$\bar x_2 = y_1$|$x\xaf2=y1$, where *y*_{i} is the mole fraction of component *i*. Equations (13) and (20) then become

It should be noted that we could have chosen the sought-after (*f*_{2}) and specified (*f*_{1}) 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.

### E. Order parameter

To distinguish between isotropic and solid phases, several translational order parameters can be used such as the bond-orientational order parameters *Q*_{6}^{43} 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, *Q*_{6} is defined as

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

where |$Y_{6m} (\vec r)$|$Y6m(r\u20d7)$ are spherical harmonics for the position vector |$\vec r$|$r\u20d7$ and the summation is over all the bonds (total number *N*_{b}) 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 *Q*_{6} 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 *Q*_{6} in terms of generality and it is defined as

where **k** = (2π*n*_{x}/*L*_{X}, 2π*n*_{y}/*L*_{Y}, 0). The values *L*_{X} and *L*_{Y} are the box lengths along the X and Y axes, respectively. The integers *n*_{x} and *n*_{y} 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 *N*^{1/2} for the liquid and solid phase, respectively.

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

where *V* is the volume of the box, *N* is the number of particles, and *V*_{r} 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 *Q*_{6} 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.

### F. Selection of pinning potential parameters *κ* and *a*

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

where *N*_{z} 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 *λ* = *Q*_{6}, |*ρ*_{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 (=*k*_{B}*T*/*κ*) should be significantly smaller than that in the unconstrained *p* (*λ*) for the stable phase (say, |$s_\lambda ^2$|$s\lambda 2$) or else *λ* fluctuations would not be sufficiently suppressed to prevent the disappearance of the metastable phase; hence

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/\beta s\upsilon 2$.

### G. IP algorithm

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

*Mf*_{1}*f*_{2}ensemble on each single-phase system using the same values of the properties*M*,*f*_{1,}and*f*_{2}, where*f*_{2}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\xaf1,x\xaf2, 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

*Mf*_{1}*f*_{2}ensemble with the pinning potential of Eq. (11) to evaluate ⟨*λ*⟩^{′}, and allow estimation of*ϕ*^{I}−*ϕ*^{II}via Eq. (13).Estimate the coexistence value of

*f*_{2}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 *f*_{2} value is attained.

## III. SIMULATION DETAILS AND RESULTS

### A. Single-component polyhedral particles

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 *NP*_{z}*T* ensemble to estimate the coexistence pressure, *P*_{coex}, between isotropic and *rotator* phases. We report our results in terms of reduced coexistence pressure given by |$P_{coex}^*$|$Pcoex*$ = β*P*_{coex}*σ*^{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*. *NP*_{z}*T* 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 × 10^{6} 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 criterion^{45} which requires ruling out overlaps between any two particles (via the separating axes theorem^{46}). For all three shapes, we have applied the method using three different order parameters *Q*_{6},^{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 *n*_{x} and *n*_{y} are chosen to maximize the contrast between liquid and *rotator* phase; in particular we used (*n*_{x}, *n*_{y}) = (22, 0), (9, 9), and (14, 0) for the COs, TOs, and RDs, respectively. The harmonic field parameters *a* and *κ* used for *Q*_{6}, |*ρ*_{k}|, and *υ* for all three shapes are listed in Table I.

Shape . | Order parameter . | a
. | κ
. |
---|---|---|---|

CO | Q_{6} | 0.2 | 80 000 |

υ | 2.06 | 400 000 | |

|ρ_{k}| | 15.0 | 10 | |

TO | Q_{6} | 0.2 | 80 000 |

υ | 2.05 | 400 000 | |

|ρ_{k}| | 25.5 | 10 | |

RD | Q_{6} | 0.2 | 80 000 |

υ | 2.07 | 400 000 | |

|ρ_{k}| | 25.5 | 10 |

Shape . | Order parameter . | a
. | κ
. |
---|---|---|---|

CO | Q_{6} | 0.2 | 80 000 |

υ | 2.06 | 400 000 | |

|ρ_{k}| | 15.0 | 10 | |

TO | Q_{6} | 0.2 | 80 000 |

υ | 2.05 | 400 000 | |

|ρ_{k}| | 25.5 | 10 | |

RD | Q_{6} | 0.2 | 80 000 |

υ | 2.07 | 400 000 | |

|ρ_{k}| | 25.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 *NP*_{z}*T* 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 *Q*_{6} 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 *Q*_{6} 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.

By using both *Q*_{6} 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 *Q*_{6} and *υ* and shows that the converged values of |$P_{coex}^*$|$Pcoex*$ obtained by using the proposed order parameter *υ* matches well with that based in *Q*_{6}. However, a MC cycle in the IP run is approximately 3–4 times longer when using *Q*_{6} 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 *Q*_{6} the cost was minimized by constraining any new calculation to only those particles whose neighbor's positions had been affected by a translation move).

. | . | . | . | Direct . | . |
---|---|---|---|---|---|

. | Order . | . | . | interfacial . | Thermodynamic . |

Shape . | parameter . | |$P_0^*$|$P0*$ . | |$P_{coex}^*$|$Pcoex*$ . | method^{4}
. | integration^{6}
. |

CO | Q_{6} | 2.7 | 2.94 | N/A | 2.9 |

3.2 | 2.94 | ||||

υ | 2.7 | 2.94 | |||

3.2 | 2.95 | ||||

TO | Q_{6} | 2.0 | 2.18 | 2.5 | 2.18 |

2.5 | 2.19 | ||||

υ | 2.0 | 2.19 | |||

2.5 | 2.18 | ||||

RD | Q_{6} | 3.2 | 3.37 | 3.6 | N/A |

3.6 | 3.37 | ||||

υ | 3.2 | 3.38 | |||

3.6 | 3.38 |

. | . | . | . | Direct . | . |
---|---|---|---|---|---|

. | Order . | . | . | interfacial . | Thermodynamic . |

Shape . | parameter . | |$P_0^*$|$P0*$ . | |$P_{coex}^*$|$Pcoex*$ . | method^{4}
. | integration^{6}
. |

CO | Q_{6} | 2.7 | 2.94 | N/A | 2.9 |

3.2 | 2.94 | ||||

υ | 2.7 | 2.94 | |||

3.2 | 2.95 | ||||

TO | Q_{6} | 2.0 | 2.18 | 2.5 | 2.18 |

2.5 | 2.19 | ||||

υ | 2.0 | 2.19 | |||

2.5 | 2.18 | ||||

RD | Q_{6} | 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 scheme^{6} 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.

### B. Mixture of cubes and spheres

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 10^{6} MC cycles for equilibration and 2 × 10^{6} 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 SGNP_{z}T ensemble). All trial moves are accepted according to the Metropolis criterion^{45} which required ruling out overlaps between any two spheres, any cube-sphere pair (via the Arvo's algorithm^{47}) and any two cubes (via the separating axes theorem^{46}).

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*/*σ* = 1^{19} and for *d*/*σ* = 1.23^{32} 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 *Q*_{6} 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).

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} } )$|$\beta (\mu 2solid\u2212\mu 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} } )$|$\beta (\mu 2solid\u2212\mu 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 1^{st} 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 × 10^{5} 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.

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(*y*_{1} − *a*^{′})^{2} where *a*^{′} is chosen to be the average composition (*y*_{1}) 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.

## IV. CONCLUSIONS AND OUTLOOK

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 *NP*_{z}*T* 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}|, *Q*_{6,} and *υ* (specific volume). We find that when using |*ρ*_{k}| the system is unable to stay in a two-phase coexistence state whereas using *Q*_{6} (a more robust translational order parameter) does allow for accurate estimates of |$P_{coex}^*$|$Pcoex*$. The use of *υ*, which unlike *Q*_{6} involves no computational cost to update, also allows for reliable estimates of |$P_{coex}^*$|$Pcoex*$ that match closely those obtained using *Q*_{6} 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 here^{6} 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 GDI^{32} can be seen to rely on the same extrapolation scheme of free energies.

## ACKNOWLEDGMENTS

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.