Mapping coexistence lines via free-energy extrapolation: application to order-disorder phase transitions of hard-core mixtures.

Mapping coexistence lines via free-energy extrapolation: Application to order-disorder phase transitions of hard-core mixtures Fernando A. Escobedo Citation: The Journal of Chemical Physics 140, 094102 (2014); doi: 10.1063/1.4866764 View online: http://dx.doi.org/10.1063/1.4866764 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/9?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Investigation of the vortex order-disorder phase transition line in Bi 2 Sr 2 Ca Cu 2 O 8 + x using ac techniques J. Appl. Phys. 103, 07C705 (2008); 10.1063/1.2833818 Surface order-disorder phase transitions and percolation J. Chem. Phys. 125, 184707 (2006); 10.1063/1.2370875 Monte Carlo simulation and phase behavior of nonadditive hard-core mixtures in two dimensions J. Chem. Phys. 117, 5780 (2002); 10.1063/1.1501126 Gas–liquid nucleation in partially miscible systems: Free-energy surfaces and structures of nuclei from density functional calculations J. Chem. Phys. 111, 5485 (1999); 10.1063/1.479809 A comprehensive study of the phase diagram of symmetrical hard-core Yukawa mixtures J. Chem. Phys. 109, 4498 (1998); 10.1063/1.477053

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.120.242.61 On: Mon, 01 Dec 2014 21:58:08

THE JOURNAL OF CHEMICAL PHYSICS 140, 094102 (2014)

Mapping coexistence lines via free-energy extrapolation: Application to order-disorder phase transitions of hard-core mixtures Fernando A. Escobedoa) School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, New York 14853, USA

(Received 30 December 2013; accepted 9 February 2014; published online 6 March 2014) In this work, a variant of the Gibbs-Duhem integration (GDI) method is proposed to trace phase coexistence lines that combines some of the advantages of the original GDI methods such as robustness in handling large system sizes, with the ability of histogram-based methods (but without using histograms) to estimate free-energies and hence avoid the need of on-the-fly corrector schemes. This is done by fitting to an appropriate polynomial function not the coexistence curve itself (as in GDI schemes) but the underlying free-energy function of each phase. The availability of a free-energy model allows the post-processing of the simulated data to obtain improved estimates of the coexistence line. The proposed method is used to elucidate the phase behavior for two non-trivial hard-core mixtures: a binary blend of spheres and cubes and a system of size-polydisperse cubes. The relative size of the spheres and cubes in the first mixture is chosen such that the resulting eutectic pressure-composition phase diagram is nearly symmetric in that the maximum solubility of cubes in the sphere-rich solid (∼20%) is comparable to the maximum solubility of spheres in the cuberich solid. In the polydisperse cube system, the solid-liquid coexistence line is mapped out for an imposed Gaussian activity distribution, which produces near-Gaussian particle-size distributions in each phase. A terminal polydispersity of 11.3% is found, beyond which the cubic solid phase would not be stable, and near which significant size fractionation between the solid and isotropic phases is predicted. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4866764] I. INTRODUCTION

Mapping phase diagrams of model materials via molecular simulation is one of the main scientific contributions that simulators have made during the last decades, efforts that have been accompanied by the development of many distinct approaches.1–33 In most of these methods, simulations are performed at specific thermodynamic states in purposely selected ensembles and can be roughly grouped as methods that allow: (i) direct determination of the coexistence conditions for both phases, (ii) tracing points along a coexistence line in a step wise manner, or (iii) connecting reference states (of known free energy) to the conditions of interest. Examples for type (i) methods are interfacial simulations and the Gibbs-ensemble method,6–12 for type (ii) methods are Grand canonical and other ensembles whose results are extrapolated via multi-histogram reweighting (MHR),12–19 or via GibbsDuhem integration (GDI),34, 35 and for type (iii) methods are various versions of thermodynamic integration.1–5 A different class of methods, type (iv), can be denoted as “density of states” methods since they aim at obtaining the underlying entropy function (or part thereof) for the multiplicity of macrostates defined by extensive properties (such as energy, volume, etc.), without necessarily simulating specific state points.20–33 Often, different methods are used to provide complementary information. While simulating phase diagrams for some systems is relatively straightforward via different

a) [email protected]

0021-9606/2014/140(9)/094102/17/$30.00

conventional methods, e.g., involving fluid phases of small molecular species, many systems still present challenges. In particular, systems that involve order-disorder phase transitions for single and specially multiple components, typically require numerous simulations, specialized implementations and significant care to map out coexistence lines efficiently (e.g., type-4 methods have shown limited applicability for such systems so far). In those situations, thermodynamic integration or the GDI method34–41 have been widely used and are typically the methods of choice. The GDI method relies on the pre-knowledge of phaseequilibrium for at least one point along the sought-after coexistence line. Those conditions provide the “initial value” to integrate a suitable Gibbs-Duhem differential equation (e.g., a Clapeyron-type of equation for a pure component). Since advancing the integration one step forward is essentially an extrapolation operation, this can be done using other methods, e.g., via the MHR method. A key connection between GDI methods and MHR methods was examined in a series of papers and cast within the framework of pseudo-ensembles.42–44 A pseudo-ensemble45–47 can be seen as embodying the generalized idea of extrapolation, and that in some situations it may be more convenient to use a “surrogate” ensemble in order to mimic a “target” ensemble (rather than simulating the target ensemble). Since an ensemble is merely a specification of thermodynamic conditions, in the context of GDI the “target” ensemble represents the coexistence conditions one integration-step forward, obtained by extrapolation of results from the current-state “surrogate” ensemble. Pseudo ensembles also provided a unifying framework to cast other methods

140, 094102-1

© 2014 AIP Publishing LLC


094102-2

Fernando A. Escobedo

