Coexistence diagrams of mixtures by molecular simulation

Pergamon

Chemical Enaineerins Science. Vol. 4% No. 16. w. 2633-2645.

1994

0009-2509(94)EOO78-5

COEXISTENCE

Department

DIAGRAMS OF MIXTURES SIMULATION

BY MOLECULAR

MANOJ MEHTA and DAVID A. KOFKE of Chemical Engineering, State University of New York at Buffalo, Buffalo, NY 14260-4200, U.S.A. (First receioed

10 January

1994; accepted

18 March

1994)

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.

1. lNTRODUCTION Molecular simulation is a computational technique for determining the thermodynamic and transport properties of systems defined by model intermolecular potentials (Allen and Tildesley, 1987). Properly executed, simulation can provide most properties of interest to an arbitrary degree of precision. The role it plays in bridging molecular theories with experiment is widely appreciated, but there is now substantial effort being put forth to develop it for use in a new role. Simulation may ultimately see its widest application as a predictive tool, much as a theory is, which can supplement limited experimental information. The application of simulation in this role means different things to different groups: biochemists and protein engineers for example apply simulation to guide them in designing new proteins, while the properties of greatest interest to chemical engineers include those that pertain to phase equilibria. In any case, success in the development of simulation as a predictive tool relies heavily on the establishment of accurate models of the intermolecular interactions. To optimize these models for the study of phase coexistence, it is esscntial that we have efficient and reliable techniques for the evaluation of coexistence diagrams. Naive approaches to the simulation of phase equilibria-to mimic experiment, for example, by observing directly the equilibration of coexisting phases in the same simulation volume-fail because the interfacial properties between the (necessarily small) phases overwhelm their individual bulk properties (Rowlinson and Widom, 1982; Gubbins, 1989). Determination of equilibrium between phases can be per-

formed using a (tedious) search for the conditions that satisfy the formal requirements of phase coexistence: equality of temperature, pressure, and all species fugacities between the phases (Adams, 1976, 1979; et al., 1986; Fincham et al., 1986; Panagiotopoulos Kotke and Glandt, 1988). Evaluation of these quantities in principle requires a separate simulation for each phase. Thus the calculation of a phase diagram by this means is an enormous undertaking, as determination of each point on the coexistence surface entails many simulations at “uninteresting” state conditions; the burden is compounded by the well-known difficulty of calculating fugacities by molecular simulation. In the late 1980s the Gibbs ensemble method was introduced (Panagiotopoulos, 1987) and refined (Panagiotopoulos et al., 1988; Smit et al., 1989) by Panagiotopoulos and coworkers, and this picture changed entirely. Panagiotopoulos very cleverly pointed out that both phases may indeed be simulated at once if they are removed from physical contact while remaining in thermodynamic and material contact (i.e. the total mass and volume are fixed). This feat is accomplished by performing the simulations of the coexisting phases simultaneously and, for the most part, independently. Fluctuations in the volumes of the two phases are coupled to ensure pressure equilibration, and particles are exchanged between them-in grand-canonical fashion-to ensure fugacity equilibration. The Gibbs ensemble has since been applied by many research groups to study coexistence (mostly vapor-liquid coexistence) in a wide variety of model systems, and it has enjoyed a spectacular degree of

2633

2634

MANOJ MEHTA and DAVID A. KOFKE

success (Panagiotopoulos, 1992). Nevertheless, the Gibbs ensemble as originally formulated suffers from certain limitations introduced by the particle exchange moves which enforce fugacity equilibration. Simulation of solid-fluid equilibria, for example, has not been accomplished with the Gibbs ensemble and, despite remarkable progress, limitations have been encountered in the study of coexistence in polymer systems (Laso et al., 1992; Mooij et al., 1992). Promising alternatives to the Gibbs ensemble certainly exist: thermodynamic-scaling Monte Carlo, a recent generalization (J. Valleau, personal communication) of density-scaling Monte Carlo (Valleau, 1991), does not rely on the simultaneous simulation of coexisting phases, and it is remarkable in the precision of the results it has provided for vapor-liquid coexistence, particularly near the critical point; block density distribution function analysis (Rovere et al., 1988, 1990, 1993) takes finite-size scaling concepts off the lattice to yield coexistence data, also using a single simulation volume and also with special attention to the critical point. These methods are still under development, and they have yet to be applied to coexistence other than vapor-liquid, or even to mixtures, although this may be only a matter of time (Marx et al., 1993). To overcome the problems of simulating complicated coexistence phenomena with the Gibbs ensemble method, the particle-exchange step must be extracted from the algorithm. We must recognize that the great contribution of Panagiotopoulos is the idea of simulating two phases that are coupled thermodynamically and materially but are otherwise independent. Particle exchange provides a particularly convenient means for introducing both thermodynamic and material coupling, but there are other ways to achieve the same while retaining the elegance and simplicity of the Gibbs ensemble technique. We may for example couple the phases materially by performing two types of volume change steps (Kolke, 1993a): (1) the usual move in which volume is exchanged between the phases in a manner that conserves the total volume, and (2) virtual-particle exchanges, in which the volumes are adjusted simultaneously to mimic the changes in density that accompany transfer of particles. Given this means for coupling the phases materially, we are free to select any other convenient method to introduce thermodynamic coupling. In recent work, we presented “Gibbs-Duhem integration” as a means to this end (Kofke, 1993a, b). The method entails simultaneous simulation of each phase in the NpT ensemble (i.e. the particle number, pressure, and temperature of each phase are imposed and do not vary), so thermal and mechanical equilibrium are enforced directly. The mechanism for equating the fugacities is the Clapeyron equation. This formula specifies how the pressure must change with temperature for two coexisting phases to remain in coexistence. Starting at a state condition for which the two phases are known to be in equilibrium, the Gibbs-Duhem integration method can be used to

