J. Phys. Chem. B 2010, 114, 11515–11524
11515
Molecular Based Modeling of Associating Fluids via Calculation of Wertheim Cluster Integrals Hye Min Kim, Andrew J. Schultz, and David A. Kofke* Department of Chemical and Biological Engineering, UniVersity at Buffalo, The State UniVersity of New York, Buffalo, New York 14260-4200 ReceiVed: April 20, 2010
We examine a virial-like treatment for the equation of state of associating fluids defined in terms of a detailed molecular model. The approach implements Wertheim’s formulation of statistical thermodynamics of a classical fluid of molecules having strong directionally dependent association interactions. We employ the theory in its fundamental form, which expresses the pressure as an expansion in two or more aggregation densities, which themselves are related to each other by other series expansions. We employ Mayer-sampling Monte Carlo simulations to evaluate the cluster integrals defining the coefficients appearing in these series, yielding a multidensity virial-like equation of state, appropriate for the molecular system for which the cluster integrals were computed. We demonstrate this approach with a well-studied Lennard-Jones + association model, considering cases of atoms having one and two binding sites, respectively, and including all clusters involving up to four atoms. It is shown for this application that the Wertheim treatment is vastly superior to the standard (single-density) virial formulation, which fails in its description of the associating-fluid equation of state at very low densities. The Wertheim formulation for associating fluids is seen to be effective up to densities where the standard virial treatment to the same order begins to fail when applied to nonassociating fluids. 1. Introduction The virial equation of state (VEOS) is an appealing way to describe the thermodynamic behavior of real gases at low density:1
Z≡
p ) 1 + B2F + B3F2 + B4F3 + B5F4 + · · · FkT
(1) where Z is the compressibility factor, p is the pressure, F is the number density, k is the Boltzmann constant, and T is the absolute temperature. For practical application the series must of course be truncated, and we use VEOSn to indicate the virial series truncated after the nth virial coefficient. The nth virial coefficient (Bn) represents deviations from ideal-gas behavior arising in particular from interactions among n molecules, and it can be expressed in terms of cluster integrals over positions and orientations of n molecules.2 These cluster integrals can be conveniently expressed as diagrams, for example3
In each diagram, the black bonds represent the Mayer f function, fij ) exp(-βu(rij)) - 1 where uij is the intermolecular potential * To whom correspondence should be addressed. E-mail: kofke@ buffalo.edu.
Figure 1. Convergence of VEOS for Lennard-Jones system at kT/ε ) 1.4. Compressibility factor Z is plotted against density reduced by the critical density; Fcσ3 of the Lennard-Jones system is 0.317.34
between molecule i and molecule j with separation distance rij. The black points represent the integral of the molecule’s position over all space coordinates r and orientations ω;2,4 each molecule’s orientation integral is normalized to unity by the factor of Ω, as indicated. The great appeal of VEOS is the rigorous connection that it makes between the thermodynamic behavior and the molecular interactions, via these cluster integrals. In this manner it provides a viable route to bulk thermodynamics from first principles. VEOS performs well in describing the thermodynamic behavior of noncondensed nonpolar molecular fluids. For example, when applied to the Lennard-Jones system, VEOS4 is nearly converged at the critical density, as shown in Figure 1.5 In fact we find that the Lennard-Jones critical temperature from VEOS6 agrees very well with the literature value, as seen in Table 1, and at higher
10.1021/jp103573k 2010 American Chemical Society Published on Web 08/12/2010
11516
J. Phys. Chem. B, Vol. 114, No. 35, 2010
Kim et al.
TABLE 1: Critical Temperature According to VEOS Truncated after Various Orders n, for the Lennard-Jones5 and GCPM water7 Models
a
n
Tc,LJ (ε/k)
Tc,GCPM (K)
3 4 5 6 7 literaturea
1.4452693 1.299140 (14) 1.29048 (11) 1.3194 (10) 1.313 (10) 1.313 (1)
792 587 529
Literature data.
642
8,34
Figure 2. Convergence of VEOS for GCPM water at 673 K. Compressibility factor Z is plotted against density reduced by the critical density; Fc of GCPM is 0.3344 g/cm3.8
temperatures it is capable of describing the equation of state up to densities approaching the triple point.5 Even more interesting, at conditions of density where it has converged, VEOS provides a description of the equation of state that is more accurate than molecular simulations of thousands of atoms, which still have measurable inaccuracies due to finite-size effects.5 However, when applied to polar molecular fluids, VEOS is much less effective in describing thermodynamic properties. We recently reported6,7 values of virial coefficients up to B5 for the Gaussian charge polarizable model (GCPM),8 and we show in Figure 2 that in using these coefficients there is a significant deviation between VEOS4 and VEOS5 well below the critical density, indicating that the series is not converged at these conditions; Table 1 likewise shows that the critical temperature itself is poorly predicted. One approach to improving this comparison is to extend the VEOS to higher order, introducing coefficients beyond B5. The difficulty of computing these coefficients increases rapidly with each increment in order, particularly when accounting for the multibody effects inherent in polarizable models such as GCPM (and it is futile to describe gases using nonpolarizable models such as SPC that are formulated to describe the liquid9). Furthermore, we anticipate that convergence of the virial series will be slow for these systems, so an inordinate number of additional terms would be needed to reach the level of performance seen in the LennardJones application. Instead it is likely that a different clusterintegral approach that explicitly accounts for the effects of association would be more effective. The general problem of modeling associating fluids has a long history.10-17 Many treatments are based on a chemical model, in which physical association is treated as a reversible chemical
reaction among differently sized aggregates, which otherwise interact in a nonspecific way that can be described using conventional equations of state. Most models employ a cubic equation of state for this purpose, with interaggregate interactions described via simple 1-fluid composition-dependent mixing rules that determine the parameters of the equation of state. Models such as these obviously have only a qualitative connection to underlying molecular details that give rise to the thermodynamic behavior, and therefore they are inherently limited in their predictive capabilities. An important class of association models is formulated differently, in that they retain much more of the molecular detail, at least in principle. The most successful among these is the statistical associating fluid theory (SAFT) model and its derivatives.18,19 SAFT is built directly from the formalism and theory of associating fluids that was developed by Wertheim starting in 1984.20-22 There are two distinct parts of Wertheim’s work that are expressed in SAFT. The first is a diagrammatic treatment much like the Mayer theory embodied in VEOS but modified to take into explicit account the presence of strong directional interactions. Wertheim generated a very powerful theory of associating fluids by manipulating these reformulated cluster integrals to take advantage of the limits that steric interactions place on association. The second element of Wertheim’s work that appears in SAFT is a thermodynamic perturbation theory that expresses the nonspecific interactions in terms of a tractable molecular-based reference system. This element permits SAFT to be expressed in a relatively simple form that is suitable for practical applications. It also necessarily leads to some glossing over of the molecular details. In the present work we adopt just the first component of Wertheim’s theory, viewing it as a molecular-based alternative to VEOS for application to associating systems. In this manner we aim to formulate an equation of state that can apply to the associating system at densities well beyond the point where VEOS fails. Such a treatment has the key feature of retaining all molecular detail in the formulation of the model. That is, given a molecular model, one can expect the equation of state emerging from the treatment to describe the model’s thermodynamic behavior very accurately over appropriate conditions of density and temperature (and moreover it can be extended to mixtures with the same level of rigor). We expect of course ultimately to suffer the same limits exhibited by VEOS itself, namely, that the treatment will begin to fail at sufficiently high densities, and that it will not be directly applicable at all to the (subcritical) liquid phase. Part of the aim of this work is to explore the conditions where this happens for some simple model systems. In the next section, we review Wertheim theory as an alternative thermodynamic method for associating fluids. We develop the equation of state based on his theory, which maintains molecular-level detail but incorporates fluid association. In Section 3, the simple associating molecular models that we use to test the treatment are presented. After we review in Section 4 the Mayer sampling method for calculating cluster integrals, we present and discuss the simulation results in Section 5. In Section 6, we conclude and summarize. 2. Wertheim Equation of State (WEOS) Wertheim introduced the key ideas of his theory first through its development for systems that exhibit a single association site. Subsequently he went on to extend the theory to systems having multiple bonding sites. We review both of these formulations here. We will use “WEOS” to refer to treatments
Molecular Based Modeling of Associating Fluids
J. Phys. Chem. B, Vol. 114, No. 35, 2010 11517 more than one F bond connected to any point are excluded. This is a nontrivial consequence of the way the series is formulated (it does not hold, for example, in VEOS) and is a manifestation of the so-called steric incompatibility effect as introduced by Andersen.23,24 Given two mutually bonded molecules, it is assumed that the interactions are such that a third molecule cannot be bonded to either without introducing steric overlap with one or the other, and the diagrams are structured such that this overlap causes the configuration to contribute zero to the integral. As a consequence, some of the diagrams in the series have values that are very small, or identically zero, and so may be excluded from the formulation. The development leads to a function denoted c(0), expressed as the following series
1 1 1 c(0) ) W1F + W02F2 + W20F20 + W03F3 + 2 2 6 1 1 W F2F + (3Wa04 + 6Wb04 + Wc04)F4 + 2 21 0 24 1 (2Wa22 + 4Wb22 + Wc22 + Wd22)F20F2 + 4 1 (2Wa40 + 4Wb40 + Wc40)F40 + · · · 8 Figure 3. (a) Mayer f function; (b) decomposition of Mayer f function into reference part (fR) and association part (F). Both plots assume that molecules are oriented in a way that permits association interaction.
based on Wertheim’s development, with 1s-WEOS representing the one-site theory and 2s-WEOS written for the two-site formulation. 2.1. One-Association-Site Theory. The initial step in the Wertheim development20 (like others that preceded it23-25) is separation of the intermolecular potential into a reference potential uR(12) that has a weak (or no) attractive component and a short-ranged and highly directional association between binding sites uAA(12). Connected to this is a decomposition of the Mayer f function that appears in the usual cluster expansion, into a reference part (fR) and an association part (F), thus
f(12) ) exp[-βuR(12)] - 1 + exp[-βuR(12)] · {exp[-βuAA(12)] - 1} ) fR(12) + eRfAA(12) ) fR(12) + F(12)
(4)
We use W to indicate a cluster integral, with subscripts ij that correspond to the exponents of the associated density terms, Fi0Fj, and a letter superscript where needed to identify topologically distinct clusters of the same order ij. The diagrammatic representation of each term is listed in Table 2. Two other functions of interest are obtained from the functional differentiation of c(0) with respect to the densities and are given by the following series:
δc(0) 1 1 ) W11 + W102F + W103F2 + W121F20 + δF(1) 2 2 1 1 b1 b2 c1 3 a1 b1 (3Wa1 04 + 3W04 + 3W04 + W04)F + (2W22 + 2W22 + 6 2 c1 d1 2 (5) 2Wb2 22 + W22 + W22)F0F + · · ·
c0(1) )
δc(0) 1 b3 ) W120F0 + W221F0F + (2Wa2 22 + 2W22 + δF0(1) 2 1 d2 a1 b1 b2 2 + Wc2 22 + W22)F0F + (2W40 + 2W40 + 2W40 + 2 3 (6) Wc1 40)F0 + · · ·
c1(1) )
(3)
We can see that how the f bond is broken into an F bond and an fR bond in Figure 3. The other difference between VEOS and 1s-WEOS is the density variables. In VEOS, only one kind of density variable, the total density (F), appears in the equation of state. In 1sWEOS, two kinds of density variable enter, one again being the total density (F), and in addition to this the monomer density (F0) appears. Derivation of 1s-WEOS starts with the activity expansion of the grand-canonical free energy. Wertheim substitutes for the Mayer f function its sum of fR and F. Resummation of the series is performed to yield a new doubly connected series in F and F0. In the new series a point has weight F if is a “monomer”, meaning it has no F-bonds joining it to another point; otherwise, it has weight F0. Given that we have only a single association site per molecule, any diagrams having
2Wb4 22
The coefficients in these series are specified in a manner similar to that for c(0), with numeric superscripts appended to indicate that these diagrams have a root point and further to identify topologically distinct diagrams obtained from differentiation of a common term in the c(0) series. These are listed also in Table 2. The second of these provides a relation between F and F0:20
F ) F0(1 + c1)
(7)
where we suppress the (1) argument on c1 (and the densities) because we consider a homogeneous system. Finally, the expression for the pressure is given by:20
11518
J. Phys. Chem. B, Vol. 114, No. 35, 2010
TABLE 2: Wertheim Coefficients and Values Calculated from MSMC for 1s-WEOS (and Which Appear Also in 2s-WEOS)a
Kim et al. Equations 6, 7, and 8 are the primary working equations of the 1s-WEOS theory. Once all of the diagrams to a given order are computed, then for a given overall density F we can compute F0 from eq 7 and the total pressure from eq 8 (it is worth noting that, if there is no association and F(12) ≡ 0, the coefficients of the terms involving F0 vanish and the equation reduces to VEOS for the reference potential). We denote 1s-WEOSn as the result of this calculation with the series 7 and 8 truncated after terms with coefficients given by diagrams containing n points. 2.2. Two-Association-Site Theory. The derivation of 2sWEOS22 is similar to 1s-WEOS, and we focus here on the differences. In 2s-WEOS, each molecule has two association sites, labeled A and B. Unlike 1s-WEOS, chain aggregates longer than dimers appear in 2s-WEOS, further admitting the possibility of ring formation. Presently we exclude ring diagrams by stipulating that the angle between sites A and B does not permit these structures to form, which is valid for the model studied here for the diagrams we consider, but in general this cannot be assumed. The Mayer f function is decomposed as
f(12) ) fR(12) + eR{fAB(12) + fBA(12) + fAB(12) · fBA(12)} (9) where the subscript AB indicates the A site on molecule 1 bonding to the B site on molecule 2, and vice versa for subscript BA. Bonding interactions are defined only between unlike sites, so there are no fAA and fBB terms in the decomposition. Here we exclude double bonding formation (a 2-molecule ring), so either (or both) fAB(12) or fBA(12) must be zero for a given configuration, and the decomposition of the Mayer f function is simplified to
f(12) ) fR(12) + eR{fAB(12) + fBA(12)} ) fR(12) + FAB(12) + FBA(12) (10) As for 1s-WEOS, we arrive at a function c(0), but in this case having diagrams that allow two association bonds to a site and with density weights assigned to the points as follows. With each point having association bonds at sites R, we associate a sum of number density ∑γ⊂(A,B-R)Fγ(1).22 For example, if the point has association bonds at site A, we weight the point as ∑γ⊂BFγ(1), which is F0 + FB. Wertheim denoted this sum of number densities as σR(1). In 2s-WEOS, possible σR(1) factors are given by a In the diagrams, the fR bond is a solid line colored as red, and the F bond is a hatched line colored as blue. Numbers in parentheses indicate the 67% confidence limits of the last digit(s) of the given value. Diagrams without root points (2nd column) should be divided by V when ascribing the tabulated values to them.
βp ) F0 -
∫ F(1)c0(1)d(1) + c(0)
1 1 1 1 ) F0 - W02F2 + W20F20 - W03F3 - (3Wa04 + 2 2 3 8 1 6Wb04 + Wc04)F4 - (2Wa22 + 4Wb22 + Wc22 + 4 1 Wd22)F20F2 + (2Wa40 + 4Wb40 + Wc40)F40 + · · · (8) 8
σ0(1) ) F0(1) σA(1) ) F0(1) + FA(1) σB(1) ) F0(1) + FB(1) ) σA(1) σAB(1) ) F0(1) + FA(1) + FB(1) + FAB(1) ) F(1)
(11) Since in our work we allow association only between site A and site B, FA (1) and FB (1) are necessarily equal, and thus so are σA(1) and σB(1). Inasmuch as the value assigned to a point depends on the type of bonds it has joining it to other points, it is necessary to maintain the distinction between FAB and FBA. However, in all cluster series that we encounter here, every diagram containing an FAB bond between a pair of points will have a
Molecular Based Modeling of Associating Fluids TABLE 3: Same as Table 1, but Showing Diagrams Appearing Only in 2s-WEOS
J. Phys. Chem. B, Vol. 114, No. 35, 2010 11519 that we accord this diagram will represent the sum of these variants. Further, we note that, when interpreted this way, all of the diagrams appearing in 1s-WEOS also appear in 2s-WEOS. For the models we use, the values of these diagrams appearing in 2s-WEOS are simple multiples of their values in 1s-WEOS; the multiplier is the number of ways that all F bonds can be exchanged for FAB or FBA bonds, and it is listed with each diagram in Table 2. For the 1sWEOS diagrams appearing in 2s-WEOS, in the subscript ij on the coefficient, i is the sum of exponents on σA and σB (which will have equal exponents), and j is the exponent on σAB (these diagrams do not multiply powers of σ0). The diagrams only used in 2s-WEOS are listed in Table 3 and include diagrams in which two F bonds join a single point. Here it should be understood that these bonds again represent FAB or FBA, but when making the substitutions it must result that both types join a given point; for example, a diagram cannot have a point joined by two FAB bonds. The values we report for these diagrams are the sum over all such substitutions, all of which have a common set of σ weights distributed among the points. The coefficient Wijk multiplies k . the term σ0i σAj/2σBj/2σAB We can now proceed with the series that appear in 2s-WEOS. The result for c(0) is expressed as
1 1 1 2 3 c(0) ) W1σAB + W02σAB + W20σAσB + W03σAB + 2 2 6 1 1 1 W σ σ σ + W120σ0σAσB + (3Wa04 + 6Wb04 + 2 21 A B AB 2 24 1 4 2 + (2Wa22 + 4Wb22 + Wc22 + Wd22)σAσABσAB + Wc04)σAB 4 1 1 (2Wa40 + 4Wb40 + Wc40)σA2σB2 + (Wa121 + 2Wb121 + 8 2 1 Wc121 + Wd121 + We121)σ0σAσBσAB + (Wa220 + Wb220 + 2 2Wc220 + Wd220)σ20σAσB + · · · (12) Again cR(1) functions are obtained from c(0) by functional differentiation. In 2s-WEOS, c(0) is differentiated with respect to σR(1) instead of density variables.22
cR(1) )
δc(0) δσΓ-R(1)
(13)
The glossing over of how the σ weights are assigned to the points obscures the process slightly, but once the variants are considered and recombined, the outcome can still be expressed using the generic F bond with σ weights assigned to the diagram rather than the points. The results are as follows, given in terms of four functions c0(1), cA(1), cB(1), and cAB(1).
δc(0) 1 2 + ) W11 + W102σAB + W103σAB δσAB(1) 2 1 1 1 b1 b2 c1 3 W σ σ + (3Wa1 04 + 3W04 + 3W04 + W04)σAB + 2 21 A B 6 1 b1 b2 c1 d1 (2Wa1 22 + 2W22 + 2W22 + W22 + W22)σAσBσAB + 2 1 a1 c1 d1 e1 (W + 2Wb1 121 + W121 + W121 + W121)σ0σAσB + · · · 2 121 (14)
c0(1) ) counterpart in the series with that bond replaced by an FBA bond, perhaps with the σ weights reshuffled among the points; and, at least in the models examined here, these clusters will have the same numerical value. Accordingly, we will introduce an economy to our presentation and represent these variants by a single diagram with the different FAB and FBA bonds represented by a single type of F bond, and the value
11520
J. Phys. Chem. B, Vol. 114, No. 35, 2010
Kim et al.
δc(0) 1 1 ) W120σA + W221σAσAB + δσB(1) 2 2 1 1 1 b3 b4 c2 W σ σ + (2Wa2 22 + 2W22 + 2W22 + W22 + 2 120 0 A 4 1 2 a1 b1 b2 c1 2 Wd2 22)σAσAB + (2W40 + 2W40 + 2W40 + W40)σAσB + 4 1 a2 b3 c2 d2 (W + Wb2 121 + W121 + W121 + W121 + 2 121 1 a1 b1 c1 c2 We2 121)σ0σAσAB + (W220 + W220 + W220 + W220 + 2 2 Wd1 220)σ0σA + · · · ) cB(1) (15)
cA(1) )
cAB(1) )
Figure 4. One-association-site model.
δc(0) 1 1 + 2Wb4 ) W2120σAσB + (Wa3 121 + δσ0(1) 2 2 121
d3 e3 a2 b2 c3 Wc3 121 + W121 + W121)σAσBσAB + (W220 + W220 + W220 + d2 Wc4 220 + W220)σ0σAσB + · · ·
(16)
The total density F(12) is composed of Figure 5. Two-association-site model.
F(12) ) F0(12) + FA(12) + FB(12) + FAB(12)
(17) where F0(12) is monomer density, FR(12) is the density of particles bonded only at site R, and FAB(12) is the density of particles bonded at sites A and site B simultaneously. Each density variable can be expressed in terms of cR(1) functions22
FA(1) ) F0(1)cA(1) ) FB(1)
(18)
FAB(1) ) F0(1)[cAB(1) + cA(1)cB(1)]
(19)
The pressure for 2s-WEOS is given in terms of cluster integrals and σR(1) factors.22
βp ) F -
it is significantly more effective than VEOS in doing so. Presently we restrict our attention to abstract models that have well-defined association conditions, and once we have established the effectiveness for these simple cases we can in subsequent work consider how the approach can be applied to more practically relevant systems such as realistic models for water. Gubbins and co-workers followed a similar path in the development of SAFT, and we can apply here the same molecular models they used in their work. In this study, we examine a one-association-site model (“1site model”) (Figure 5) and a two-association-site model (“2site model”) (Figure 6). The potentials are defined by sum of the Lennard-Jones (LJ) potential and site-site association.
u(12) ) uLJ(12) +
∫ [σABc0(1) + σBcA(1) + σAcB(1) +
1 1 2 σ0cAB(1)]d(1) + c(0) ) F - W02σAB - W20σAσB 2 2 1 1 W σ3 - W21σAσBσAB - W120σ0σAσB - (3Wa04 + 3 03 AB 8 3 4 - (2Wa22 + 4Wb22 + Wc22 + 6Wb04 + Wc04)σAB 4 3 2 - (2Wa40 + 4Wb40 + Wc40)σA2σB2 Wd22)σAσABσAB 8 3 a (W + 2Wb121 + Wc121 + Wd121 + We121)σ0σAσBσAB 2 121 3 a (W + Wb220 + 2Wc220 + Wd220)σ20σAσB + · · · (20) 2 220 As with the 1-site model, given the total density F and values for all cluster integrals, we can compute the other densities via 18 and 19, the σ from 11, and then the pressure from 20. 3. Association Models Our aim in this work is to develop and demonstrate the capability of WEOS to describe associating systems when applied using explicit calculation of the cluster integrals that appear in its formulation, and in particular to examine whether
∑ ∑ uAB(12) A
(21)
B
The LJ potential is given by
[( σr ) - ( σr ) ]
uLJ(12) ) 4ε
12
6
(22)
where r is the separation distance between two molecules and ε and σ are the LJ energy and size parameters. For the strong site-site association, a conical square-well potential ua(12) is used.
ua(12) )
{
-εa if r < reff, θA1 < θeff, and θB2 < θeff 0, otherwise (23)
where θRn is the maximum angle between a line from the center of molecule 1 to molecule 2 and a line from a site on molecule n to the center of molecule n (see Figures 4 and 5); reff is the maximum separation distance for association, and θeff is the maximum angle to make a bond between a site on molecule 1 and a site on molecule 2. In this work the association energy εa
Molecular Based Modeling of Associating Fluids
J. Phys. Chem. B, Vol. 114, No. 35, 2010 11521
is set as -16ε, and the parameters reff and θeff are set to 1.0σ and 27°, respectively. 4. Methods 4.1. Mayer Sampling. Mayer sampling Monte Carlo (MSMC) simulation26 is used to evaluate cluster integrals appearing in VEOS and WEOS. MSMC employs methods developed for free energy calculation, and in particular in this work we employ overlap sampling in the optimum form given by Bennett.27-29 Through overlap sampling, a ratio of the cluster integral of interest (the target system) to a known cluster integral (the reference system) is evaluated.
Γ(T) ) Γ0
〈γ/π〉π /〈γOS /π〉π 〈γ0 /π0〉π0 /〈γOS /π0〉π0
(24)
where Γ is a cluster integral that we want to calculate, γ is the integrand or sum of integrands in Γ, and π is the distribution function for weighting the configurations; the subscript 0 represents a value for the reference system, and the angle brackets specify an ensemble average weighted by the subscripted distribution functions. For an importance-sampling formulation, π is set as the absolute value of γ. γOS represents the overlap function between the target and reference system as
γOS )
|γ0 ||γ| R|γ0 | + |γ|
(25)
where R is a parameter optimized for the convergence of the calculation. For all n-particle clusters appearing in 1s-WEOS (i.e., those in Table 2), the reference system is defined as the sum of clusters for the nth virial coefficient of hard spheres of diameter 1.5σ. Thus in this reference, associating configurations are given no extra weight. Diagrams appearing also in 2s-WEOS (i.e., those in Table 3) are handled differently, because ignoring the association bonds in the reference is considerably more inefficient. Instead in these cases we use a reference formed as an associating chain. Each F-bond in the target is replaced by a hard association bond (Fhard) that has the effect of forcing the joined pair to be associated (in other words, a bond that is zero if the association condition in eq 23 is not satisfied, and unity if it is). If a point in the target is not connected by any F bonds (this is encountered only in the diagrams in W121), then it is joined to the end of the chain by a hard-sphere (diameter 1.0σ) f bond. All remaining pairs of points are then joined by 1.0σ hard-sphere e bonds. The reference chain diagrams before the addition of the e bonds can be evaluated analytically: for a chain of n molecules, the integral over positions and orientations of molecules 2 through n is 2 × ((4/3)πσ3)n-1 × ((1/2)(1 - cos
θeff))2(n-1); for the W121 reference another factor of (4/3)πσ3 is multiplied for the point joined by the hard-sphere f bond. We use this cluster as a reference for a simple Mayer-sampling calculation that gives us the desired reference, with the e bonds included. The values of these hard-chain reference systems are listed in Table 4. We calculated virial coefficients and Wertheim coefficients up to the fourth order for both one- and two-association-site models. One billion MSMC samples were conducted in calculating clusters for the 1-site model, while 1010 samples were used for the 2-site model. For error analysis, 10 independent simulations of these lengths were performed for two- and threemolecule coefficients, and 20 independent runs were used for the four-molecule coefficients. Thus the total number of MSMC samples used to determine the each coefficient ranged from 1010 to 2 × 1011, depending on the number of molecules and binding sites. We collected these results at two temperatures, kT/ε equal to 1.5 and 2.0, respectively. The ratio of the amount of sampling of the reference system to the target system varied depending on the diagramssthe simulation effort can be optimized by adjusting this ratio such that the increment in uncertainty from the two contributions is equal. In 1s-WEOS (i.e., those in Table 2), most diagrams were computed with 4 to 19% of total sampling in the reference system. In most diagram calculations, the main contribution to the uncertainty came from the target system; thus most sampling was performed in the target system to reduce the uncertainty. However, if all of the points in a 1sa WEOS diagram are joined by an F bond (e.g., W40 ), then between 38 to 71% of total sampling occurred in the reference system. This is explained by the fact that that the reference system needs a long time to make one particle bind with another successfully. In the calculation of diagrams involved in 2sWEOS (i.e., those in Table 3), most diagrams were computed with 57 to 88% of total sampling in the reference system. In the target system, the Lennard-Jones potential does not allow an atom within a chain to approach to within about 0.96σ of the next atom in the chain. However, in the chain reference system, the adjacent atoms have no such repulsion, and so they can overlap substantially; thus, the volume of phase space of the chain reference system is considerably larger than the target system. Consequently the main contribution to the uncertainty came from the chain reference system and the sampling spent much time in the reference system compared to the target system. Some simple tweaking of the reference could have been done to make this more efficient (e.g., imposing a minimum separation distance), but we did not pursue such modifications in this work. 4.2. Isothermal-Isobaric (NPT) Ensemble Simulation. To provide a basis for gauging the performance of the VEOS and WEOS models, we collected equation of state data via conventional Monte Carlo (MC) simulation in the isothermal-isobaric (NPT) ensemble. The simulated systems for 1-site and 2-site models consisted of N ) 256 molecules for
TABLE 4: Cluster Integrals and Values for Hard Chain Reference System Used to Compute Coefficients in 2s-WEOSa
a
b
The black hatched bond represents an Fhard bond, and the black dashed bond represents an e-bond for hard spheres of diameter 1.0σ. Analytic results. c Mayer sampling results.
11522
J. Phys. Chem. B, Vol. 114, No. 35, 2010
pressures less than 0.05 and 1024 molecules for pressures greater than 0.05. The truncation distance of the intermolecular potential was set as 6.0σ, and we applied standard long-range corrections. We performed 2 × 1011 MC steps for each simulation. Ten independent simulations were performed for error analysis, with each simulation initialized to an fcc lattice. The initial configurations were equilibrated for 10% of total MC steps for the 1-site model and 20% for the 2-site model. In our NPT data, the standard deviation of the mean is less than about 0.1%. To aid in the convergence of the simulation averages, we introduced MC trials that implement the unbonding-bonding (UB) association bias algorithm.30 One of the main reasons for difficulty in simulating the associating system is that the region of the interaction between a site on molecule 1 and a site on molecule 2 that satisfy the criteria of the association is quite small. Therefore, the probability of making a bond with conventional trial moves is very low. Further, once the association has occurred, the probability to break this binding is also low, because of the accompanying energy penalty. The UB algorithm assigns an association-bias interaction region to each molecule, and with it the generation and breaking of association interactions occur with enhanced probability. In our simulations, the bonding region for the UB bias algorithm was assigned as a cone with the angle between the orientations at 27°.31 We explicitly exclude multiple bonding interactions within the NPT simulation. If two molecules approach in the same association site area simultaneously, the energy is taken to be infinite, and we reject this trial. Also any volume change trial that results in the formation or deletion of a new association bond is rejected (this rejection of any trial that causes the change of association is applied to not only expansion of the box but also compression because of the reversibility). Trials are also included in which molecular aggregates are rigidly translated or rotated; such moves are rejected in the rare cases in which they generate a new association bond. 4.3. SAFT. It is worthwhile to compare WEOS and the simulation data to the established SAFT treatment. In our application of SAFT, the compressibility factor of the 1-site model is given by32
Z1-site ) ZLJ + Zassociation ) (1 + B2,LJF + B3,LJF2 + B4,LJF3 + · · · ) + ∂X0 1 1 (26) F ∂F T,N X0 2
( )(
)
where the fraction of monomers X0 ) (-1 + (1 + 4F∆)1/2)/ 2F∆, with
∆ ) π(1 - cos θeff)2[exp(-εa /kT) - 1]
∫0r
eff
gLJ(r)r2dr
The contribution ZLJ is just the compressibility factor for the Lennard-Jones core, and we choose here to describe it using VEOS7. In this respect the SAFT treatment is identical to the WEOS formulation. The association contribution is where the difference lies. A required input here is the radial distribution function of Lennard-Jones potential, gLJ(r). In routine applications of SAFT this is sometimes handled using perturbation about hard spheres or by looking up tabulated values. In the present case we evaluated it (or more precisely, the integral in which it appears) via MC simulation of the nonassociating LJ
Kim et al.
Figure 6. Pressure of the 1-site model with association strength εa ) -16ε at reduced temperature (a) 1.5 and (b) 2.0 from VEOS (black lines), WEOS (red lines), NPT simulation (blue circles), and SAFT (black triangles). All quantities are in units of Lennard-Jones σ and ε.
system. For these simulations the number of atoms N was 256 for simulating systems of density less than 0.1 and 1024 for systems of density greater than 0.1. The truncation distance of the potential was set to 6σ, and we applied standard long-range corrections. The derivative appearing in eq 26 was evaluated numerically. Formulas for the 2-site model are very similar32
Z2-site ) (1 + B2,LJF + B3,LJF2 + B4,LJF3 + · · · ) + ∂XA 1 1 2F ∂F T,N XA 2
( )(
)
(27)
where the fraction of particles not bonded at site A is XA ) (-1 + (1 + 4F∆)1/2)/2F∆ and the fraction of monomers is X0 ) XA2. 5. Results and Discussion 5.1. 1-Site Model. Figure 6, parts a and b, displays the pressure of the 1-site model as a function of total density at temperatures 1.5 and 2.0, respectively. Shown in the plot are pressures given by VEOS and WEOS, in comparison to NPT MC simulation data. The simulation results for T ) 1.5 were not collected beyond the density 0.1. In this regime the simulations could not maintain the gas phase density, finding
Molecular Based Modeling of Associating Fluids instead that it would condense to the liquid, indicating that the fluid is in a metastable or unstable state beyond the density 0.1. The isotherm for T ) 2.0 appears to be supercritical, and we collected data up to a density of 0.6. Also, although we performed calculations for VEOS4 we do not include them in the plots because the uncertainty of B4 of the association site models is more than 100% of the value, and the plot of VEOS4 is diverging highly from other VEOS. It is clear in both figures that there are significant deviations between the VEOS and the simulation data at all but the lowest densities. In contrast, the application of even the lowest-order WEOS shows much better convergence than all of the VEOS. There is moreover a very good agreement between the WEOS and the NPT simulation results. At temperature 1.5, WEOS is converging at third order, and the addition of the fourth-order contribution is very small. Also, the precision of WEOS is such that the standard deviation in the calculated pressure is less than 1%. The figure also shows that pressures calculated by SAFT agree well with WEOS and the NPT simulation data. This is not too surprising, given the common origins with WEOS in the Wertheim formalism. We find that VEOS at T ) 2.0 is slightly more effective than that observed at the lower temperature. This might be expected, inasmuch as the increase in temperature diminishes the importance of association. Still, overall, the conclusions are unchanged from the lower-temperature results. For densities below 0.05, the pressures from WEOS agree very well with VEOS, which appears to be converged. However, beyond that density the VEOS show significant deviation from each other and from the NPT data, while WEOS is again proving effective to a much higher density. WEOS4 agrees well with the NPT simulation to the densities of about 0.5. Beyond density 0.3, the contribution of the fourth-order term in WEOS is quite important. Also we find in this range of densities that the NPT data agree with WEOS4 slightly better than the SAFT results, although upon approaching a density of about 0.6 WEOS4 loses accuracy and SAFT provides the better comparison to simulation. Overall, the results tell us that the thermodynamic perturbation theory underlying SAFT, and the density expansion inherent in WEOS, are both good treatments for the conditions considered here and, more importantly, that the Wertheim formalism is effective at capturing the complications due to association. Wertheim’s treatment also allows for the prediction of the degree of association, which can be characterized as the fraction of monomers (F0/F), computed via eqs 6 and 7. In Figure 7, we compare the monomer fraction from WEOS to averages observed in the NPT simulation. The agreement between the Wertheim theory and the NPT simulation is excellent, and also SAFT data agree well with the WEOS approach. 5.2. 2-Site Model. Turning now to the 2-site model, in Figure 8 we again see that all of the VEOS show significant deviations compared to WEOS. At T ) 1.5, the fourth order of WEOS still makes an important contribution. Again, the NPT data are not available because spontaneous condensation suggests that the fluid is in a metastable or unstable state beyond the density 0.02. In Figure 8b, we have isotherms for T ) 2.0, and it appears that this temperature is near the critical value. The overall picture is as found in the other cases. NPT simulation results seem to agree better with WEOS3 than WEOS4 at high density region, but this must certainly be a coincidence given that WEOS4 describes the simulation data better than WEOS3 at lower density. For WEOS to predict the true behavior beyond density 0.3, higher order terms than the fourth order are required. The
J. Phys. Chem. B, Vol. 114, No. 35, 2010 11523
Figure 7. Fraction of monomers of the 1-site model with association strength εa ) -16ε at reduced temperatures 1.5 and 2.0 from WEOS (red lines), NPT simulation (blue circles), and SAFT (black triangles). All quantities are in units of Lennard-Jones σ and ε.
Figure 8. Pressure of 2-site model with association strength εAB ) -16ε at reduced temperature (a) 1.5 and (b) from VEOS (black lines), WEOS (red lines), NPT simulation (blue circles), and SAFT (black triangles). All quantities are in units of Lennard-Jones σ and ε.
isotherms moreover suggest that 0.3 is roughly the critical density, which in previous studies was found to be about the upper limit of density for which VEOS4 is effective in describing nonassociating fluids.33 Given that VEOS is used to describe the reference system, one cannot expect anything better from WEOS. What is interesting then is that WEOS is effective in treating the effects of association, allowing the cluster series to work to the point where the nonassociation characterization displays its limits.
11524
J. Phys. Chem. B, Vol. 114, No. 35, 2010
Figure 9. Fraction of monomers of the 2-site model with association strength εa ) -16ε at reduced temperatures 1.5 and 2.0 from WEOS (red lines) and NPT simulation (blue circles) and SAFT (black triangles). All quantities are in units of Lennard-Jones σ and ε.
In Figure 9, the fraction of monomers for the 2-site model using WEOS4 is showing excellent agreement with the NPT simulation and SAFT results. 6. Conclusions Strong intermolecular association significantly reduces the effectiveness of the VEOS. An extension in the same spirit as VEOS but which accounts for the effects of association is found in Wertheim’s multidensity formalism. The approach permits the formulation of an equation of statesin the form of several coupled algebraic equationssfor a detailed molecular model that is accurate to the model, within the limits of convergence of the series. The model coefficients are expressed in the form of graphs representing cluster integrals that can be computed by the MSMC method. In this work we presented the graphs needed to compute the coefficients for 1-site and 2-site association models, up to fourth order in the density, and we computed these coefficients for a simple association model. In this application WEOS displays a much better convergence than the conventional VEOS, and a comparison to conventional Monte Carlo simulation results for the same model further confirms the effectiveness of the Wertheim approach. The results suggest that the Wertheim treatment permits the computational cluster-integral modeling approach to remain effective for densities up to those where VEOS begins to fail when applied to nonassociating systems. If this is the case, it suggests that it will be practical to characterize associating systems up to densities that encompass the critical point. Applications to realistic models such as GCPM water will require extension of the Wertheim formalism to multibody, nonpairwise-additive molecular models. Acknowledgment. This work is supported by the National Science Foundation under Grant CBET-0854340. Calculations were performed using resources from the University at Buffalo Center for Computational Research. References and Notes (1) Mason, E. A.; Spurling, T. H. The Virial Equation of State. In The International Encyclopedia of Physical Chemistry and Chemical Physics, Topic 10: The Fluid State; Pergamon Press: Elmsford, NY, 1969; Vol. 2. (2) McQuarrie, D. A. Statistical Mechanics; University Science Books: Herndon, VA, 1976.
Kim et al. (3) Mayer, J. E.; Mayer, M. G. Statistical mechanics; Wiley: New York, 1940. (4) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids, 3rd ed.; Academic Press: New York, 2006. (5) Schultz, A. J.; Kofke, D. A. Sixth, seventh and eighth virial coefficients of the Lennard-Jones model. Mol. Phys. 2009, 107 (21), 2309– 2318. (6) Benjamin, K. M.; Schultz, A. J.; Kofke, D. A. Virial Coefficients of Polarizable Water: Applications to Thermodynamic Properties and Molecular Clustering. J. Phys. Chem. C 2007, 111 (43), 16021–16027. (7) Benjamin, K. M.; Schultz, A. J.; Kofke, D. A. Fourth and fifth virial coefficients of polarizable water. J Phys. Chem. B 2009, 113 (22), 7810–7815. (8) Paricaud, P.; Predota, M.; Chialvo, A. A.; Cummings, P. T. From dimer to condensed phases at extreme conditions: Accurate predictions of the properties of water by a Gaussian charge polarizable model. J. Chem. Phys. 2005, 122 (24), 244511. (9) Benjamin, K. M.; Schultz, A. J.; Kofke, D. A. Gas-Phase Molecular Clustering of TIP4P and SPC/E Water Models from Higher-Order Virial Coefficients. Ind. Eng. Chem. Res. 2006, 45 (16), 5566–5573. (10) Heidemann, R. A.; Prausnitz, J. M. A van der Waals-type equation of state for fluids with associating molecules. Proc. Natl. Acad. Sci. U.S.A. 1976, 73 (6), 1773–6. (11) Green, D. G.; Jackson, G. Theory of phase equilibria for model aqueous solutions of chain molecules: water + alkane mixtures. J. Chem. Soc., Faraday Trans. 1992, 88 (10), 1395–409. (12) Visco, D. P., Jr.; Kofke, D. A. Modeling the Monte Carlo simulation of associating fluids. J. Chem. Phys. 1999, 110 (12), 5493–5502. (13) Economou, I. G.; Donohue, M. D. Chemical, quasi-chemical and perturbation theories for associating fluids. AIChE J. 1991, 37 (12), 1875– 94. (14) Anderko, A. Association and semiempirical equations of state. J. Chem. Soc., Faraday Trans. 1990, 86 (16), 2823–30. (15) Economou, I. G.; Donohue, M. D. Equations of state for hydrogen bonding systems. Fluid Phase Equilib. 1996, 116, 518–529. (16) Visco, D. P.; Kofke, D. A. Improved thermodynamic equation of state for hydrogen fluoride. Ind. Eng. Chem. Res. 1999, 38 (1), 4125–4129. (17) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria; Prentice-Hall: Englewood Cliffs, NJ, 1986. (18) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. SAFT: equation-of-state solution model for associating fluids. Fluid Phase Equilib. 1989, 52, 31–8. (19) Muller, E. A.; Gubbins, K. E. Molecular-based equations of state for associating fluids: A review of SAFT and related approaches. Ind. Eng. Chem. Res. 2001, 40 (10), 2193–2211. (20) Wertheim, M. Fluids with highly directional attractive forces. I. Statistical thermodynamics. J. Stat. Phys. 1984, 35, 19–34. (21) Wertheim, M. Fluids with highly directional attractive forces. II. Thermodynamic perturbation theory and integral equations. J. Stat. Phys. 1984, 35, 35–47. (22) Wertheim, M., Fluids with highly directional attractive forces. III. Multiple attraction sites. J. Stat. Phys. 1986, 42, 459–491. (23) Andersen, H. C. Cluster expansions for hydrogen-bonded fluids I. Molecular association in dilute gases. J. Chem. Phys. 1973, 59, 4717. (24) Andersen, H. C. Cluster expansions for hydrogen bonded fluids. II. Dense liquids. J. Chem. Phys. 1974, 61, 4985. (25) Lockett, A. M., III. A theory of homogeneous condensation from small nuclei. I. Modified Mayer theory of physical clusters. J. Chem. Phys. 1980, 72 (9), 4822–31. (26) Singh, J. K.; Kofke, D. A. Mayer Sampling: Calculation of Cluster Integrals using Free-Energy Perturbation Methods. Phys. ReV. Lett. 2004, 92 (22), 220601/1–220601/4. (27) Kofke, D. A. Free energy methods in molecular simulation. Fluid Phase Equilib. 2005, 228-229, 41–48. (28) Bennett, C. H. Efficient Estimation of Free Energy Differences from Monte Carlo Data. J. Comput. Phys. 1976, 22, 245–268. (29) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: New York, 1996. (30) Wierzchowski, S.; Kofke, D. A. A general-purpose biasing scheme for Monte Carlo simulation of associating fluids. J. Chem. Phys. 2001, 114 (20), 8752–8762. (31) Tsangaris, D. M.; de Pablo, J. J. Bond-bias simulation of phase equilibria for strongly associating fluids. J. Chem. Phys. 1994, 101 (2), 1477– 89. (32) Ghonasgi, D.; Chapman, W. G. Theory and simulation for associating hard chain fluids. Mol. Phys. 1994, 83 (1), 145–58. (33) Schultz, A. J.; Kofke, D. A. Virial coefficients of model alkanes. J. Chem. Phys. 2010, submitted. (34) Perez-Pellitero, J.; Ungerer, P.; Orkoulas, G.; Mackie, A. D. Critical point estimation of the Lennard-Jones pure fluid and binary mixtures. J. Chem. Phys. 2006, 125 (5), 054515/1-054515/9.
JP103573K