that have used single-point extrapolations as part of the process to target coexistence states.48–51 In some applications of the MHR approach, one does not map the coexistence line in a stepwise manner but simulates many points (at different conditions) at once and then combines the results from all of them to calculate the coexistence conditions as a post-processing step (making first sure that results match a known phase-equilibrium point or that at least one simulation has been able to ergodically sample configurations from both phases, e.g., near a vapor-liquid critical point or aided by multicanonical or “flat” histogram methods12–14, 51–54 ). The chosen points for such simulations need not be at coexistence but for the method to work, they should be sufficiently near coexistence (to sample the relevant parts of phase space to allow extrapolation to coexistence conditions) and near to each other (at least in pairs) to connect overlapping regions of phase space. Such schemes could then be seen more properly as “refinement” methods if not used to generate the set of near-coexistence points to begin with; some reported applications take the original set of points from former publications, or from preliminary runs which may exploit an approximate functional form of the coexistence line12 (like the Clausius-Clapeyron linear relation between pressure and reciprocal temperature). In this context, this work is more interested in methods that can do both, generating points near a coexistence line and the more straightforward task of refining the initial estimates. For this reason, we only focus on methods that map out a coexistence line in a sequential pointwise manner, with every new point resulting from the extrapolation (via GDI or MHR) of previous points. One of the shortcomings of the original GDI implementations is that they were sensitive to the accumulation of errors arising not only from the integration scheme adopted but also from the simulated averages that always needed to correspond to coexistence conditions: each prediction step was complemented by an on-the-fly corrector scheme. The coexistenceline free-energy difference integration (CFDI) method of Meijer and Azhar55 got around this issue by only implementing predictor steps where the simulation points need only be near coexistence; coexistence conditions could be estimated (without requiring a new run) as a post-simulation step by obtaining first the free-energy differences between successive pairs of points (on each phase) using Bennett’s method.56 The CFDI is in fact isomorphic to a method where the postprocessing extrapolation to coexistence is done via a MHR method that uses only the two nearest simulation points (for each phase) at a time.57 Indeed, the “weighting” of different simulation points associated with MHR methods requires the estimation of free energy differences between points. In this context, the original GDI method has one key advantage over the CFDI or MHR methods: it is less sensitive to system size. It is a known issue of free-energy difference and histogrambased methods that the maximum “step-size” between simulation points to allow proper estimation of free-energy differences (related to the extent of histogram overlap) decreases substantially with system size. Note that the need of overlapping histograms is also crucial, albeit implicit, for effective operation of other methods such as replica-exchange53, 57 (used sometimes in conjunction with methods to map coex-

J. Chem. Phys. 140, 094102 (2014)

istence lines). Van’t Hof et al.58 proposed a method where multiple simulations are performed (and the results combined via MHR) to advance a typical GDI integration step; it gets around the potential problem in CFDI of non-overlapping histograms between pairs of points, but at the cost of several intermediate additional runs. Furthermore, larger sizes also imply harder convergence and noisier histogram “tails;” in such a case, one can only expect the first moments and at most 2nd moments of the underlying distribution functions (histograms) to be statistically meaningful. GDI is more robust in that respect because, whatever integration scheme one adopts, it essentially amounts to assuming that a smooth function, typically a polynomial constructed with first-moment data, fits a portion of the coexistence line (a relation between the saturation properties) being mapped out. In this work we propose a variant of the GDI method that combines some of the advantages outlined above of the original GDI methods (robustness in handling large system sizes), with the ability of CFDI and MHR to estimate freeenergies and hence avoid the GDI need of on-the-fly corrector schemes to simulate at the exact coexistence conditions. This is done by fitting to a smooth polynomial function not the coexistence curve itself (as in GDI schemes) but the underlying free-energy function of each phase (i.e., the logarithm of the ensemble partition function). This method, which we term Free-Energy Extrapolation “FENEX” method, is shown to naturally reduce to some limiting cases of MHR, CFDI, and GDI. FENEX then addresses the issue of both finding an approximate set of coexistence conditions and refining it in the same spirit of the original GDI methods but, in the spirit of MHR, without additional simulations for on-the-fly corrections. The method is applied to the mapping of pressurecomposition phase diagrams for two non-trivial systems: (1) a binary mixture of hard spheres and cubes and (2) a polydisperse mixture of hard cubes. Elucidating the phase behavior of colloidal suspensions containing polyhedral particles like cubes is of both theoretical and practical importance.59–64 On the one hand, cubes form an unusual ordered mesophase61–63 whose tolerance to guest components or to size polydispersity is still not fully understood;63 on the other hand, the growing ability to synthesize cubes of any desirable (nano)size and surface chemistry has fueled interest in preparing suspensions and layered assemblies for various potential uses. II. THEORY A. Basic definitions

Consider a system where the ensemble properties M, f1 and f2 are fixed. For simplicity we only consider two thermodynamic fields “f” although the generalization to more variables is straightforward. Assume that the probability of any macrostate in the selected ensemble can be written as (X1 , X2 |f1 , f2 ) = W (X1 , X2 ) exp(−f1 X1 − f2 X2 )/Q(f1 , f2 ), (1) where W is thedensity of states and the partition function is Q(f1 , f2 ) =

W (X1 , X2 ) exp (−f1 X1 − f2 X2 ) dX1 dX2 , (2)


094102-3


J. Chem. Phys. 140, 094102 (2014)

and where Xi is the extensive property conjugate to field fi . If the associated free energy (φ being the intensive property and = φM the extensive property) is related to the partition function via = φM = − ln Q(f1 , f2 ).

taking logarithms we get dX1 dX2 (X1 , X2 |f1,B , f2,B ) ln(Q/QB ) = ln × exp(−X1 f1 − X2 f2 ) .

(3)

Then the fundamental thermodynamic equation for this M f1 f2 ensemble is d = X¯ 1 df1 + X¯ 2 df2 ,

(4)

dφ = x¯1 df1 + x¯2 df2 ,

(5)

x¯i = X¯ i /M = (∂φ/∂fi )fj ,j =i .

(6)

where

It is assumed that at any simulated point one is able to obtain the first and second moments of the main properties of interest: x¯1 , x¯2 , cov11 , cov22 , cov12 , where x¯i = xi is the average of xi and covij = covji is the covariance of variables xi and Xj :

(11)

If the probability density obtained at point B can be approximated via a bivariate Gaussian distribution, i.e. (X1 , X2 |f1,B , f2,B ) ≈