trace out the phase diagram directly and efRciently. Each simulation determines a point on the coexistence curve, and particle exchanges are never required. The purpose of this report is to demonstrate the application of the Gibbs-Duhem integration method to mixtures, The tamiliar form of the Clapeyron equation is derived in the p-Tplane, but it is relatively easy to extend this formula to describe coexistence lines in other planes. In the next section, we review the thermodynamic formalism for mixtures, and from it we develop Clapeyron-type formulas for the coexistence lines in composition-dependent planes. Two formulas are presented: one prescribes simulations in a “semigrand” representation, while the other requires “osmotic” simulations. Each approach will be seen to have particular advantages depending on the type of system being simulated. In Section 3, we describe how the two formulas are implemented in simulation, and in Section 4 we demonstrate application with three well-studied model mixtures. We conclude with a summary and suggestions for further study. 2. FORMA-M

A perhaps unfamiliar aspect of the development that unfolds below is the treatment of the species chemical potentials (or, more precisely, the fugacities) as independent variables. This is of course a perfectly valid means of defining the thermodynamic state of the system; in fact, researchers who deal with adsorption or membrane equilibria find that such an approach arises naturally. Such a point of view is easily and routinely accommodated in Monte Carlo simulation (Allen and Tildesley, 1987), and we adopt it because it is the appropriate way to extend our previous work to mixtures_ For ease of presentation, we will specifically discuss application to a mixture of two components; extension to mixtures of more species will be discussed briefly in the concluding section. We begin by considering the Gibbs-Duhem equation (Denbigh, 1971) for a binary mixture, which may be written in the form xrdlnfi

= h.d/I+Zdlnp

-xadlnfz.

(2.1)

Here p is the pressure, and 2 is the compressibility factor, pv/RT, where v the molar volume; xI is the mole fraction of species i, with fugacity fi; h, is the residual molar enthalpy, defined as the enthalpy above an ideal gas at the same temperature; finally, B = l/R T, where R is the gas constant and T is the absolute temperature. The fugacity fraction (GrifRths and Wheeler, 1970; Koflce and Glandt, 1988; Kofke, 1991) of species 2

fz

52=fl +f1

(2.2)

is a convenient quantity that varies from zero to unity as the mixture composition goes from pure species 1 to pure species 2. Its introduction permits the Gibbs-Duhem equation to be written in the alternate

Coexistence

diagrams

of mixtures

but equivalent form x2

dln(f,+f,)=h,dB+Zdlnp-

tat1

-

52 -

d52

P2)

(2.3) we note that this formula is symmetric with respect to interchange of the species labels. Both r2 and the sum (fl +f2) are field-type intensive variables (Griffiths and Wheeler, 1970), and must therefore have the same values in two coexisting phases. For reasons that will become clear below, we shall designate eq. (2.1) as the “osmotic” form of the Gibbs-Duhem equation, while eq. (2.3) will be referred to as the “semigrand” form. Clapeyron-type formulas are developed by considering changes in the intensive “field” variables that maintain equilibrium between coexisting phases (Denbigh, 1971). For example, if two phases a and y are in mutual equilibrium, then change in temperature, pressure, and species-2 fugacity must occur in such a manner that the change in the species-l fugacity is the same in both phases. Thus, according to eq. (2.1) ($-$)dp+($-5)

by molecular

2635

simulation

ing the behavior of these formulas as one of the components becomes infinitely dilute. As this limit is approached, the fugacity of the dilute component obeys Henry’s law:J, = nixi, while the abundant species follows Raoult’s law: fi =$ xJ (Hi is the Henry’s constant, and fro is the fugacity of pure component j at the temperature of the mixture and at its saturation pressure). We consider first the limiting behaviour as x2 + 0. Equation (2.5) is easy to describe in this limit since 42 + Hz/p, while eq. (2.7) requires us to note that fugacity fraction l2 vanishes as x2 H2/fF _ Then

=PWH: -

