David A. Kofke

Direct evaluation of phase coexistence by molecular simulation via integration along the saturation line

Thermodynamic integration along a path that coincides with the saturation line is proposed as an efficient means for evaluation of phase equilibria by molecular simulation. The technique allows coexistence to be determined by just one simulation, without ever attempting or performing particle insertions. Prior knowledge of one coexistence point is required to start the procedure. Integration then advances from this state according to the Clapeyron formula—a first‐order ordinary differential equation that prescribes how the pressure must change with temperature to maintain coexistence. The method is unusual in the context of thermodynamic integration in that the path is not known at the outset of the process; results from each simulation determine the course that the integration subsequently takes. Predictor–corrector methods among standard numerical techniques are shown to be particularly well suited for this type of integration. A typical integration step along the saturation line proceeds as follows: An...

Gibbs-Duhem integration: A new method for direct evaluation of phase coexistence by molecular simulation

A method that combines the best elements of thermodynamic integration and the Gibbs ensemble technique is proposed for the direct evaluation of phase equilibria by molecular simulation. Given the conditions of coexistence at a single state point, simultaneous but independent NPT simulations of each phase are performed in succession along the saturation line. In each simulation, the pressure is adjusted to satisfy chemical potential equality according to the Gibbs-Duhem equation. Each coexistence point is determined by just one simulation, and particle insertions are never performed or attempted. Vapourliquid coexistence for the Lennard-Jones model is evaluated, and extensions are discussed.

THERMODYNAMIC AND STRUCTURAL PROPERTIES OF MODEL SYSTEMS AT SOLID-FLUID COEXISTENCE. II: MELTING AND SUBLIMATION OF THE LENNARD-JONES SYSTEM

Properties of coexisting fluid and solid phases of the Lennard-Jones model potential are reported. The melting line, representing equilibrium between the solid and liquid phases, is computed by performing a thermodynamic integration along a path that follows the coexistence line itself, according to the Gibbs-Duhem integration method recently introduced by us. Monte Carlo simulations are conducted to compute the quantities needed by the method; two system sizes, with particles numbering approximately 100 and 500, respectively, are examined. The starting point for the integration is the soft sphere system obtained in the limit of high temperature. The procedure terminates at the solid-liquid-vapour triple point, and thus the entire melting line is determined. Our results for the melting pressure are about 10% higher than those reported previously in the literature, while the triple point temperature we compute is in accord with the established value. An analysis of the approach to the high-temperature limi...

Monte Carlo simulation of multicomponent equilibria in a semigrand canonical ensemble

A general formalism and methodology are presented for the Monte Carlo simulation of equilibria in multicomponent systems, and are applied to the study of phase equilibrium in a model binary mixture, and to phase and chemical equilibrium in the ternary mixture Br2−Cl2−BrCl. Very good agreement with available experimental data is obtained. The formalism is based upon an ensemble in which the ratios of all fugacities to a reference fugacity are imposed. Simulations of mixtures performed in this ensemble fluctuate in composition while keeping constant the total number of particles. The reference fugacity is computed by integration of simulation averages of the composition for varying values of the fugacity ratios. Several features of the approach are: (1) it can be readily applied to mixtures of any number of components, even to polydisperse systems, with little additional computational effort; (2) the chemical potential of only one species at one state need be evaluated by conventional techniques, regardless...

On the acceptance probability of replica-exchange Monte Carlo trials

An analysis is presented of the average probability of accepting an exchange trial in the parallel-tempering Monte Carlo molecular simulation method. Arguments are given that this quantity should be related to the entropy difference between the phases, and results from simulations of a simple Lennard-Jones system are presented to support this argument qualitatively. Another analysis based on the energy distributions of a replica pair is presented, and an exact expression for the trial-move acceptance probability in terms of the overlap of these distributions is derived. A more detailed expression is presented using an approximation of constant heat capacity, and an asymptotic form for this result, good for large system sizes, is reported. The detailed analyses are in quantitative agreement with the simulation data. It is further shown that treatment of the energy distributions as Gaussians is an inappropriate way to analyze the acceptance probability.

THERMODYNAMIC AND STRUCTURAL PROPERTIES OF MODEL SYSTEMS AT SOLID-FLUID COEXISTENCE. I: FCC AND BCC SOFT SPHERES