1 1 / 2 2π covX11 ,B covX22 ,B θ 1 (X1 − X¯ 1 )2 (X2 − X¯ 2 )2 × exp − + 2θ covX11 ,B covX22 ,B

(X1 − X¯ 1 )(X2 − X¯ 2 ) , (12) −2 covX11 ,B covX22 ,B / covX12 ,B

(7) covij = xi Xj − xi Xj = Xi xj − Xi xj , from x¯i = W (X1 , X2 )xi exp (−f1 X1 − f2 X2 ) dX1 dX2 /Q (f1 , f2 ) it is straightforward to show that

where θ = 1 − cov2X12 /(covX11 covX22 ) and covXij = Xi Xj − Xi Xj , then the integral on the right-hand side of Eq. (11) is essentially a definition of the characteristic function (or moment generating function) of a bivariate Normal distribution and has a known explicit analytical solution.65 Using that result and that the left-hand side of Eq. (11) is [according to Eq. (3)] just −M(φ−φ A ) leads to

−covij = (∂ x¯i /∂fj )fk ,k=j = ∂ 2 φ/∂fi ∂fj .

−M(φ − φB ) = −X¯ 1,B f1 − X¯ 2,B f2 + 1/2covX11 ,B (f1 )2

(8)

For such an ensemble, if two phases I and II (used therefrom as superscripts) are at coexistence, then it must be true that not only f1I = f1I I and f2I = f2I I but also that φ I = φ II . To give a concrete example for this notation, in the isothermal-isobaric NPT ensemble one would have that M = N (number of particles), f1 = β(= 1/kT, where k is Boltzmann’s constant), f2 = βP (P = pressure), and φ = βμ (where μ is chemical potential), so that x¯1 = u(molar energy), 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. In implementing the method for other ensembles, variables φ, fi , and x¯i can be readily identified by comparing to Eq. (2) the particular ensemble partition function, or by comparing to Eq. (5) a suitable version of the Gibbs-Duhem equation: yi d(βμi ) = udβ + vd(βP ), (9) i

where yi is the mole fraction of component i. B. A second order Taylor expansion of φ is equivalent to assuming that is a Gaussian distribution

Applying Eq. (1) at conditions ( f1 , f2 ) and ( f1,B , f2,B ) and eliminating from these the common function W(X1 ,X2 ): Q(X1 , X2 |f1 , f2 ) = QB (X1 , X2 |f1,B , f2,B ) exp (−X1 f1 −X2 f2 ),

(10)

where for brevity we wrote f1 = f1 − f1,B , f2 = f2 − f2,B , Q = Q( f1 ,f2 ), and QB = Q( f1,B ,f2,B ). Integrating both sides and

+ 1/2covX22 ,B (f2 )2 + covX12 ,B f1 f2 . Or since covXij = Mcovij : φ − φB = x¯1,B f1 + x¯2,B f2 − 1/2cov11,B (f1 )2 − 1/2cov22,B (f2 )2 − cov12,B f1 f2 . (13) On account of the definitions of Eqs. (6) and (8), it is clear that Eq. (13) is the 2nd order Taylor expansion of φ around point B. By partial differentiation of Eqs. (13) and (5), it is easy to see that in this case, x¯1 = x¯1,A − cov11,B f1 − cov12,B f2 ,

(14)

x¯2 = x¯2,B − cov22,B f2 − cov12,B f1 .

(15)

This connection between a 2nd order Taylor expansion of φ [Eq. (13)] and a Gaussian distribution for the probability density of macrostates [Eq. (12)] was already pointed out in Ref. 42 but not explicitly shown. C. Generalized polynomial extrapolation to phase coexistence

If Eq. (13) is applied to each of two phases I and II so that properties from each phase at point B (x¯1,B , x¯2,B , cov11,B , cov22,B, cov12,B ) are known from simulations in the ensembles MI f1,B f2,B (for phase I) and MII f1,B f2,B (for phase II) then if f1 and f2 are to correspond to phase coexistence, by setting φ II = φ I we can write


094102-4


J. Chem. Phys. 140, 094102 (2014)

II I I f1 − 1/2 covI11,B − x¯1,B − covI11,B (f1 )2 φBI I − φBI + x¯1,B II f2 = I I . I I + cov12,B − covI12,B f1 + 1/2 covI22,B − x¯2,B − x¯2,B − covI22,B f2

This equation allows predicting the coexisting f2 value for a specified f1 if the difference φBI I − φBI is known (e.g., if point B was previously determined to be at coexistence so that φBI = φBI I ). Note that Eq. (16) is quadratic in f2 and hence solvable via the quadratic formula. If point B was already at coexistence, then in the limit when f1 → 0 and f2 → 0, Eq. (16) reduces to the conventional Clapeyron-type of differential equation which is the usual starting point for GDI:35

df2 df1

coex

II I x¯1,B − x¯1,B . = − II I x¯2,B − x¯2,B

(17)

This underscores the idea that Eq. (16) can be seen as a recipe for stepwise integration of Eq. (17). Now, the analysis of Sec. II B demonstrates that the quadratic expansion of Eq. (13) is a suitable approximation to the local behavior of φ ( f1 , f2 ) (and hence of the partition function) provided that the actual histogram at state point B approaches a Gaussian distribution (i.e., for large systems). The fact that the second-moment covij data may be nosier and hence less reliable than the first moments x¯i data, suggests that the coefficients of a 2nd order expansion of φ could also be obtained instead by fitting data from more than a single state [in contrast to Eqs. (13) and (16)]. Hence one could have φ = c10 f1 + c01 f2 + c20 (f1 )2 + c02 (f2 )2 + c11 f1 f2 ,

(18)

where the c’s are constant coefficients, which could be found by, e.g., fitting exactly the x¯i data from the most recent point (B) and minimizing the deviation of covij data from that point and of x¯i data from a previous point. Further, a change in variables may be indicated if it were known that the transformations g1 = g1 (f1 ),

g2 = g2 (f2 )

(19)