(--> x2

#.a

=f?WHi Z’-

(2.4)

the subscript u indicates a change along the saturation line, and the fugacity coefficient @2 =f2/(px2); similar formulas for changes at constant pressure or constant species fugacity are easy to derive from eq. (2.4). A equivalent development may be advanced from the semigrand form of the Gibbs-Duhem equation. In this case, we consider variations in temperature, pressure, and fugacity fraction that keep changes in the sum (ft +p2) equal between phases. Thus, (h:-&‘)d/3+(Z’+Z@)dlnp

and, once more, variations at constant temperature obey

x2

I

p..¶= t2(1 -x;2;(5:

* - Z@).

(2.7)

Before proceeding to describe how eqs (2.5) and (2.7) are used to evaluate phase diagrams by molecular simulation, we will conclude this section by consider-

(2.8b)

In the latter formula we have used the fact thatfP has the same value in the liquid and vapor. The limiting formulas as component 1 becomes dilute can likewise be expressed in terms of the Henry’s constant of species 1 in otherwise pure species 2, thus (2.9a)

(x2 --, 1) =

_f,"WH: - l/HI) Z'- iv

(2.9b)

where 4: is the fugacity coefficient of pure species 2, and has the same value in both phases. 3. SIMULATION

(->

Z’

H: -H:

where the superscripts indicate values in phase 1 or g (liquid or gas), as appropriate. In particular, changes that maintain coexistence at constant temperature are governed by a Clapeyron differential equation, which may take the following form

p

- l/HZ)

#&Hi Z’ - H: Z") d In A = 0

(2.8a)

(x2 - 0) dlnp

dlnp

-(+$)

a In

l/H:)

Z’ - ZB

METHOD

Each Clapeyron formula implies a specific simulation methodology for evaluating coexistence diagrams. Equ&on (2.5) prescribes simulation in an osmotic ensemble, for which the independent variables include temperature, pressure, and species-2 fugacity; the remaining variable-which must be extensive-is the mole number of species 1, n, (of course, since we cannot simulate systems of macroscopic size, this takes the form of the number of molecules of species 1, N,). Equation (2.7) specifies simulations in a semigrand ensemble (Brian0 and Glandt, 1984), which has as independent variables the temperature, pressure, fugacity fraction, and total number of molecules N. The goal of the procedure is to determine the coexistence properties of a model fluid over some range of composition. An “initial condition” is required to begin the process, that is, the coexistence properties must be known at one composition--one of the pure fluids may in fact serve in this role. Given this information, either of the Clapeyron differential equations may be integrated numerically to trace out the corresponding coexistence line; predictor-corrector methods are an appropriate integration scheme for this purpose (Kofke, 1993a). Each step of the integration proceeds as follows. The independent vari-

MANOJ MEHTA and DAVID A. KOFKE

2636

able--f or < 2, depending on Clapeyron formula used to govern the process-is incremented, and a “predictor” pressure is computed from the initial condition and/or prior simulation data. Simultaneous Monte Carlo simulations of the two coexisting phases commence at the specified values of T, p andf, or t2. As the simulation proceeds, thermodynamic averages are recorded and the “integrand” of the appropriate Clapeyron formula is estimated. This estimate is used to “correct” the simulation pressure while the simulation is in progress; appropriate predictor/corrector formulae are presented in Table 1. In this table each formula describes how the pressure (designated here as y) in a given simulation (subscripted i + 1) is determined from y and F [the right-hand side of eq. (2.5) or (2.7)] of previous simulations (subscripted i, i - 1, etc.), corresponding to _fz or r 2 that differ by h. For example, if h = 0.05 and the current simulation (i + 1) is conducted at lz = 0.5, then yi and Fi are In(p) and F(&, p) for
identities, just as they sample positions within the simulation box. The mole fractions in each phase are then computed as averages, and they enter into the integration as specified by eq. (2.7). The procedure for performing simulations in the semigrand ensemble is well documented (Kofke and Glandt, 1987, 1988; Kofke, 1991), and it may be implemented here without significant complication; details are reviewed in the appendix. In the osmotic implementation, simulation proceeds according to eq. (2.5). Again, two systems are simulated simultaneously, but now at fixed T, p, f2, and N? . While Nf is imposed and does not vary, insertion and deletion of species-2 molecules must be regularly attempted in each simulation volume, and as a result N,” does fluctuate during a simulation. The average value of N; will be determined by, among other things, the imposed value of Nf; because they are both extensive variables, the average N; will increase proportionately with NT at fixed T, p, and f2 _In particular, if it is desired to map the coexistence diagram from pure species 1 to pure species 2, one must take care to decrease N; as the series progresses (i.e. from one simulation to the next) so that the system size does not become unwieldy. In both the osmotic and semigrand approaches, the two simulated systems are not necessarily coupled materially: volume changes and partide insertions (for osmotic simulations) occur independently in each phase. This feature contrasts with the Gibbs ensemble method, which couples such changes so that the total volume and particle number of the composite system is fixed. In some situations, it is desirable that the phases be so coupled. As mentioned in the Introduction and described further in the appendix, material coupling can be introduced artifically, without performing particle exchanges, by properly coupling the volume changes. The effect of such coupling is to force the densities of the two systems to remain on opposite sides of an arbitrarily chosen value. This constraint prevents either phase from unilaterally condensing or vaporizing, and thereby defeating the method. Henry’s constants are required to connect the integration to the pure-component limits. For the semigrand simulations, only the ratio of the solute Henry’s constant to the solvent fugacity is required (we designate the solvent component which is not at infinite

Table 1. Predictor-corrector Simulation in series

First

formulae Corrector formula

Predictor formula

h Yi+l=Y*+-(Fi+l+F*)

Yr+ 1 = Yr + hK

2

Y;+I

=

Y,-I

+

h Yl+i=~i-l+-(Fi+l+4F;+F,_~)

ZhFi

3

Y,+~ = yi-,

All subsequent

Y,+I=~,+$(SSF~--~F~-~

-+- 2hFi

+ 37F,_,

yi+i

- 9F,-,I

= yi

h

++9Fi+l

yi+, = yi + $(9F,+i

+ 19F, -

5Fi_I

+ F,_2)

+ 19Fi - 5F,_,

+ Fi-Z)

Coexistence diagrams of mixtures by molecular simulation dilution). This ratio can be determined in a puresolvent simulation by performing trial identity changes: one of the solvent molecules is converted into a solute molecule and the “exchange energy” Au is noted, the particle is then converted back before the simulation proceeds. The ensemble average of the Boltxmann factor of the exchange energy gives the desired ratio Hi f~ = -t

.

(3.1)

The Henry’s constants needed for the osmotic series must be evaluated with trial insertions. A molecule of the dilute species is inserted into an otherwise pure solvent, and the energy Au of this test particle is noted. Then Henry’s constant is given by the average Hi=kT

N+l v

exp ( - Au/k T)

-1

(3.2) > where V is the system volume and N is the number of solvent molecules (which fluctuates according to the imposed value of fz if the solvent is species 2, but is fixed if the solvent is species 1). Each simulation approach, taken alone, provides an incomplete set of thermodynamic data. The semigrand method for example does not give directly the fugacity of either species. To determine them, we must return to the semigrand Gibbs-Duhem equation, eq. (2.3). Knowing the sum (f, + fi) at the initial coexistence point, we may evaluate it at any other state along the coexistence line by integration of eq. (2.3). For this purpose it is best to reduce the integral to a onedimensional form (

ln(~+f2)=ln(~~+f~)+~~~~(~),.. x9-52 -

52u

-52)

1

(3.3)

dr2

where the 0 superscript indicates a value at the initial coexistence point. The derivative in the integrand is given by eq. (2.7), while 2 and xz may be taken as the simulation averages from either phase a. Knowledge of both the fugacity sum and the fugacity fraction allows either fugacity to be computed easily. A similar analysis-using eqs (2.1) and (2.5tmay be applied to evaluate the species-l fugacity at any point in an osmotic series. Each method for evaluating coexistence diagrams of mixtures has certain advantages depending on the system under study. The semigrand method never requires particle insertions, but it does demand that particles be able to change species identities. The osmotic approach does not require such identity changes, but is must be possible to insert particles of one of the species. The semigrand method thus is best suited for mixtures of species that are not too dissimilar, although each component may be complicated enough to preclude insertion; many simple systems fit this description, along with some types of polymer melts. With the osmotic method, the species may be ES

49:16-G

2637

quite dissimilar, but as Tong as one of them is insertable, the method can be applied; a good example is a polymer in solution. A more subtle measure of the relative merits of the two methods relates to their stability properties: do errors grow or diminish as the series proceeds? There are several ways to address this question (Hamming, 1973; Finlayson, 1980), but we will examine only one here. If we assume that the only error in the results arises from an incorrect initial pressure, we can examine the growth of this error as the integration series proceeds. An incorrect initial pressure would imply that the pure-phase fugacityf: is not the same in both phases (choosingfs = 0 for the initial case)

fP’ = (1 + &)fpJ

(3.4)

where the error E vanishes as the initial pressure approaches its correct saturation value. The integration procedure is constructed to keep the ratio of fugacities (f//fiu) constant; if the initial pressure is chosen correctly, this constant is unity, and all subsequent simulations will follow the true coexistence line. Otherwise, at a given point in the integration series, we will be simulating two phases at the same T, p. andfi (or &), but with fugacities & related by eq. (3.4). We would like to know to what extent this causes the pressure to deviate from the correct coexistence value. Assuming a small difference, we linearize the dependence offi on p, and find for the osmotic case

dp -= P

C(+I)’ - Wxt)‘lo

-dp

0 P

0

(-m,)‘-

WXI)”

.

(3.5)

Here dp is the deviation of saturation pressure from its correct value, and the other symbols are as defined above; the 0 subscript indicates a value at the initial condition. Assuming that the integration starts from pure 1, as it progresses the mole fraction of species 1 decreases and the error in pressure--on a percentage basis--becomes smaller. In the limit xi, XT -, 0 (i.e. pure 2), dpjp + 0 as well. Thus, even though the series begins with an incorrect pressure at one end, the correct saturation pressure is approached at the other. Again, this analysis assumes that the only error in the procedure is that due to an incorrect initial pressure; unaccounted sources of error include the statistical errors of the simulation averages, and the error inherent in the predictor/corrector integration procedure. However, as each simulation provides, in a sense, the “initial condition” for subsequent simulations, the analysis is useful as an indicator of the inherent stabity of the method. A similar analysis for the semigrand method shows (3.6) Here, unlike the osmotic case, the error does not go to zero in the limit of pure 2. The formula can be used to guide the integration procedure: one should when possible integrate in a direction in which (Z’ - ZB) is increasing. Of course, this information is not always

2638

MANOJ MEHTA and DAVID A. KOFKE

Table 2. Model parameters for the mixtures examined in this study

0.2

0.4

0.6

0.8

Liquid-phasemole fraction.x2 Fie. 1. Prouaeation of error in saturation messure. comp&d by ap’pl$ng Gibbs-Duhem integration ‘methodb to the van der Waals model fluid. Model parameters and temperature are selected roughly to mimic the isotherm for mixture II defined in Section 4 of the text. Series begins with a 5% error in the pressure for pure species 1, and proceeds to the right.

known a priori, but it can be monitored during the series if error in the initial pressure is a concern. To demonstrate the above results in error propagation, we applied the Gibbs-Duhem integration technique in the van der Waals model fluid, for which all properties may be determined analytically. We introduced a 5% error in the initial pressure, and proceeded with the integration for each method. The integration step was chosen small enough that the error it introduced could be considered negligible. Figure 1 shows that as the integration proceeds, the percentage error in the saturation pressure approaches zero using the osmotic method, and varies slowly (in this case, decreasing) using the semigrand method. This analysis confirms the stability of both methods with respect to a small error in the initial pressure. Incidentally, a similar analysis for a pure-fluid Gibbs-Duhem integration (pressure vs. temperature) also results in eq. (3.6). Thus, as the critical point is approached, AZ = (Z’ - ZB) approaches zero, and the error in the pressure diverges. The analysis indicates that, in general, it is best to choose the initial condition at high temperature and to integrate toward smaller values of T. 4. APPLICATIONS

We applied both Gibbs-Duhem integration methods to evaluate a single isotherm for each of three Lennard-Jones (LJ) binary mixtures. The LJ parameters for each mixture are presented in Table 2. The components of mixture I are identical, but the cross-interaction energy parameter is weaker than either pure-component value; the components of mixture II have the same energy parameters, but are of different size; mixture III has both the energy and size parameters different, and one of the components is supercritical. The coexistence properties of these mixtures have been studied previously using the Gibbs

Mixture

=11

Cl2

~22

El1

81,

I II

1.0 1.0 1.0

1.0 0.769 0.768

1.0 1.0 1.0

0.75

III

1.0 0.885 0.884

822 1.0

1.0

1.0

0.713

0.597

ensemble (Panagiotopoulos et al., 1988), so we have data with which to verify our calculations. We present our results as projections in three planes: px, P_<~, and PJ~, respectively; the first of these is the usual form in which mixture coexistence data are presented, while the latter two are of interest because they show the path of integration taken by the semigrand and osmotic methods, respectively. Details of the simulation may be found in the appendix, Mixture I is symmetric with respect to interchange of the species labels, and this simplification is reflected in the pressure-composition diagram. Nevertheless, to test the accuracy of the technique we used it to compute the entire phase diagram. The results are shown in Fig. 2. The pressure-composition diagram exhibits the required symmetry, and is in good agreement with the Gibbs ensemble data. The p<, and the p-f’ plots begin to demonstrate a point we wish to make concerning the relative advantages of the semigrand and osmotic approaches. Integration of the Clapeyron equations can be done most easily when the integrand is smooth and slow-varying. In such an instance, a relatively large integration step may be used, or equivalently, the error is smaller for a given step size. Moreover, a large change in the magnitude of the integrand implies a marked change in the slope of the coexistence line. To negotiate this transition, the integration step size should be modified accordingly. An adaptive integration algorithm can be introduced for this purpose, but it would sometimes be preferable to avoid the problem altogether. Figure 2 shows that for mixture I both integration pathways have significant curvature, but that the osmotic pathway drops precipitously after reaching a maximum. The integration step is not decreased to compensate for this change in the path, and as consequence the calculations lose accuracy. The semigrand pathway does not suffer from a similar feature, so reliable results can be obtained over the entire range of composition. The species in mixture II differ only in their size parameters. The pressure
Coexistence diagrams of mixtures by molecular simulation

2639

0.08

0.06

0:o

0:2

014

0:6

0.0

0.2 0.4 0.6 0.8 Mole fraction of component2, x2

1.0

0.0

0.4 0.2 0.6 0.8 Fugacity fraction of component2, &

1.0

0.00

0.02 0.04 0.06 Fugachy of component2, 5fzanf

l:o

0:8

Mole fraction of component2,x2

0.06 0.0

0.8 0.2 0.4 0.6 Fugacity fraction of component2,cr

1.0

I....

0.10: (c)

0.08

0.08-

0.06

0.00

I”.’

0.01 0.02 0.04 0.03 Fugacityof component2, Sf+,,’

Fig. 2. Vapor-liquid

coexistence diagrams for mixture I at (~,T/E,~) = 1.15. (00 0) semigrand series (lines connect the points as a guide to the eye); (+ + +) osmotic series; (x x x ) Gibbs-nsemble data of Panagiotopoulos et al. (1988). (a) pressure-composition diagram; (b) semigrand integration pathway; (c) osmotic integration pathway [( x ) are Gibbs-ensemble data for the liquid, ( + ) are for the vapor]. reduced temperature

1993a). Unlike the semigrand method, an osmotic series cannot connect to pure species 2. In this limit, NT must be zero, and there remains no extensive variable to specify the system size. Mixture III posed several difficulties. Figure 4 indicates that the advantages of the two pathways are reversed from that seen in mixture I: the osmotic path is straight, while the semigrand path rises steeply after a period of slow increase. Moreover, one of the species is supercritical, so the coexistence line terminates at