Coexistence between the fluid and solid phases of systems modelled by the inverse-power potential φ(r) = e(σ/r) n is studied as a function of the potential softness s ≡ 1/n. Beginning from the fcc fluid coexistence point of hard spheres (s = 0), the Gibbs-Duhem integration molecular simulation technique is applied to trace out the solid-fluid transition pressure (reciprocal temperature) in increments of 0·01 in s, reaching a maximum s = 0·33. System sizes of the order of 500 spheres are used to model each phase, and a systematic study of finite-size effects is not attempted; thus results for s nearing its maximum are only tentative. Significant disagreement is seen with the results of early studies of inverse-power systems (for n = 12, 9, 6, and 4), while confirmation of more recent data (for n = 12 and 6) is found. Freezing into a bcc crystal is also investigated, and it is estimated that for softness s > 0·16 the bcc phase is the stable one at freezing. Coexistence between the fcc and bcc phases is not ...

Appropriate methods to combine forward and reverse free-energy perturbation averages

We consider the accuracy of several methods for combining forward and reverse free-energy perturbation averages for two systems (labeled 0 and 1). The practice of direct averaging of these measurements is argued as not reliable. Instead, methods are considered of the form β(A1−A0)=−ln[〈w(u)exp(−βu/2)〉0/〈w(u)exp(+βu/2)〉1], where A is the free energy, β=1/kT is the reciprocal temperature, u=U1−U0 is the difference in configurational energy, w(u) is a weighting function, and the angle brackets indicate an ensemble average performed on the system indicated by the subscript. Choices are considered in which w(u)=1 and 1/cosh[(u−C)/2]; the latter being Bennett’s method where C is a parameter that can be selected arbitrarily, and may be used to optimize the precision of the calculation. We examine the methods in several applications: calculation of the pressure of a square-well fluid by perturbing the volume, the chemical potential of a high-density Lennard-Jones system, and the chemical potential of a model for ...

Surface tension and vapor–liquid phase coexistence of the square-well fluid

Vapor–liquid interfacial tension of square-well ~SW! fluids is calculated using three different methods viz., molecular dynamics ~MD! with collision-based virial evaluation, Monte Carlo with virial computed by volume perturbation, and Binder’s density-distribution method in conjunction with grand-canonical transition-matrix Monte Carlo ~GC-TMMC!. Three values of the SW attractive well range parameter were studied: l51.5, 1.75, and 2.0, respectively. The results from MD and GC-TMMC methods are in very good mutual agreement, while the volume-perturbation method yields data of unacceptable quality. The results are compared with predictions from the statistical associating fluid theory ~SAFT!, and SAFT is shown to give a good estimate for the systems studied. Liquid and vapor coexistence densities and saturation pressure are determined from analysis of GC-TMMC data and the results are found to agree very well with the established literature data.

Coexistence diagrams of mixtures by molecular simulation

Abstract We present extensions of the Gibbs—Duhem integration technique that permit its application to mixtures. The Gibbs—Duhem integration method provides a means to evaluate the equilibrium properties of coexisting phases efficiently and reliably by molecular simulation. In particular, the method prescribes how a sequence of simulations can be performed at state conditions that lie on the coexistence surface, thus it is particularly well suited for the mapping of complete phase diagrams. When applied to pure substances, the Clapeyron equation guides the choice of state conditions as the series proceeds. We present two extensions, based respectively on “semigrand” and “osmotic” forms of a generalized Clapeyron equation that we also present. Both methods require sampling of compositions during the simulation, but they differ in how this is accomplished. Simulations based on the semigrand method have steps in which molecules attempt to change their species identity, keeping the total number of molecules fixed; the osmotic simulations have moves in which insertions or deletions are attempted for one of the species, while the number of molecules of the other species remains always fixed. Each approach is seen to have certain advantages that depend on the nature of the mixture being studied. Both methods are demonstrated with application to three prototype binaries, and extension to multicomponent mixtures is discussed.

Selection of temperature intervals for parallel-tempering simulations

A heuristic is developed specifying that temperatures in replica-exchange simulations should be spaced such that about 20% of the phase-swap attempts are accepted. The result is found to be independent of the heat capacity, suggesting that it may be applied generally despite being based on an assumption of (piecewise-) constant heat capacity.