allow for a smoother, slower varying φ (g1 , g2 ) function, potentially leading to more accurate low-order expansions on g1 and g2 , or to capture (or resolve) a non-analytic/divergent behavior of φ or fi at some limiting conditions. Examples of these situations may occur in using log P as opposed to P (as in Clapeyron’s equation) or the fugacity in lieu of a chemical potential. Note that since the (user specified) g functions of Eq. (19) are known, so are their inverse functions: f1 = f1 (g1 ),

f2 = f2 (g2 ),

(20)

and the first and second derivatives of φ ( f1 , f2 ) found in simulations [Eqs. (6) and (8)] are readily relatable to derivatives of φ (g1 , g2 ) [see formula (A3) in Appendix A]. Of course, one could not expect a 2nd order polynomial function to be valid over a wide range of ( f1 , f2 ) or (g1 , g2 ) away from ( f1,B , f2,B ), but our analysis makes it plausible that

(16)

polynomial functions are suitable choices to describe the behavior of φ. Generally, if additional information (about the first and second moments of ) is known for several points, it is then justifiable to use a higher order polynomial expansion of φ, such that its coefficients can be found by best-fitting procedures. To account for a generalized situation, we will assume that for each phase J we have found suitable coefficients for a polynomial model, φJ =

p p−n

J cmn g1m g2n ,

(21)

n=0 m=0 J where p is a preset integer, c00 = φBJ , and gi = gi – gi,B , so that the phase-equilibrium equation equivalent to Eq. (16) is m p I I I g1 − cm0 φBI I − φBI + m=1 cm0 g2 = p p−n , (22) m n−1 I I I − n=1 m=0 cmn − cmn g1 g2

which is a polynomial of order p in the unknown g2 and hence analytically solvable if p ≤ 4; alternatively, the solution can often be obtained by successively iterating Eq. (22), starting with an initial guess of g2 = 0 on the right-hand side. With the new value of g2 , the sought after value of f2 is obtained from f2 = f2 (g2 ) [Eq. (20)]. To find the fitting constants c in Eq. (21) one needs to first find the numerical values of the first and second derivatives of φ (g1 ,g2 ) based on those for φ ( f1 , f2 ) (i.e., in terms of the x¯i and covij measured in simulations). Details of different choices of p and the resulting formulas to evaluate the c constants are given in Appendix A. Note that while a larger value of p may seem more appealing as the associated polynomial can incorporate more simulated data, lower orders are probably preferable for forward extrapolation, to avoid the spurious divergence that high order terms can lead to, or to better capture a rapidly changing φ near the latest point. Regardless of the particular φ model adopted, the key idea of FENEX, as illustrated in Fig. 1(a) is to find local fits to the φ J ( f1 , f2 ) surface for the two phases (I and II) that resolve the two-phase coexistence region.

D. FENEX recipe to map out coexistence lines

Given simulations in the Mf1 f2 ensemble, the goal is to map the phase diagram in the f1 − f2 plane by stepping over prescribed values of f1 so that the corresponding coexistence values of f2 are to be calculated. The two phases are denoted by I and II and is assumed that the mapping will be performed for fixed MI and MII . It is assumed that for any simulation one conducts equal-length runs to obtain average values for the properties x¯1 , x¯2 , cov11 , cov22 , cov12 (and associated errorbars) on each phase.


094102-5


J. Chem. Phys. 140, 094102 (2014)

(1) At the start up one assumes that at point “B” one knows the values of f1,B and f2,B that lead to phase coexistence or near coexistence of phases I and II (and hence that φBI = φBI I or that the difference φBI I − φBI is known), at which conditions simulations were conducted for both phases so that x¯1,B , x¯2,B , cov11,B , cov22,B , cov12,B are known for each phase. Then if f1 = f1 – f1,B is given, the corresponding value of g1 = g1 – g1,B is obtained [Eq. (19)] and the value of g2 = g2 – g2,B for which coexistence is attained is estimated from Eq. (22), noting that, given the limited data available, the highest order polynomial for Eq. (21) that can be used is p = 2 [i.e., Eq. (16) if gi = fi ]. However, if covij,B values are not known to sufficient accuracy, one may adopt p = 1 so that the usual linear extrapolation formula is used [Eq. (17)]. (2) Re-label point B (and all properties thereof) as point A and perform simulations at conditions for the new point “B” ( f1,B = f1 , f2,B = f2 ) for the two phases I and II to obtain the corresponding values for x¯1,B , x¯2,B , cov11,B , cov22,B , cov12,B . (3) Calculate the free energy differences between the last two states simulated for each phase and hence estimate φBI and φBI I . With 1 = g1,A − g1,B , 2 = g2,A − g2,B , and c00 = 0 we get J m n cmn 1 2 , − φBJ − φAJ = p

p−n

(23)

n=0 m=0

FIG. 1. (a) Depiction of free-energy surfaces (φ) for two phases (in blue and red) and how a FENEX step in effect generates pieces of their surfaces around points A and B (patches shown enclosed by dotted and dashed lines) which allow extrapolation toward nearby coexistence conditions (i.e., point C at which φCI = φCI I ). (b) and (c) illustrate the free-energy extrapolation for a given phase based on two simulation points A and B showing, for simplicity, only single-variable functions [in reality they should be plotted as 3D surfaces on 2 variables: X1 and X2 in (b) and f1 and f2 in (c)]. Panel (b) shows that the probability density distributions need not overlap. In the φ vs. fi curve in (c), simulations do not directly find φ A or φ B , but they provide data for the slope (∂φ/∂fi ) and curvature (∂ 2 φ/∂fi 2 ) at those points; these data are depicted by short thick segments. For point A it is assumed that φ A has already been set, but the position of point B is indeterminate as illustrated by the green segment that can be arbitrarily shifted up or down. A polynomial function fitting all data at A and B (including slopes and curvatures) produces the dashed curve, which allows both estimating φ B and extrapolating toward a new point C. For these plots we assumed that g = f.