0.08

Fig. 3. Vapor-liquid coexistence diagrams for mixture II at reduced temperature (k, T/s1 I) = 1.15. Symbols are as in Fig. 2, with the following additions: (0 0 0) “reverse” semigrand series starting from pure component 2; (A A A) osmotic series conducted with double the number of particles (totalllng 512) and half as large an integration step as used in the (osmotic) series represented by the filled diamonds.

a critical point. As the transition weakens upon approach of the critical, the difference in the properties of the two phases diminishes. This has three consequences. First, the phases become more likely to unilaterally vaporize or condense, so the artificial material coupling becomes essential. Second, the integrand of the Clapeyron equation-which is the ratio of differences between the phases-becomes increasingly susceptible to statistical errors in the simulation averages. Third, as implied by eqs (3.5) and (3.6). the methods begin to lose stability, and errors in the pressure tend to grow. Figure 5 demonstrates this trend. This figure presents the denominators of the

MANOJ MEBTA and DAVID A. KOFKE

2640

0.0

0.2

0.0

0.4

0.6

1.0

0.8

Mole fraction of component2. x2

010

0.2

0.4

0.6

,+b.” +,t

0.15-

1.0

0.8

ax0

: (c) 0.20 -

Fig. 5. Inherent stability indicator of the semigrand and osmotic integration series for mixture III. Plotted are the denominators of the right-hand side of eq. (3.5) (right ordinate) and eq. (3.6) (left ordinate) as computed in the simulations. Both integration series proceed from left to right in the figure.