which is applied separately to each phase J (= I and II). (4) At this stage all relevant properties of at least the previous two points A and B and for both phases I and II are known, which are used to fit the c constants of model (21). Then if f1 = f1 – f1,B (and g1 = g1 –g1,B ) is given for a new point, the value of g2 and hence f2 = f2 (g2 ) for which coexistence (φ II = φ I ) is expected is obtained from Eq. (22). At this stage one goes back to step 2 and iterate steps 2–4 until the target number of steps has been completed. (5) As a post processing task, coexistence properties can be refined from the “near-coexistence” simulated data. Indeed, at any given point B, the simulations were performed at the same values of ( f1,B , f2,B ) in each phase that were “preliminary” estimates of coexistence conditions; improved estimates (marked with an asterisk) can be obtained by extrapolating the results obtained at ( f1,B , f2,B ) using a variant (with g1 = 0) of Eq. (22):

∗ g2,B − g2,B = g2∗ =

−

p n=1

Then x¯1∗,J = (dg1 /df1 )f1,B ∗ x¯2∗,J = (dg2 /df2 )f2,B

φBI I − φBI . (24) II I c0n (g2∗ )n−1 − c0n

p−1

J c1n (g2∗ )n ,

(25)

J nc0n (g2∗ )n−1 ,

(26)

n=0

p n=1

∗ ∗ = f2 (g2,B ) [Eq. (20)]. The last two relations corwhere f2,B respond to ∂φ/∂fi evaluated from Eq. (21) for g1 = 0 and


094102-6


g2 ∗ , and are applied to each phase (J = I and II). Note that the polynomial fits used for Eqs. (24)–(26) could be the same as those used when originally marching from point B forward, or be ones that result from combining more or even all simulation data (as in MHR). On the other hand, if g2∗ is very small and σ ij,B data are accurate, then one could simply use first order corrections as in Eqs. (14) and (15). E. Remarks on the FENEX scheme

Note that one could in principle use different polynomial fits for different operations. For example, the polynomial used for extrapolating forward from point B (i.e., for steps 1 and 4) could be more heavily weighed by the data at the most recent point than those from previous (most distant) points, while the polynomial used to estimate the free-energy difference between two points (in step 3) could weigh equally data from both such points. For step 3, Bennett’s method (which gives an identical prescription to MHR with two histograms) could be applied so as to numerically solve for φ (for each phase J) the equation J (X1 , X2 |f1,A , f2,A )dX1 dX2 1 + exp(−X1 1 − X2 2 − M J φ J ) J (X1 , X2 |f1,B , f2,B )dX1 dX2 , (27) = 1 + exp X1 1 + X2 2 + M J φ J where φ J = φBJ − φAJ and 1 = f1,A − f1,B , 2 = f2,A − f2,B . Application of Bennett’s method with the actual histograms (as with the CFDI and MHR methods) would require these to have significant overlap to produce good estimates of free energy differences φ J . In contrast, the simple Eq. (23) could still produce sensible results even if such histograms do not overlap by exploiting the additional conditions imposed by the continuity and smoothness of the assumed φ J ( f1 , f2 ) polynomial function. This idea is illustrated in Figs. 1(b) and 1(c) for a one dimensional function φ( fi ). Note that in the limit of large system size the functions would tend to delta functions, in which case Eq. (27) reduces J J J J + x¯1,B )1 + 1/2(x¯2,A + x¯2,B )2 which is to −φ J = 1/2(x¯1,A the same result that Eq. (23) would give for a model with p = 2, gi = fi that fits exactly first-derivative data at points A and B. The FENEX scheme is readily generalizable to cases where more variables enter into the ensemble formulation (e.g., for multiple components). Consider for concreteness the case when M and 3 field properties f0 , f1 , and f2 define the ensemble [e.g., the “M f0 f1 f2 ” ensemble with the fundamental equations: dφ = x¯0 df0 + x¯1 df1 + x¯2 df2 and φM = −ln Q( f0 , f1 , f2 ) in lieu of Eqs. (5) and (3)]. As long as one is mapping the phase diagram in the f1 − f2 plane (where values of f1 are preset and f2 values are calculated) while f0 is always constant, then scheme above would apply with little modification. III. MODELS AND SIMULATION DETAILS

The FENEX method was first implemented in an NPT ensemble [with notation described in paragraph after Eq. (8)]

J. Chem. Phys. 140, 094102 (2014)

to validate our implementations by reproducing the known vapor-liquid coexistence line of the pure van der Waals fluid.34 For the mixture systems of interest, however, simulations used semigrand-type of ensembles only.

A. Isothermal, isobaric, semigrand ensemble (SGNPT)

Both a binary mixture of hard spheres + hard cubes and a size-polydisperse mixture of cubes were simulated in suitable versions of the SGNPT ensemble.36, 66, 67 In this ensemble, the total number of particles N is fixed at N = 1000 (except for the sphere-rich fcc phases for which N = 864) and the temperature is fixed at a reduced value of T∗ = (βε)−1 = 1 (where ε is an arbitrary parameter), noting that T is irrelevant in the hard-core systems being simulated. The composition of the system is not fixed but allowed to fluctuate in response to specified values of the chemical potential differences βμi = β(μi − μr ) for each species i, where r is a reference component. The dimensionless osmotic pressure is also specified in the ensemble, given in reduced units as P∗ = βPl3 , where l = 1 is the characteristic length of a reference component in the simulated systems. At each state point in the integration, both phases (in putative coexistence) are simulated using a cubic simulation box. Runs involved at least 106 MC cycles for equilibration and 3 × 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 4 volume move attempts. Volume moves for the polydisperse cube system involved the coupled rescaling of particle sizes and box length.67 Swap moves involved picking randomly two particles of different species and attempting to swap their positions (while keeping the orientation for the cubes). All trial moves are accepted according to the Metropolis criterion57 which in our cases required ruling out overlaps between any two spheres, any cube-sphere pair (via the Arvo’s algorithm68 ) and any two cubes (via the separating axes theorem69 ). 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). For the cubes + spheres system the two components are of similar size so that the acceptance rate of mutations was in the range 1%–20% (lower at higher densities). For the polydisperse cubes, a new trial size is chosen by increasing or decreasing the original size by a random amount (with a preset maximum value to achieve a 30% acceptance rate); note that in this case there is an unbound number of components as sizes occur over a continuous range. Assuming that no overlap is incurred in a trial mutation, the acceptance probability if given by Pacc (old → new)

μi,new −β μi,old , = min 1, exp β i

i

(28)

where the sums are over all particles present in the new or old states. The size of translation, rotation, and volume perturbations was adjusted so as to get an acceptance probability of 0.4, 0.4, and 0.3, respectively.


094102-7


J. Chem. Phys. 140, 094102 (2014)

B. Mixtures of cubes and spheres

In this case we have two components 1 and 2 corresponding to the cubes (of edge length σ ) and spheres (of diameter d), respectively, whose difference in chemical potential needs to be specified. The fundamental thermodynamic equation for this SGNPT ensemble can be obtained by rearranging Eq. (9) setting the configurational energy u = 0 and letting ξ = β(μ2 − μ1 ) to obtain d(βμ2 ) = vd(βP ) + y1 dξ,

(29)

where y1 is the mole fraction of component 1. The formal connection with the notation of Sec. II is established by comparing Eqs. (29) and (5) letting M = N, φ = βμ2 , and either: (a) f1 = βP, x¯1 = v and f2 = ξ , x¯2 = y1 or (b) f1 = ξ , x¯1 = y1 and f2 = βP, x¯2 = v. Choice (a) is indicated if steps in βP are to be prescribed, while choice (b) applies when steps in ξ are to be prescribed. We adopted choice (a) for this study and in reduced pressure units P∗ = βPσ 3 , the step size was P∗ = 0.4 for most points (with smaller values near the purecomponent states). In total 48 state points were simulated (for all phases involved), each point requiring from 3 to 7 h in a R 5500, 3.06 GHz processor node. Xeon As validation for the FENEX predictions, selected interfacial simulations were conducted for this system. In interfacial simulations the two coexisting phases and the intervening interfaces are present in the same box (much elongated along one direction to facilitate the formation of distinct bulk regions) at a given pressure, composition, and total number of particles. Details of these simulations are given in Refs. 63 and 64. C. Size polydisperse mixture of hard cubes

In this ensemble we fix the temperature, the pressure, total number of particles (N), and the spectrum of differences in chemical potentials with respect to a virtual species (say 1). We rearrange the fundamental thermodynamic equation (9) for this case to be (for u = 0), d(βμ1 ) = vd(βP ) − yi dξi , (30) i

where as before ξ i = βμi − βμ1 . We follow Ref. 67 and choose a Gaussian spectrum of fugacities so that, ξi = −(σi − σ1 )2 /2υ,

(31)

where σ is the edge length of particle i and υ is the single variable controlling the width of the distribution (the monodisperse case would be recovered for υ → 0). Equation (30) can then be rewritten as d(βμ1 ) = vd(βP ) + m2 d(1/υ). Provided that we define m2 =

1 2

(32)

immaterial to produce the βP vs. 1/υ (or υ) coexistence curve but one must ensure that the slope (df2 /df1 )coex of the line being traced remains finite. In our case, most steps were in βP [choice (a)], which was particularly suitable for larger polydispersity since the coexistence (1/υ) was a non-monotonic function of βP. In the reduced pressure units P ∗ = βP σ13 , step sizes were not constant because P∗ tended to diverge as polydispersity approached a limiting value67 (see Sec. IV B); instead pressure values were increased by a numerical factor that started at 1.05 (near the monodisperse limit) and gradually increased to 1.3 for larger pressures. In total 64 state points were simulated (for both phases), each point requiring R 5500, 3.06 GHz processor node. from 2 to 7 h in a Xeon It is important to point out that imposing an activity distribution as with Eq. (31) will not capture the experimentally relevant phase behavior encountered when fixing the composition of the parent distribution or of one of the phases (e.g., in generating cloud or shadow curves). The limitations of the approach adopted here are well known,70–73 and while alternative methods to simulate polydisperse behavior exists,72–76 these typically require iterative schemes that are not readily cast into the FENEX approach. D. FENEX formulas

In applying FENEX, several of the polynomial fits described in Appendix A were tested and found to give consistent results. For both systems studied and for the simulation moves, lengths, and parameters used, we found that second moment data were accurate enough so that exact-fit (as opposed to best-fit) polynomials were usable. Except at the starting state point in the process (see below), we set gi = fi and employed the constants from formulas (A11) and (A18), and (A19) for marching forward (steps 1 and 4 in Sec. II D) and for extrapolating to coexistence [post processing Eqs. (24)– (26) in step 5 of Sec. II D], while formula (A26) was used for evaluating the free-energy differences (step 3 in Sec. II D). Coexistence line mapping was started at the known twophase coexisting conditions for the pure components. For the size-polydisperse cubes, we followed Ref. 67 to obtain the initial slope needed to start the extrapolation (over g1 = 1/f1 = υ) at the known conditions for the monodisperse cubes (υ → 0). For the sphere + cube system, specialized formulas are needed to start up the process at the infinite-dilution conditions of one component; these are derived in Appendix B. The two order-disorder coexistence lines, one started from the pure spheres and the other from pure cubes, met at the eutectic point (see Sec. IV A). IV. RESULTS A. Binary mixture of cubes and spheres

yi (σi − σ1 ) . 2

(33)

i

Then the formal connection with the notation of Sec. II [Eq. (5)] is established by setting M = N, φ = βμ1 and either: (a) f1 = βP, x¯1 = v, f2 = 1/υ, x¯2 = m2 or (b) f1 = 1/υ, x¯1 = m2 , f2 = βP, x¯2 = v. The choice of (a) or (b) is often

Hard spheres are known to exhibit a first-order phase transition from an isotropic phase at low pressures to a crystalline fcc phase at high pressures and an order-disorder transition pressure “ODP” of βPd3 = 11.57, where d is the sphere diameter.57 Monodisperse cubes59–63 exhibit a transition between an isotropic phase and a cubic mesophase at an ODP


094102-8