:

o

5. CONCLUDING

0.10 -

0.05 -

OiKl

0% o.b5 o.io Fugacityof compoaent2, Bf20,,

end of the simulation range. Corresponding to this outcome is a diverging of the measured saturation properties as given by the two techniques. The difference is consistent and reproducible, and does not seem to result from statistical errors in the simulation data. It survives, for example, changes in system size, integration step size, and simulation duration (data not included in the figure). According to Fig. 5, this discrepancy should be resolved in favor of the osmotic method. We note that the scatter in the Gibbs-ensemble data in this region is comparable to the difference between the osmotic and semigrand results.

o.io

Fig. 4. Vapor-liquid coexistence diagrams for mixture III at reduced temperature (k, T/e, ,) = 0.928. (0 0 0) semigrand series; filled symbols describe several osmotic series, which differ in system and integration step size: diamonds (216 total particles, 0.008 step in fugacity); triangles (512 particles, 0.016 step); circles (512 particles, 0.008 step). Gibbs-ensemble data are represented with crossa as in Figs 2 and 3.

right-hand side of each of these equations, plotted as a function of the saturation pressure; the critical point is found toward the right-hand side of the figure. The approach of these quantities to zero here is an indication of instability of the integration procedure. It is clear from the figure that the semigrand method suffers most from this effect. The osmotic stability indicator has just begun its descent toward zero at the

REMARKS

We have presented two approaches for implementing the Gibbs-Duhem integration method to evaluate the coexistence properties of mixtures. Each method has its particular merits, and the choice of which to use must be based on the nature of the mixture being studied. The main consideration is the means by which compositions can be most easily sampled: if identity changes can be routinely accepted, then the semigrand approach is best, while if insertion/deletion trials are viable for at least one of the components, then the osmotic approach can be considered. If approach to a critical point is involved, the inherent stability of each method should be examined. Finally, the shape of the projection of the coexistence surface in the integration plane can affect the ease with which the integration is implemented. Although an intelligent integrator with an adaptive stepsize algorithm can easily navigate any coexistence line, a straighter line generally entails fewer simulations. Unfortunately, the shape of this line cannot in general be assessed prior to performing the integration series, so this consideration cannot be applied to select an integration method; experience with the technique may eventually provide guidelines in this area. Both the semigrand and osmotic integration methods can be extended to mixtures that contain