near βPσ 3 ≈ 6.2, where σ is the cube edge length; in such a mesophase the particle positions lie on average in a cubic lattice62 but exhibit liquid-like fluctuations and mobility.61–63 The phase behavior of spheres and cubes for the case that the diameter of the former is equal to the side length of the latter (i.e., d = σ ) has been mapped out earlier via interfacial simulations.64 It was found that such a mixture exhibits the expected eutectic type of behavior (seen in hard-sphere mixtures for some range of size disparity40, 77–79 ) at pressures in between the ODP for the pure components, with minimal inter-species mixing at high pressures where sphere-rich and cube-rich crystal phases ensue. The eutectic composition was near the pure hard sphere limit, and while the cube-rich cubic phase was able to dissolve a significant amount of spheres at the eutectic pressure (up to 20% on a mole basis), the sphererich crystal was almost pure, unable to dissolve cubes.64 This asymmetry can be explained by the larger volume of the cubes relative to spheres, which makes it impossible to accommodate cubes into the spheres’ fcc-lattice sites. While one could

FIG. 2. Results for mixture of cubes + spheres for case when pure components have equal ODP (d = 1.23 σ ): (a) Pressure-composition phase diagram where results from FENEX are shown by filled circles and those from interfacial runs by open squares. (b) Pressure vs. volume fraction for the different phases in coexisting with another phase (shown within parentheses). S = sphere-rich fcc solid, I = isotropic liquid , and C = cube-rich cubic ordered phase.

J. Chem. Phys. 140, 094102 (2014)

infer that making the sizes of both species more “similar” would enhance mutual solubility, it is non-trivial to determine what size ratio of components would be optimal. It has been hypothesized that the disparity of pure-component ODP values is a key determinant of ordered-phase compatibility.64 The ODP can be seen as marking the turning point where packing entropy takes over as the dominant entropic force determining the structure of the system; i.e., the ODP value captures the readiness of a pure component to order. We hence chose d/σ = 1.23 so that the both pure-component ODPs match. Using the FENEX method, the resulting phase diagram presented in Fig. 2 shows not only a more symmetric eutectic behavior (with a eutectic composition of ∼42% cubes) but also a sphere-rich fcc crystal that can now dissolve a large amount of cubes (up to 20%) without sacrificing the solubility of spheres in the cube-rich solid seen when d/σ = 1.64 Sample snapshots are shown in Fig. 3, where one can see that the cubes in the sphere-rich fcc phase are not orientationally ordered. Results for other particle size ratios and for other mixtures would be needed to see if indeed the equality of pure-component ODPs has merit as a “symmetrizing” or “compatibilizing” rule. For the cube-rich ordered phase, the continuous change from a cubic mesophase into a cubic crystal is not evident

FIG. 3. Snapshots of mixture of cubes + spheres (d = 1.23 σ ) at coexisting conditions for the eutectic pressure (P∗ ≈ 8.96) from Fig. 2: (a) S = sphererich fcc solid, (b) I = isotropic liquid (center), and (c) C = cube-rich cubic ordered phase; spheres shown in red and cubes in yellow.


094102-9


J. Chem. Phys. 140, 094102 (2014)

FIG. 4. Translational (Qglobal ) and orientational (P4 ) order parameters for the cube-rich cubic phase at coexistence with the isotropic phase (below the eutectic pressure P∗ ≈ 8.96) or with the fcc sphere-rich crystal (above the eutectic pressure).

in the phase diagram of Fig. 2. However, Fig. 4 shows how the cubatic orientational order parameter60 P4 and the global translational order parameter62 Qglobal vary for the range of pressures where the ordered cube-rich phase was simulated. While the trends are somewhat noisy, one can see that while P4 changes smoothly with pressure, Qglobal exhibits two different regimes where it changes very gradually, separated by a region of rapid change between P∗ = 7.6 and P∗ = 8. Note that 6.2 < P∗ ≤ 7.6 coincides with the pressure range where mesophase behavior was detected in pure cubes.61 While the lower translational order observed for P∗ ≤ 7.6 could reflect the “disordering” effect caused by the spherical guest “impurities,” finite size effects are also playing a role in these observations. It was shown in the Ref. 62 that rotations of the cubatic order director inside the simulation box (which are more prevalent for small systems) affect translational order and the free-energy that the system can achieve. Also, our simulations showed that for 8 ≤ P∗ < 10 a system could produce a lowQglobal (∼0.3) state or a high-Qglobal (∼0.5) state depending on the starting conditions, but with the latter consistently resulting for long enough runs. This potential trapping of the ordered cubic phase in low-Qglobal ordered states is also responsible for the small differences observed in Fig. 2 between the simulated coexistence data from FENEX and the direct interfacial simulations. The latter tended to only produce the low-Qglobal ordered states, likely due to their sensitivity to interfacial fluctuations and to any slight incommensurability of the lattice unit cell to the (fixed) box cross section. One of the advantages of FENEX over conventional GDI is that the former also produces data for the ensemble freeenergy function φ along the coexistence line. In our case, the values of φ (=βμ2 here) may perhaps be less interesting than other derived properties such as, e.g., the molar entropy of each phase, which (for u = 0) can be obtained from: s / k = βP v − (y1 βμ1 + y2 βμ2 ), or using the symbols defined in Eq. (29) and thereafter s ∗ = βP v + y1 ζ − φ.

(34)

Figure 5 shows plots of this reduced entropy (values relative to an arbitrary reference value) for the isotropic and cubic phases at equilibrium. First note that, as expected, the isotropic phase has larger entropy than the cube-rich solid

FIG. 5. Relative entropy for the coexisting isotropic (I) and cube-rich ordered phase (C) in the sphere + cube system (for pressures from the pure cubes’ ODP to the eutectic pressure). Dashed and dotted lines correspond to the ideal mixing entropy –[yi ln(yi ) + (1 − yi )ln(1 − yi )] for each phase shifted to coincide with the lowest pressure point.

phase; albeit the difference is small, it increases with pressure. Note also that the isotropic-phase entropy expectedly tracks closely the changes in ideal mixing entropy solely due to composition changes. Starting from the ODP for the pure cubes, increasing pressure would be expected to have two opposing effects on each phase: To increase the mixing entropy (associated with increasing amounts on the spheres in these saturated phases, see Fig. 2) and to reduce the packing entropy (associated with a reduced free volume). Figure 5 shows that these two effects give rise to a small maximum before which entropy first slightly increases and thereafter decreases with pressure. B. Size polydisperse mixture of hard cubes