Coexistence diagrams of mixtures by molecular simulation more than two species. The starting point is the generalization of the Gibbs-Duhem formulas presented in eqs (2.1) and (2.3); for a mixture of c components, these become, respectively, xldlnfi=h,dB+ZdInp-

tx.dlnf.

(5.1)

n=2

dln

!l=1 1 k_&

[

=h,dB+Zdlnp-

s=i” i

pd<,,.

(5.2)

The form of eq. (5.2) maintains a symmetry among all the components. We note that the sum on the righthand side of this formula is written over all species; however, since the fugacity fractions must sum to unity, the de. are not all independent [e.g. if c = 2, eq. (2.3) results upon realization that dtl = - dr2]. The main complicating feature of multicomponent mixtures is the proliferation of potential integration pathways that accompanies the addition of more species. Once one decides how the fugacities or fugacity fractions are to vary’ as the integration proceeds, one can go on to derive an appropriate Clapcyron formula. The methods presented in this report may then be applied directly. The choice of whether to apply Gibbs-ensemble Monte Carlo or Gibbs-Duhem integration will be influenced by several factors. Both methods--in their basic form-are equally easy to implement from a programming standpoint, If the system is such that Gibbs-ensemble particle transfers can be accepted with sufficient frequency (even as few as 1% of attempted transfers is enough), then the computer time required for the two methods is comparable (per coexistence point). Gibbs ensemble is certainly the method of choice in this instance, because it has the very important advantage of being able to provide coexistence data at arbitrary points on the phase diagram. As the systems or coexistence phenomena become more complex (involving polymers or solids, for example), the issue becomes less clearcut. If it can) be applied at all, Gibbs-ensemble Monte Carlo requires longer sampling or special insertion techniques in these cases. It may still be the preferred approach if arbitrary coexistence points are desired, or if a “nearby” coexistence datum is not easily established for initiating a Gibbs-Duhem series. However, if one is interested more in determining the complete phase behavior of a model substance, or if particle exchange is simply not viable, the advantage leans toward Gibbs-Duhem integration. The choice between the two for specific cases will certainly become clearer as we continue to learn more about both methods. Acknowledgements-Financial support for this work has been provided by the National Science Foundation, under grant CTS-9212682, and under the Presidential Young Investigator program. REFERENCES Adams, D. J., 1976, Calculating the low-temperature line by Monte Carlo. Molec. Phys. 32, 647457.

vapaur

2641

Adams, D. J., 1979, Calculating the low-temperature vapour line by Monte Carlo. Mokc. Phys. 37, 211-221. Allen, M. P. and Tildesley, D. J., 1987, Computer Simulation of Liquids. Clarendon Press, Oxford. Briano. J. G. and Glandt, E. D., 1984, Statistical thermodynamics of polydisperse fluids. J. them. Phys. SO. 33363343. Denbigh, K., 1971, Principles of Chemical Equilibrium. Cambridge University, Cambridge. Fincham, D.. Quirke, N. and Tilde&y, D. J., 1986, Computer simulation of molecular liquid mixtures. I. A diatomic Lennard-Jones model mixture for COdC,H,. J. them. Phys. 84,4535-4546. Finlayson, B. A., 1980, Nonlinear Analysis in Chemical Engineering. McGraw-Hill, New York. Griffiths, R. B. and Wheeler, 3. C., 1970, Critical points in multicomponent systems. Phys. Reu. R 2. 1047-1064. Gubbins, K. E., 1989, The role of computer simulation in studying phase equilibria. Molec. Sim. 2, 223-252. Hamming, R. W., 1973, Numerical Methodsfor Scientists and Engineers. Dover, New York. Kofke, D. A., 1991, Solid-fluid coexistence in binary hard sphere mixtures by semigrand Monte Carlo simulation. Molec. Sim. 7, 285-302. Kofke, D. A., 1993a, Direct evaluation of phase coexistence by molecular simulation via integration along the saturation line. J. them. Phys. 98,414%4162. Kotlce, D. A., 1993b, Gibbs-Duhem integration: a new method for direct evaluation of phase coexistence by molecular simulation. Molec. Phys. 78, 1331-1336. Koike, D. A. and Gland& E. D., 1987, Nearly monodisperse fluids. I. Monte Carlo simulations of Lennard-Jones particles in a semigrand ensemble. J. them. Phys. 87, 4881-wO. Kofke, D. A., Glandt, E. D., 1988, Monte Carlo simulation of multicomponent equilibria in a semigrand canonical ensemble. M&c. Phys. 64, 1105-1031. Laso, M.. de Pablo, J. J. and Suter, U. W., 1992, Simulation of phase equilibria for chain molecules. J. &em. Phys. 97, 2817-2819. Marx, D., Nielaba, P. and Binder, K., 1993, Path-integral Monte Carlo study of a model adsorbate with internal quantum states. Phys. Rev. EJ 47, 7788. Mooij, G. C. A. M., Frenkel, D. and Smit, B., 1992, Direct simulation of phase equilibria of chain molecules. J. Phys.: Condens. Matter. 4. L255-L259. Panagiotopoulos, A. Z., 1987, Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble. M&c. Phys. 61, 813-826. Panagiotopoulos, A. Z., 1992, Direct determination of fluid phase equilibria by simulation in the Gibbs ensemble: a review. Molec. Sim. 9, l-23. Panagiotopoulos, A. Z., Quirke, N., Stapleton, M. and Tildesley, D. J.. 1988, Phase equilibria by simulation in the Gibbs ensemble: alternative derivation, generalization and appiieation to mixture and membrane equilibria. Mo[ec. Phys. 63, 527-545. Panagiotopoulos, A. Z., Suter, U. W. and Reid, R. C, 1986, Phase diagrams of non-ideal fluid mixtures from Monte Carlo simulation. Ind. Engng Chem. Fundam. 25, 525535. Revere, M., Heermann, D. W. and Binder, K., 1990, The gasliquid transition of the two-dimensional LennardJones fluid. J. Phys.: Condens. Matter. 2, 7009-7032. Revere, M., Heermann, D. W. and Binder. K., 1988, Block density distribution function analysis of two-dimensional Lennard-Jones fluids. Europhys. Z.&t. 6, 585-590. Revere, M., Nielaba, P. and Binder, K., 1993, Simulation studies of gas-liquid transitions in two dimensions via a subsystem-block-density distribution analysis. 2. Phys. 90, 2 15-228. Rowlinson, J. S. and Widom, B.. 1982, Molecular Theory of CapilImity. Oxford University Press, Oxford. Smit, B., de Smedt, P. and Frenkel, D., 1989, Computer simulations in the Gibbs ensemble. Molec. Phys. 68, 93 l-950.

MANOJ MEHTA and DAVID A. KOFKE

2642

Valleau, J. P, 1991, Density-scaling: a new Monte Carlo technique in statistical mechanics. J. Camp. Phys. %, 193-216. APPENDIX In this section we provide details of the simulations performed in this study. We first review the method used to couple the volume changes, then we consider each integration technique, semigrand and osmotic. Volu?nQcoupling As discussed in the Introduction, the volumes of the two phases were coupled lo prevent unilateral condensation or vaporization of either of the phases (Kofke, 1993a). In lieu of independent sampling of V’ and p, we sample two new quantities, t, and tr. These define X and Y X = [l + exp(-

t,)]-‘;

Y=[l

+exp(-rty)]-l (AI)

which in turn define the densities of the two phases (Smit et al., 1989) 1-X @=1-_yP

PI +;

where p is a constant that lies between the expected average densities of the coexisting phases, but is otherwise arbitrary; in each simulation p was selected as the mean of the equilibrium densities of the preceding simulation in the series. Given N” for each phase, pa = NafVm gives directly the volumes. Equation (Al) ensures that X and Y lie between zero and unity (their use as surrogates for X and Y is not entirely necessary as X and Y normally would not sample values outside their domain of definition anyway). The transformation from Y’ and Y0 to tX and tr as the sampling variables introduces a Jacobian which must be considered when accepting proposed volume changes (Kofke. 1993a); it may he written J@

tr)