The effect of size polydispersity on hard and soft spheres has been investigated in several studies via simulations.66, 67, 70, 71, 74–76 To the best of our knowledge, only one previous study has considered size-polydisperse cubes;63 in that study, however, only the case of quenched polydispersity was considered wherein a single-phase scenario was enforced and the polydispersity was discretized and fixed (i.e., particles in the system kept their originally prescribed size) in an isobaric-isothermal ensemble. In contrast, we consider here order-disorder, two-phase states for a system with continuous polydispersity wherein particles sample different sizes constantly over a continuous range of values (SGNPT ensemble). A key advantage of such a SGNPT set up is that not all the relevant species need to be present in the simulation box at all times as the concentrations of the species are attained in an average sense (thus allowing the efficient simulation of systems with an “infinite” number of components).66, 67 Our FENEX results (Figs. 6–9) for polydisperse cubes bear some important similarities to those for spherically shaped particles.66, 67, 70, 71 In the P∗ vs υ plane (Fig. 6), υ first increases monotonically with pressure, but then it turns back toward zero while P∗ seems to diverge. As with hard spheres, this behavior is an artifact of the way how the activity spectrum is imposed which causes the average size of the particles to continuously shrinks as P∗ is increased. This effect can be partially removed by plotting βPυ 3 as shown in the


094102-10


J. Chem. Phys. 140, 094102 (2014)

FIG. 6. Order-disorder phase coexistence pressure as a function of the variance of the imposed activity distribution υ. In the inset the pressure is rescaled to βPυ 3 .

inset of Fig. 6, which shows that this renormalized pressure will approach a well-defined limiting value. A more physical osmotic pressure can be defined by rescaling the pressure by the average size of the cubes in one of the phases (chosen here to be that for the solid phase: σ C ) rather than the arbitrary σ 1 = 1 which is only meaningful near monodispersity. Likewise υ loses significance as a metric for polydispersity and hence a different, experimentally relatable quantity must be used instead such as the size dispersity “s:” s 2 = σ 2 / σ 2 − 1.

(35)

A plot of βP σ 3C vs. s (Fig. 7) shows that, like in hard spheres, a so-called “terminal” polydispersity is approached beyond which the ordered phase is no longer stable; cubes have a higher threshold (s ∼ 11.3%) compared to spheres (s ∼ 5.7%),67, 70 suggesting that the cubic order in our systems is rather robust and resilient to size imperfections. The corresponding terminal polydispersity of the isotropic phase is ∼14.3%, also larger than that found in hard spheres, but unlike the solid phase situation, this value is more dependent on the details of the imposed activity distribution.70 To simulate beyond this terminal polydispersity and possibly connect it back with the monodisperse closed-packed solid at infinite pressures, one would need to employ a distribution of chemi-

FIG. 7. Coexistence pressure of the ordered (C) and isotropic (I) phases as a function of the width of the composition distribution (s).

FIG. 8. Phase diagram in the plane of particle packing fraction and size dispersity for the ordered (C) and isotropic (I) phases. Tie lines are drawn as straight lines.

cal potential differences that goes beyond a second order term in σ ,70 a task that lies beyond the scope of this work (primarily intended to illustrate the uses of FENEX). Figure 8 shows a phase diagram expressed in terms of volume packing fractions of each phase vs. dispersity, quantities that are experimentally measurable. Interestingly and in contrast to the case of hard spheres, the volume fraction of the isotropic phase remains more or less constant over the range of coexistence conditions simulated; it appears that any increase in volume fraction that could be expected by the increase in pressure (see Fig. 7) is compensated by the increase in the fraction of the smaller particles. Indeed, fractionation is apparent between the two phases once the polydispersity increases beyond ∼5%. This is more readily seen in Fig. 9 where the composition distributions for the two phases is plotted at different points along the coexistence line: as pressure increases, the average size gets smaller and the width of the size distribution wider in the isotropic phase than in the ordered phase. Note also in Fig. 9 that the shape of the composition distributions is nearly Gaussian for all cases but the ordered phase at the highest pressure. Altogether, these results suggest that a highly polydisperse suspension of cubes can be fractionated by increasing their concentration to induce phase-separation of an ordered phase with a rather narrow size distribution (

A theoretical model for lipid mixtures, phase transitions, and phase diagrams.

Theoretical and computer simulation study of phase coexistence of nonadditive hard-disk mixtures.

Fluid-solid coexistence from two-phase simulations: binary colloidal mixtures and square well systems.

Investigation of phase transitions of lipids and lipid mixtures by sensitivity differential scanning calorimetry.

Estimation of lateral species separation from phase transitions in nonideal two-dimensional lipid mixtures.

Quantum phase transitions of atom-molecule Bose mixtures in a double-well potential.

Solid-solid phase transitions via melting in metals.

Deconvolution of High Dimensional Mixtures via Boosting, with Application to Diffusion-Weighted MRI of Human Brain.

Modeling nonequilibrium dynamics of phase transitions at the nanoscale: Application to spin-crossover.

Phase Coexistence in a Dynamic Phase Diagram.

On the determination of phase boundaries via thermodynamic integration across coexistence regions.

Binding site selection analysis of protein-DNA interactions via solid phase sequencing of oligonucleotide mixtures.

Phase Transitions in Dipalmitoylphosphatidylcholine Monolayers.

Phase transitions in plant cuticles.

Tuning bio-inspired skin-core structure of nascent fiber via interplay of polymer phase transitions.

Phase transitions in polymer monolayers: Application of the Clapeyron equation to PEO in PPO-PEO Langmuir films.

Photoinduced crystal-to-liquid phase transitions of azobenzene derivatives and their application in photolithography processes through a solid-liquid patterning.

Dynamic phase coexistence in glass-forming liquids.

Autophagy mediates phase transitions from cell death to life.

Using trading strategies to detect phase transitions in financial markets.

Critical behaviours of contact near phase transitions.

Erratum: Phase transitions via selective elemental vacancy engineering in complex oxide thin films.

Enhanced electron transport in Nb-doped TiO2 nanoparticles via pressure-induced phase transitions.

Phase transitions via selective elemental vacancy engineering in complex oxide thin films.