x9

=

N’N’

(p’- p)(p - p’) PP’P’(P’

-

P”)

.

(A3)

The quantities X and Y have a simple physical interpretation (Smit et nl., 1989). Consider a composite system of overall density p, comprising two systems of density p’ and p@, respectively. The liquid contains n’moles and occupies a volume Y’, with nr and e” defined likewise for the gas. X and Y then represent the fraction of the total particles and volume, respectively, that are in the liquid: X = n’/(n’+ n’); Y = u’/(u’+ u”). In performing the coupling, we are in effect letting the two simulation volumes represent larger systems that exchange matter when X varies (virtual particle exchange), and volume when Y varies. Except near a critical point, we have found that the density of each phase rarely, if ever, approaches the overall density p which it cannot unilaterally cross. So for the most part the coupling has little apparent effect on the sampling. Its role increases as the critical point is approached, and in fact here it serves to delineate, or even define, the two phases. Although artificial, such a construction is necessary for small systems. A system in coexistence exhibits two maxima in the distribution of densities, and these correspond to the two “phases”. If between the peaks there is a region of density for which the distribution is essentially zero (i.e. in a macroscopic system at any T < T,, or in a small system well removed from Tc), the phases may be defined unambigiously. Otherwise, we must introduce a somewhat arbitrary cutoff in density to separate and thereby define them. The introduction and choice of this cutoff density could complicate extrapolation to an infinite system, but we have not yet examined this notion. When applied to the study of other types of coexistence-those in which the density is not the primary distinguishing feature of the phases-the use of a cutoff may have a more deleterious effect. Coexisting solids and liquids for example may have overlapping density distributions, and

preventing each from sampling its full range of density could adversely affect the ensemble averages. Since the coupling mimics the sampling of densities used by the Gibbs ensemble, these considerations are relevant to those simulations as well.

In the semigrand method, each phase is simulated at constant T, P, N’ and t2. A flow diagram for a semigrand simulation algorithm is presented in Fig. Al. We sampled particle configurations in the usual way (Allen and Tildesley, 1987), and coupled volume sampling was performed as described above. The composition sampling required in a semigrand simulation was accomplished via particle identity changes: a box was chosen at random, and a molecule was picked randomly from that box; the identity of that particle was then changed, i.e. if it was a molecule of species 1, it was changed to 2, and uice versa. The proposed change was accepted with the probability given in Fig. Al (Kotke and Glandt, 1988). For most of the Gibbs-Duhem integration series, the integration step size AC2 was 0.05. Exceptions are as follows: in mixtures I and II we doubled and quadrupled this step to examine its effect on the results; in mixture III we used this step for tr --z0.8, then decreased it to 0.01 to integrate c2 from 0.8 to the critical value (approximately 0.98). In mixtures I and II, all simulations modeled each phase with 108 particles; this number was used also in the study of mixture III for c2 -c 0.8, but it was increased to 216 particles per phase in some simulations to describe the phases better as the critical point was approached. We also studied integration of the inverse derivative, &$Jalnp, for t2 ) 0.8 in mixture III; the results were consistent with the other semigrand simulations for this system. The initial configuration for each simulatibn in the series was taken as the final configuration of the previous run. As depicted in Fig. Al, each simulation was divided into cycles, where a cycle represented either: (a) N attempted particle displacements (where N is the total number of particles in both phases: N = N’ + Ng); (b) one coupled volume change; or(c) N attempted particle identity changes; the choice of (a), (h) or (c) was made with equal probability in a given cycle. In (a) and (c). the particles were chosen at random (as opposed to sequentially, or in some other deterministic way), while in (b) the choice between the two types of volume move (tX or t,,) was made with equal probability. Most simulations sampled 10,000 cycles, after an initial equilibration phase lasting 5ooO cycles; a few simulations of 20,000 production cycles were also performed with no significant change in the results. Osmotic ensemble In the osmotic ensemble, each phase is simulated at constant T, P, N; and f2. The flow diagram is presented in Fig. AZ. Osmotic simulations proceed much as semigrand ones do, with sampling of configurations, compositions, and densities. However, compositjon sampling is accomplished via insertion and deletion of species-2 particles, rather than the identity exchange used in semigrand simulation; the number of particles of species 1 is held fixed. It is advantageous to select as species 2 the component that is most easily deleted or inserted, as this leads to greater acceptance of these moves and hence better sampling. Care must be exercised in particle insertion/deletion. These changes, unlike the identity change moves in the semigrand simulations, affect the density and thus could undermine the purpose of the volume coupling outlined above. To bypass this problem we adjusted the volume when at tempting the insertion/deletion in such a way that the density remained unchanged. Thus a particle inscrtion/deletion was performed as follows: a box was chosen at random: it was decided with equal probability whether an insertion or deletion was to be attempted, and the volume change needed to keep the density of that box unchanged

4

START

initial&

I I

I attempt

I I

I I I I I

0

a coupled volume

exchange. perturb

h

n

or ty

e

C

accept with proWy (Ii)

I

Y C

I

I I

I e

I I 1

updateblocksum

e orobab i’ w (n)min[l.exp(-SACI)]

m = -1 for identity change from species 2 to species 1,

andm = +I otherwise

Fig. Al. Flow diagram describing the essential features of a simulation performed in the semigrand ensemble. Formulas (a)-(c) represent the acceptance probabilities for the Monte Carlo moves indicated in the figure. AU and A Vare. the change in internal energy and volume respectively, accompanying the move; V’ and N’ are the volume and particle number in phase a. p” = N’/Va, and the subscripts “old” and “new” indicate, respectively, values before and after the attempted move. The fugacity fraction of component 2 is 5 2, and p is the “overall” density chosen as part of the volume-coupling scheme; “min” indicates the minimum of the two arguments. Other variables arc as defined in the text. “Initialization” implies reading in the initial configuration, computing energies, etc; “relaxation cycles” proceed just as the cycles described explicitly in the figure (enclosed in the dashed-line box), except they do not contribute to the simulation averages; “blocksums” are used to estimate the statistical error in the simulation averages. Many variations of the algorithm can be used to produce correct results. For details, one should consult the excellent text of Allen and Tildesley (1987).

2644

MANOJ MEHTA and DAVID A. KOFKE

I

I I I I

I

I I I

3

I

altempt insuliddcl4on

;

of species 2 @ consuuu densit ad’uslvolumc

I I I

accept with pmWy (c)

I I I

I

NO

I I I I I

updateblocksum

anempts=N?

I I I

0

n e

C Y C

cS I I I

I

# cycles=blocksize 1 apply come=xOrto refine pressureestimate

I

for &leticm of species 2

Fig. A2. Flow diagrams describing the essential features of a simulations performed in the osmotic ensemble. Description follows that for Fig. Al. A, is the fugacity of species 2, and N; is the number of molecules of species 2 in phase = before the attempted insertion/deletion.

I e

Coexistence diagrams of mixtures by molecular simulation was then calculated. This hybrid move was accepted with probability given in Fig. A2. Each series was started the same way as in the semigrand runs, i.e. from the coexistence state.of one of the pure species; in all osmotic simulations, this was pure species 1. To propagate the series, steps in the fugacity of the species 2 were taken, starting withJz = 0. The entire range of composition cannot be spanned if we take-as in the semigrand series-the final configuration of one run as the initial configuration of the next. Because NT is fixed, species 1 can become dilute in a given simulation only through the addition of species-2 particles. Thus, as pure 2 is approached the system will become unacceptably large if no special precautions are taken (in fact, it becomes overspecified, with no independent extensive variable in the limit of pure 2). Accordingly, a “pre-processing” of the final configuration of a particular run was performed before using it to begin the next simulation. In this step, both phases were “re-scaled” by removing species-l particles from each. Specifically, the treatment proceeded as follows. First we set a target total number of particles N,.rD (summing both species in both phases). If the total number of particles at the end of a simulation did not exceed (10 + N,,,,), no treatment was applied. Otherwise, “surplus” particles of species 1 were deleted from each phase randomly, such that the number of species-l

2645

particles that remained in each phase was xtPN,,,s, where xf is the mole fraction of 1 in phase CLof the present configuration (this serves as an estimate of the mole fraction in the subsequent simulation), and U+ is the target fraction of particles in phase CI (we chose. 0 = 0.5 in each phase). To equilibrate the configuration, we then performed a short (5000-cycle) simulation in which no particle insertions or deletions were attempted. The resulting configuration served as the initial configuration of the next simulation of the series (these are the “process cycles” referred to in Fig. A2). The integration step size. A/, was on the order of 0.001 to 0.02, depending on the mixture. It was selected using knowledge of the pure-2 fugacity, but it could also have been chosen using Henry’s_constant data for species 2 in pure 1. With our choice, 10 to 20 simulations were required to span most of the composition range. Simulations were broken up into cycles, as in the semigrand series, except particle inse.rtionjdeletions steps were performed in lieu of the semigrand identity changes. Again, 10,000cycles were averaged, beyond an equilibration phase of 5000 cycles (not including the pre-processing cycles). Some runs averaged 20,000 cycles with no change in the results. The target number of particle was 216 (108 per phase) in most simulation; some simulations with a target of 512 total particles were performed for mixtures II and III.