Ind. Eng. Chem. Res. 1997, 36, 4041-4051
4041
Using Molecular Orbital Calculations To Describe the Phase Behavior of Hydrogen-Bonding Fluids† Jeffrey P. Wolbach and Stanley I. Sandler* Center for Molecular and Engineering Thermodynamics, Department of Chemical Engineering, University of Delaware, Newark, Delaware 19716
We have used Hartree-Fock theory and density functional theory to compute the enthalpy and entropy changes of dimerization for water, methanol, and the family of carboxylic acids. These results are used in a physical equation of state, the statistical associating fluid theory (SAFT), in order to model the phase behavior of these hydrogen-bonding compounds. A procedure has been developed to relate the calculated enthalpy and entropy changes to the association parameters in SAFT using only low-pressure data, as well as to relate molar volumes from molecular orbital calculations to the segment size and chain length parameters in SAFT. By doing so, the SAFT model is reduced to a three-parameter equation of state for associating fluids. The modified equation of state is shown to be as accurate as the original SAFT model for correlating pure-component vapor-liquid equilibrium data with fewer adjustable parameters. Introduction Molecular association has an important effect on the phase behavior of pure fluids and mixtures. As a result, much effort has been expended on modeling the effects of association on fluid phase equilibria. These efforts can be separated into two classes: chemical theories (Anderko, 1989a,b; Ikonomou and Donohue, 1986) that postulate the formation of distinct molecular aggregates in solution and physical theories that model the association as being the result of a strong, specific molecular interaction. Both methods suffer from the use of a large number of ill-defined adjustable parameters in order to obtain good fits of experimental phase behavior data. If one were able to determine these association parameters a priori, such methods would be more useful, and possibly even completely predictive. This paper discusses the use of association parameters computed from ab initio molecular orbital calculations in one physical equation of state, the statistical associating fluid theory (SAFT). Here we use the results of our molecular orbital calculations in fitting pure component vapor-liquid equilibrium (VLE) data for water, methanol, and the family of primary carboxylic acids. We show that using the information we have computed from molecular orbital calculations allows a reduction in the number of adjustable parameters in the SAFT equation of state with little or no loss in accuracy. Indeed by using such information, the SAFT model becomes a three-parameter equation of state for associating fluids. Description of the Equation of State The equation of state we have chosen to use in this study is statistical associating fluid theory (SAFT) developed by Chapman and co-workers (Chapman et al., 1989, 1990) from the first-order thermodynamic perturbation theory (TPT1) of Wertheim (1984a,b, 1986a,b). SAFT has been shown to yield accurate correlations for the properties of associating and nonassociating fluids
over a wide range of molecular sizes. (Huang and Radosz, 1990; Economou and Donohue, 1991). In this model pure fluids are treated as chains of equal-sized spherical segments interacting with a square-well potential. Two types of bonds between the segments are superimposed, covalent-like bonds to form chain molecules and association bonds. The residual Helmholtz free energy per mole of molecules in the SAFT model is given by
ares ) aseg + achain + aassoc
The segment Helmholtz free energy per mole of molecules includes hard sphere and dispersion energies hs disp aseg ) maseg 0 ) m(a0 + a0 )
S0888-5885(96)00725-7 CCC: $14.00
(2)
In eq 2, m is the number of segments per molecule, the first pure component parameter, and aseg is the seg0 mental Helmholtz energy per mole of segments. The hard sphere energy per segment is obtained from Carnahan and Starling (1969):
ahs 4η - 3η2 0 ) RT (1 - η)2
(3)
The reduced density (η) is given by
η ) τFmv0
(4)
where τ is a constant equal to 0.74048, F is the superficial molar density of molecules, and v0 is the segment molar volume. The temperature dependence of the segment molar volume in eq 4 is taken from the work of Chen and Kreglewski (1977) on the basis of the Barker-Henderson (1967) result for the square-well potential,
[
( )]
v0 ) v∞ 1 - 0.12 exp * Author to whom correspondence should be addressed. E-mail:
[email protected]. † Contributed from Computational Chemistry and Its Industrial Applications, American Institute of Chemical Engineers Annual Meeting, November 14, 1996.
(1)
-3u0 kT
3
(5)
The temperature-independent volume of an isolated segment, v∞, is the second pure-component parameter. Equation 5 also contains the third pure-component © 1997 American Chemical Society
4042 Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 Table 1. Bonding Schemes Used in This Work
parameter, the temperature-independent square-well depth of nonspecific segment-segment interactions (u0/ k). The dispersion energy per segment is derived from a power series initially fitted by Alder et al. (1972) to molecular dynamics data for the square-well fluid:
adisp 0 RT
[ ][ ] u
)
i
∑i ∑j Dij kT
η
[
τ
(6)
]
(7)
where e/k is a constant that has been correlated with Pitzer’s acentric factor and the critical temperature (Chen and Kreglewski, 1977; Kreglewski, 1984) for various molecules. However, in this work the energy parameter is for segments, not molecules, and a value of 10 for e/k was assumed for all the compounds studied (Huang and Radosz, 1990). The Helmholtz free energy per mole of molecules due to chain formation is (Chapman et al., 1990):
[ ]
1 1- η 2 achain ) (1 - m) ln RT (1 - η)3
(8)
In order to describe association, molecules are assigned bonding sites, and the interactions between these sites are modeled using the square-well potential. This introduces two additional parameters for pure associating fluids. The well-depth of the interaction between sites A and B is described by the parameter AB/k, and the well-width is characterized by the dimensionless parameter κAB. In the SAFT model, these association parameters are for the interaction between bonding sites on an isolated pair of molecules. The contribution to the residual Helmholtz free energy due to association is calculated from
aassoc ) RT
∑ A
[
ln XA -
]
XA
XA ) [1 +
FXB∆AB]-1 ∑ B
(10)
j
The Dij are universal constants, and the values used are those fitted by Chen and Kreglewski (1977) to PVT, internal energy, and second virial coefficient data for argon. The temperature-dependent segment energy, u/k from eq 6, is given by
e u u0 ) 1+ k k kT
molecule, XA is the mole fraction of the original molecules that are not bonded at site A, and the summation is over all the bonding sites of a molecule. The fractions of unbonded sites are determined from
1 + M 2 2
(9)
where M is the number of bonding sites assigned to a
The quantity ∆AB is a measure of the association strength and may be approximated (Chapman et al., 1988) as
[ ( ) ]
∆AB ) g(d)seg exp
AB - 1 (x2v∞κAB) kT
(11)
Here g(d)seg is the segmental radial distribution function, taken to be the hard sphere radial distribution function, evaluated at contact. Before using eq 9 to calculate the energy due to association, one must specify the allowed site-site interactions and the values of ∆ needed to determine the unbonded fractions XA. The resulting expressions can be simplified by assuming the equivalence of specific bonds. Table 1 presents the association schemes used in this work, including bonding approximations and expressions for XA. The association scheme nomenclature used is that of Huang and Radosz (1990). Description of the Molecular Orbital Calculations We have performed ab initio molecular orbital calculations to determine the enthalpy, internal energy, entropy, Gibbs free energy, and heat capacity changes for a number of hydrogen-bonding dimerization “reactions”. These calculations are for the association of isolated pairs of molecules and have been performed using two different calculational methods: the computationally inexpensive Hartree-Fock (H-F) method with the small 6-31g(d,p) basis set and the more rigorous density functional theory (DFT) method with the larger 6-31++g(2d,p) basis set. For the DFT calculations we used Becke’s three-parameter functional with the gradient correction to the correlation of Lee, Yang, and Parr (Becke, 1992a,b,1993; Lee et al., 1988). Details of these calculations and a comparison of our results with experimentally determined estimates of the thermodynamic property changes on association were presented earlier (Wolbach and Sandler, 1997a,b). Our results for the compounds considered here are presented in Table 2. The dimers of water and methanol are linear, as these were determined to be thermodynami-
Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 4043 Table 2. Thermodynamic Parameters of Dimerization Determined by Molecular Orbital Calculations H-F
DFT
compound
∆H (kcal/mol)
∆S (cal/mol‚K)
∆CP0 (298 K) (cal/mol‚K)
Keq (298 K)
∆H (kcal/mol)
∆S (cal/mol‚K)
∆CP0 (298 K) (cal/mol‚K)
Keq (298 K)
water methanol formic acid acetic acid propionic acid
-3.929 -3.877 -13.557 -13.924 -13.720
-18.324 -19.937 -34.515 -34.052 -34.311
2.48 3.07 2.22 2.26 2.27
7.480 × 10-2 3.05 × 10-2 248 582 362
-2.993 -3.130 -13.497 -13.958 -13.742
-21.087 -22.379 -37.290 -36.097 -37.205
1.97 2.75 1.39 1.46 1.50
3.847 × 10-3 2.53 × 10-3 55.5 220 87.5
ally favored over cyclic dimers at ambient conditions, while the three acid dimers are cyclic. Procedure There is no direct link between our calculated thermodynamic changes and the association parameters of SAFT, as that equation of state does not use the enthalpy and entropy changes of association as parameters. Thus we had to develop a separate procedure to relate the association parameters of SAFT to the results of our quantum-mechanical calculations. From our calculations we obtain values for the equilibrium constant (Keq) from which we can determine the composition of the vapor phase, and in particular the mole fraction of monomers. In deriving the expressions for the Helmholtz free energy of association in SAFT, it was assumed that the activity of a particular bonding site was unaffected by the bonding states of the other sites on the molecule. This allows us to interpret the values of XA from the association term of SAFT as the probabilities that sites are not bonded. By applying simple probabilistic arguments, we use the values of XA to determine the fraction of the original molecules that have formed a particular number of bonds. This in turn allows us to find the composition of the vapor phase and also the mole fraction of monomers. We then equate the mole fraction of monomers calculated from our quantum-mechanical results to that calculated from the values of the XA in the association term of SAFT. Since the XA are functions of the SAFT parameters AB/k and κAB through ∆AB, we can then derive relationships between AB/k and κAB and our calculated values ∆H, ∆S, and ∆CP. The relationships thus derived will be specific for each association scheme. We begin with the 1A association scheme, and the definition of the equilibrium constant
(
Keq ) exp -
xjφjP/Pref ∆G ) ) RT (x1φ1P/Pref)(xj-1φj-1P/Pref) xjφj 1 atm (12) x1φ1xj-1φj-1 P
)
In eq 12, the reference pressure is 1 atm, and the system pressure (P) will also be expressed in atmospheres. For compounds described by the 1A association scheme, only monomers and dimers form so that eq 12 becomes
Keq )
x2φ2 1 atm (1 - x1)φ2 1 atm ) P x12φ12 P x12φ12
x1 )
-1 + x1 + 4(PKeq/1 atm) (2PKeq/1 atm)
n1 NT
(14)
where n1 is the number of monomers and NT is the total number of molecules after association. From SAFT, we obtain the fraction of the original molecules that have not formed bonds as
f0 ) XA )
-1 + x1 + 4F∆AB n1 ) N0 2F∆AB
(15)
where N0 is the total number of molecules before association. The quantities x1 and f0 are not equivalent. However, for a monomer-dimer equilibrium N0 ) n1 + 2n2 and x1 is related to f0 by
( ) ( ) ( )
n1 NT n1 x1 x1 ) ) ) f0 ) n1 + 2n2 n1 n2 x1 + 2x2 2 - x1 +2 NT NT (16) We next assume that at low pressures the vapor phase can be treated as an ideal-gas mixture, so that after one has accounted for the reduction in the apparent number of molecules due to association
F)
P N0 RT NT
(17)
It is important to note that the ideal-gas mixture assumption will only be applied at very low pressures in order to use our quantum-mechanical results to fix the values of the SAFT association parameters and that we are not making this assumption over the entire vapor-liquid coexistence range. The SAFT equation will account for any vapor phase nonidealities at higher pressures. In eq 17, F is the superficial molar density (N0/V), and N0/NT is
n1 N0 N0 1 ) ) x1 ) 2 - x1 NT NT n1 f0
(18)
Inserting eqs 16-18 into eq 15 leads to the following expression for the mole fraction of monomers on the basis of the value of ∆AB from SAFT:
(13)
As we will be applying our procedure to evaluate the SAFT parameters to low pressure data (and then making predictions at high pressures), we assume in this parameter determination step that the fugacity coefficients of the monomer and dimer are each equal to unity so that the mole fraction of monomers is
)
-1 + ) xSAFT 1
x ( ) ( ) 1+2
P∆AB RT
P∆AB RT
(19)
A comparison of eqs 14 and 19 reveals that
[ ( ) ]
2RTKeq AB ) ∆AB ) ghs(d)x2v∞κAB exp -1 1 atm kT
(20)
4044 Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997
where for low-temperature vapor phase densities, it may be assumed that ghs(d) is equal to 1. The temperature dependence of the equilibrium constant is
( )
Keq ) exp -
(
{
( )})
∆H ∆S ∆CP T0 - T T exp + + + ln RT R R T T0
Keq )
xj
1 atm P
xj-1x1
(22)
The mole fraction of an n-mer (after association) may be expressed as
(
)
x1PKeq xj-1x1PKeq ) xj-2 1 atm 1 atm
2
) ... )
(
x1
)
x1PKeq 1 atm
j-1
)
nj (23) NT
Next, the mass (mole) balance leads to the following expression for x1: ∞
1)
∑ j)1
∞
xj )
∑ j)1
x1
( ) x1PKeq
j-1
1 atm
x1 1 - (x1PKeq/1 atm)
NT
n1 ∞
)
)
(
f x1 ) 1 +
PKeq
)
1 atm
n1
-1
) NT
(24)
(1 + 2F∆AB) - x1 + 4F∆AB 2(F∆AB)2
n1 N0 (25)
)
The number of molecules present before association can be determined from ∞
N0 )
ini ∑ i)1
(26)
leading to the following relationship between f0 and x1
j
j)1
nj
NT
∞
x1
( ) ]
) ∞
jxj ∑jx1 ∑ j)1 j)1 x1
x1PKeq
) j-1
1 atm
) x12 (27)
x1
(1 - {x1PKeq/1 atm})2
We again treat the vapor phase as an ideal-gas mixture and express the density as in eq 17, with the ratio N0/ NT for the two-site association model being
N0 n1 N0 1 1 ) ) x1 ) NT NT n1 f0 x1
(28)
Inserting the results of eqs 17, 27, and 28 into eq 25, we derive the mole fraction of monomers from SAFT:
(
)
∆AB RT
) 1+P xSAFT 1
-1
(29)
Comparison of eqs 29 and 24 reveals that for the 2B association model
∆AB )
RTKeq 1 atm
(30)
Once again, no unique relationship between the association parameters of SAFT and the thermodynamic parameters calculated from quantum mechanics may be derived from a calculation at a single temperature. The expressions in eqs 20 and 30 are identical to those obtained by Economou and Donahue (1991) for one- and two-site association models, given the assumption that the vapor phase is an ideal-gas mixture. Economou and Donahue (1992) extended their analysis to three- and four-site association models by assuming that eq 30 would also be correct in those cases. This is not correct as eqs 20 and 30 do not show that ∆AB ) (RTKeq/1 atm), but rather only that ∆AB is proportional to (RTKeq/1 atm). In their study this distinction was not critical, as they treated the equilibrium constant as an adjustable parameter that would absorb any constant of proportionality. However, here the equilibrium constant is fixed a priori, so we cannot make this assumption. Instead, we start with eq 12, and assuming that the fugacity coefficients of all species are equal to one, write
x2 PKeq ) 2 1 atm x
For the 2B association model, the mole fraction of unbonded molecules from SAFT is
f0 ) XAXB ) (XA)2 )
∞
x1 )
[
(21)
with T0 being the temperature at which the thermodynamic parameters of association have been determined using quantum mechanics. It is not possible to use eqs 20 and 21 to derive a unique relationship between the SAFT parameters AB/k and κAB and the quantum-mechanically determined ∆H, ∆S, and ∆CP from a calculation at a single temperature. Instead, one can either use eqs 20 and 21 at two different temperatures to eliminate AB/k and κAB in terms of the quantum-mechanically determined thermodynamic properties or use these equations at a single temperature to interrelate AB/k and κAB, thus eliminating one of those parameters. The second procedure is equivalent to allowing ∆H and ∆S to vary, but constraining them to reproduce the calculated equilibrium constant at a single temperature. For the 2B association scheme, we again begin with eq 12 and assume that the vapor phase fugacity coefficients of all species are equal to one
xj )
f0 )
jnj ∑ j)1
∆G ) RT
() ∑( ) n1
(31)
1
We then proceed to calculate the mole fraction of monomers and dimers in the vapor phase from the values of XA in the association term of SAFT for the 3B and 4C association schemes. This procedure is somewhat involved, and the details of the calculations are presented as an appendix. The final expression for the 3B model is
∆AB )
RTKeq [1 + XA + 2(XA)2](2XA - 1)XA 1 atm 2(3XA - 1)2 (32)
This equation is recursive, as XA is a function of ∆AB
Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 4045 Table 3. Comparison of Calculated Volumes with Experimental Solid Molar Volumesa experimental solid H-F DFT molar volumea molar volume molar volume (cm3/mol) (cm3/mol) (cm3/mol)
compound water methanol formic acid acetic acid propionic acid hydrogen ethane
18.857 32.696 30.281 59.963 62.478 43.327
15.907 30.683 31.254 45.373 59.161 10.832 39.456
17.744 32.345 33.673 47.806 61.903 10.688 40.209
a
From Data Compilation Tables of Properties of Pure Compounds; Danner, R. P., Daubert, T. E., Eds.; Design Institute for Physical Property Data, American Institute of Chemical Engineers: New York, 1985.
the Monte-Carlo integration. In all cases, the standard deviation of the volume distribution was less than 1% of the mean, indicating that our calculations were relatively precise. The optimum values of the segmental energy (u0/k) and either v∞ or m are determined by fitting VLE data. Values of the quantum-mechanically determined Keq at the two lowest temperatures for which data are available were used to determine both SAFT association parameters a priori. A constrained (m g 1, u0/k g 0) two-variable simplex search (routine DBCPOL, IMSL MATH/LIBRARY, version 1.0, 1987) was then used to minimize the objective function:
y)
(see Table 1). In the limit of little association (XA f 1) eq 32 simplifies to
∆AB)
1 RTKeq 2 1 atm
(33)
For the 4C model, the final expression is
∆AB)
RTKeq [1 + (XA)2](XA + 1) 1 atm 16
(34)
which in the limit of XA f 1 simplifies to
∆AB)
1 RTKeq 4 1 atm
(35)
In this work we have used eqs 33 and 35 when performing fits using the 3B or 4C association models. Results In addition to using the results of our quantummechanical calculations to eliminate at least one of the association parameters in the SAFT equation without a loss of accuracy, we show below we are also able to give physical meaning to the segment volume (v∞). The package we used for the quantum mechanics calculations (Gaussian Inc., 1993, 1995) can be used to calculate the molar volume of a molecule using Monte-Carlo integration to determine the extent of its electron density. These calculations have been performed using both the H-F and DFT methods, and our results are presented in Table 3. We interpret our calculated volumes to be the volume per molecule in the close-packed solid, as the results of our calculations correspond well to the solid molar volumes of the compounds studied, with the exception of acetic acid (whose reported solid volume we suspect may be in error). Note, however, that the volume parameter in SAFT (v∞) is the volume of a single, isolated segment. For chain molecules of equal-sized segments, the relationship between the volume of a mole of isolated segments and the close-packed (solid) molar volume is
molar volume ) x2mv∞
(36)
In this manner, we are able to constrain the product of the parameters m and v∞ to be equal to our calculated molar volume, thereby eliminating one of them. Due to random errors inherent in Monte-Carlo integration, we took our molar volumes for each molecule to be the average of 25 calculations, and the option volume ) tight was used to select a large number of points for
|Pcalc - Pexp i i |
∑i
|F1,calc - F1,exp | i i
+ Pexp i
∑i
F1,exp i
(37)
However, if the optimized parameters resulted in an average error in vapor pressure of greater than 5%, we instead used the value of Keq at the lowest temperature for which data were available to form a constraint equation on the values of AB/k and κAB, thereby eliminating only one of the association parameters. A constrained (m g 1, u0/k g 0, AB/k g 0) three-variable simplex search was then performed, again using eq 37 as the objective function. The initial guesses for this second optimization were the optimum parameters determined by the first optimization. The data sets used in performing the optimizations were from the same source and over the same temperature ranges as used by Huang and Radosz (1991). Water The first molecule we considered was water. If one were to assign a bonding site to each of the lone pairs and hydrogen-bonding protons present, water would be described by the 4C association model. However, recent experimental evidence (Wei et al., 1991) suggests that in clusters of liquid water only three sites per molecule are bonded, indicating that water should be treated with the 3B association model. We have examined both of these association models in fitting VLE data and have used the results of both our H-F and DFT calculations in performing these fits. Initially we used values of the quantum-mechanically determined Keq at two different temperatures to determine both of the SAFT association parameters a priori. However, regardless of the association model and quantum-mechanical results used, we obtained average errors in the calculated vapor pressure that were too large. Instead, we then used a single value of Keq to establish a constraint equation on the values of AB/k and κAB. Consequently, all of our fits to VLE data involved three adjustable parameters (u0/k, m or v∞, and AB/k or κAB). For comparison, Huang and Radosz (1991) also modeled water with the 3B association model and assigned an arbitrary, fixed segment volume of 10 cm3/ mol. This resulted in a fit involving four adjustable parameters. The best-fit parameters and the average errors in the calculated vapor pressure and liquid density for both of our methods are presented in Table 4. This table also includes the average errors and parameters from Huang and Radosz. For both the 3B and 4C association models, our fits using the results of the H-F calculations are better than those of the DFT calculations. Reduction in the average errors may be obtained using the DFT results if one allows the number of segments in the molecule to be
4046 Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 Table 4. Results for Water and Methanol method
association model
v∞ (cm3/mol)
m
u0/k (K)
AB/k (K)
κAB × 102
% error Pvap
% error Fliq
watera
original SAFT H-F H-F DFT DFT
3B 3B 4C 3B 4C
10.00 8.80 8.00 12.55 12.55
1.179 1.278 1.406 1.000 1.000
528.17 385.12 212.91 615.94 546.63
1810 2286 1809 1627 1237
1.593 3.091 9.109 1.098 2.203
1.22 2.57 1.74 2.70 1.94
3.15 3.32 2.89 4.75 3.69
methanolb
original SAFT H-F H-F DFT DFT
2B 2B 3B 2B 3B
12.00 8.10 8.10 9.50 9.00
1.776 2.679 2.679 2.401 2.535
216.13 235.90 222.57 306.74 288.60
2714 2003 1930 1729 1470
4.856 6.783 4.472 1.134 1.641
0.83 3.08 1.23 7.98 6.88
0.88 0.61 0.64 0.79 0.66
a Data from Wylen, G. J.; Sonntag, R. E. Fundamentals of Classical Thermodynamics, S. I. Version, 2nd ed.; John Wiley & Sons: New York, 1978. Data used over the temperature range 283-613 K. b Data from Thermodynamic Tables for non-Hydrocarbons; Thermodynamic Research Center, Texas A&M University: College Station, TX. Data used over the temperature range 293-393 K.
deviating from the measured liquid density as the critical point is approached. Methanol
Figure 1. Errors in the calculated (a) vapor pressure and (b) liquid density for water using the original SAFT model and our modifications to the original SAFT model. The original SAFT model assumes association according to the 3B scheme and uses four parameters. Our modifications assume association according to the 4C scheme and use three parameters.
less than one. However, as this is unrealistic, we do not report such results here. For the results of both the H-F and DFT calculations, the 4C association model yields a better fit of the data. Using the 4C model and either set of quantum-mechanical results, we have been able to achieve fits of comparable accuracy to the original SAFT with a reduction in the number of adjustable parameters. Also, we have been able to give the parameter v∞ better physical significance, instead of considering it to be an adjustable parameter. Figures 1(a) and 1(b) present the results of the best fits using the 4C association model for both of our threeparameter methods, and of the four-parameter original SAFT model. The three methods perform comparably for the calculated vapor pressures over most of the temperature range, with the fit using the results of our H-F calculations producing the smallest errors as the critical point is approached. The three methods perform nearly identically for the calculated liquid density over the entire temperature range, with all of the methods
Methanol contains one hydrogen-bonding proton and two lone electron pairs. If one includes all of these as potential bonding sites, methanol would be described by the 3B association model. Alternatively, one could postulate that it is not possible to bond to both lone electron pairs simultaneously and then treat methanol with the 2B association model. We have tested both models in fitting VLE data and have used the results of both our H-F and DFT calculations. As was the case for water, using quantum-mechanically calculated values of Keq at the two lowest temperatures for which data were available resulted in average errors in the vapor pressure which were too large regardless of the method and model used. By using the value of Keq at the lowest temperature for which data were available, we eliminated one of the association parameters. As a result, our fits to VLE data involved three adjustable parameters. For comparison, in the original SAFT model, methanol was treated with the 2B association model, assigned an arbitrary, fixed segment volume of 12 cm3/mol, and used four adjustable parameters. The average errors in the calculated vapor pressure and liquid density of our fits, along with that from Huang and Radosz, are presented in Table 4. The optimized parameter values for all the methods are also included in Table 4. The fit we obtain using the 3B association model and the results of our H-F calculations is comparable to the fit of the original SAFT model, even with the elimination of one adjustable parameter. As was the case for water, the more intuitive association model produces better results. Since methanol and water are similar compounds, it is reassuring that the optimum segment volumes and segment energies (u0/k) in our best fits for these compounds are similar. We were unable to obtain a fit of the vapor pressure with average errors of less than 5% using the DFT results with either association model. Figures 2(a) and 2(b) show that while we were able to fit the liquid density well, we were unable to correctly predict the curvature of the vapor pressure curve, leading to large errors. The shape of the error curve for the DFT results in Figure 2(a) is similar to the error curve for the DFT results in Figure 1(a) for water. However, in the case of water the maximum error is about 6%, while for methanol it is over 10%, and the bowing is more pronounced. The degree of bowing in the vapor pressure error curve appears to be inversely proportional to the magnitude of the equilibrium constant.
Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 4047
of the -COOH unit was calculated as the average of 1/2 the volume of H2 subtracted from the volume of formic acid and 1/2 the volume of ethane subtracted from the volume of acetic acid. This served to “remove” the volume of the non-COOH atoms from the acids. We obtained a value of v∞ for the acids of 19.8 cm3/mol from the results of our DFT calculations, and 18.2 cm3/mol from the results of our H-F calculations. We have used both our H-F and DFT results in fitting the experimental data, with nearly equal success. The average errors in calculated vapor pressure and liquid density from both methods, along with the errors from Huang and Radosz, are presented in Table 5, while our best-fit parameters for both calculational methods are presented in Table 6. Both methods show a decrease in the value of the well-depth for segmentsegment interactions (u0/k) with increase in molecular size. This may be an asymptotic approach to a limiting value for methyl group interaction. Also, both methods show consistent values of the well-depth of association (OH/k) for all the acids, except for the smallest of the group, formic acid. We next present a more detailed analysis for some specific compounds. 1. Formic Acid. In our fits of formic acid, setting both SAFT association parameters a priori resulted in vapor pressure curves which were too steeply sloped. As a result, our fits involved three adjustable parameters: m, u0/k, and AB/k or κAB. Figure 3(a) displays the deviations from the experimental vapor pressure obtained using our H-F and DFT results, along with the errors obtained in the original SAFT work. All of the methods show good agreement with the experimental data and display similar trends in the error with temperature. However, we have obtained fits of equal quality to those reported by Huang and Radosz while eliminating one, and possibly two, adjustable parameters (it is unclear whether v∞ was treated as an adjustable parameter in their correlations using SAFT). Figure 3(b) compares the deviations from the experimental liquid densities using the method here and that of Huang and Radosz. While all of the methods are in excellent agreement with the experimental data, our correlations using either the H-F or DFT results are more accurate than that of Huang and Radosz and, in fact, are within the experimental accuracy. 2. Acetic Acid. As was the case with formic acid, our fits involved three adjustable parameters. The deviations from the experimental vapor pressure are shown in Figure 4(a). Although much of the error could
Figure 2. Errors in the calculated (a) vapor pressure and (b) liquid density for methanol using the original SAFT model and our modifications to the original SAFT model. The original SAFT model assumes association according to the 2B scheme and uses four parameters. Our modifications assume association according to the 3B scheme and use three parameters.
Acids One of the restrictions of TPT1 (Wertheim, 1984a,b, 1986a,b) is that pairs of molecules may only singly bond, so in order to account for the formation of cyclic dimers in the vapor phase of the carboxylic acids, we must treat acids using the 1A association model. This assumption restricts both the vapor and liquid phases to be composed of only monomers and dimers and is the same choice made by Huang and Radosz in the original SAFT model. Since we were performing fits for an entire family of primary carboxylic acids, we felt it would be advantageous to use a single value of the segment volume (v∞) for the entire family of acids. The common segment volume was chosen to be that of the bonding functionality in carboxylic acids, the -COOH unit. The volume
Table 5. Average Errors of the Original SAFT and Our Modifications for Acids acid
T range (K)
formic
293-393
acetic
283-573
propionic
313-463
1-butanoic
333-493
1-pentanoic
353-483
a
TX.
model
% error Pvap
% error Fliq
original SAFT H-F DFT original SAFT H-F DFT
0.62 0.76 0.62 1.90 0.86 0.94
original SAFT H-F DFT original SAFT H-F DFT original SAFT H-F DFT
0.35 0.87 0.75 0.20 0.76 0.69 0.32 0.85 0.68
acid
T range (K)
0.49 0.27 0.10 0.65 1.58 0.63
1-hexanoic
372-504
1-heptanoic
385-522
0.09 0.67 0.87 0.85 1.19 0.72 0.98 1.25 0.40
1-octanoic
399-540
1-nonanoic
411-557
1-decanoic
423-572
model
% error Pvap
original SAFT H-F DFT original SAFT H-F DFT
0.42 1.09 0.73 0.47 1.02 0.52
original SAFT H-F DFT original SAFT H-F DFT original SAFT H-F DFT
0.56 1.18 0.58 0.46 1.19 0.58 0.55 1.06 0.75
Data from Thermodynamic Tables for non-Hydrocarbons; Thermodynamic Research Center, Texas A&M University: College Station,
4048 Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 Table 6. Best Fit Parameters for Our Modifications to SAFT for Acids H-F acid
v∞ (cm3/mol)
formic acetic propionic 1-butanoic 1-pentanoic 1-hexanoic 1-heptanoic 1-octanoic 1-nonanoic 1-decanoic
18.20 18.20 18.20 18.20 18.20 18.20 18.20 18.20 18.20 18.20
DFT
m
u0/k (K)
AB/k (K)
1.176 1.684 2.192 2.700 3.208 3.716 4.224 4.732 5.240 5.748
369.86 305.97 282.51 272.87 268.97 268.27 264.97 262.84 261.03 259.15
4160 5751 6123 6217 6165 6239 6168 6108 6064 6050
κΑΒ × 103
v∞ (cm3/mol)
m
u0/k (K)
AB/k (K)
κΑΒ × 103
582.3 3.465 0.855 0.580 0.631 0.485 0.565 0.635 0.684 0.685
19.75 19.75 19.75 19.75 19.75 19.75 19.75 19.75 19.75 19.75
1.099 1.582 2.065 2.548 3.031 3.514 3.997 4.480 4.963 5.446
394.33 321.28 306.75 300.23 296.65 296.42 291.44 287.84 284.81 282.13
3771 5490 5395 5604 5593 5853 5846 5847 5984 6313
464.5 4.157 1.715 0.728 0.649 0.281 0.267 0.252 0.172 0.077
Figure 3. Errors in the calculated (a) vapor pressure and (b) liquid density for formic acid using the original SAFT model and our modifications to the original SAFT model. All models assume association according to the 1A scheme. The original SAFT model uses four parameters, while our modifications use three parameters.
be attributed to scattering in the experimental data, our correlation is better over most of the temperature range than that of Huang and Radosz with one less parameter, regardless of whether we used the H-F or DFT results in the SAFT model. The errors in the calculated liquid density for all of the methods are displayed in Figure 4(b). The original SAFT and our modification perform almost identically, with the lack of a distinctive trend in the errors most likely being the result of scattering in the experimental data. All of the correlation results deviate markedly from the experimental data as the critical point is approached, indicating that the SAFT model with parameters fit to low-pressure data does not correctly predict the critical point. However, as with the vapor pressure results, no significant loss of accuracy over the entire temperature range has occurred as a result of one adjustable parameter being specified a priori. 3. Propionic Acid and Larger. Propionic acid is the largest acid for which we have performed molecular orbital calculations due to limitations on CPU time and physical disk storage (for a discussion of the scaling of
Figure 4. Errors in the calculated (a) vapor pressure and (b) liquid density for acetic acid using the original SAFT model and our modifications to the original SAFT model. All models assume association according to the 1A scheme. The original SAFT model uses four parameters, while our modifications use three parameters.
CPU and disk storage with increasing molecular size, see Wolbach and Sandler, 1997a). In order to apply our procedure to acids larger than propionic, we assumed that the thermodynamic parameters of association for the larger acids were equal to those for propionic acid. In addition, we took the number of segments (m) to be linear with carbon number, resulting in fits which involved only two adjustable parameters (u0/k, and κAB or AB/k). This represents a reduction of at least two and perhaps three parameters from the original SAFT model. The results in Table 5 show that it is possible to fit VLE data for the highly-associating acids to high accuracy with a reduction in the number of adjustable parameters regardless of whether we use the results of our H-F or DFT calculations in performing the fits. While the DFT results are slightly more accurate, the H-F procedure may be preferred due to its lesser computational requirements. Discussion Although we have been able to achieve results comparable to the original SAFT using fewer adjustable
Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 4049
parameters in all cases but one, we were unable to eliminate both of the association parameters in SAFT a priori, which was an initial goal of this work. One reason for this may be the inherent inaccuracy of the SAFT model due to its approximate nature. A second reason is that what we have calculated from quantum mechanics and what is needed in the equation of state are not exactly the same. We have calculated thermodynamic parameters strictly for monomer-dimer transitions. For molecules that possess more than one bonding site, SAFT assumes that all the possible linear associations have the same association parameter. As a result, the best-fit association parameters in SAFT actually represent some average of the parameters for all of the possible linear association schemes. For molecules that have one bonding site (i.e., acids), only dimers are allowed to form. However, there is experimental evidence that the acid crystalline phase is composed of linear chains (Jo¨nsson, 1971; Bertagnolli and Hertz, 1978), and it is believed that the liquid phase of acids also contains linear chains, which cannot be accounted for by the monomer-dimer model. The bestfit association parameters in SAFT may thus have to deviate from the strictly monomer-dimer values we have calculated in order to account for these other species. Conclusions This work shows the applicability of quantummechanically derived thermodynamic parameters of association in reducing the number of parameters necessary to describe the phase behavior of hydrogenbonding fluids. By using the results of our quantummechanical calculations, we have been able to fit purecomponent VLE data for a number of associating fluids using the SAFT model with a reduction in the number of freely-adjustable parameters and little or no loss of accuracy. The results of our H-F calculations are applicable to all the compounds considered in this study, while the more computationally-intensive DFT calculations are more accurate for the family of carboxylic acids. While it may seem troubling that the more rigorous DFT calculations yield inferior results to the H-F calculations in some cases, one should recall that this study considered a single equation of state and used only the parameters for monomer-dimer transitions. Study of other equations of state, or use of higher n-mer parameters, may yield different results. While we have not been able to use the results of our quantum-mechanical calculations to eliminate both of the association parameters, we believe this difficulty is related to the assumptions in the SAFT association model. A barrier to using more complex (and correct) association schemes is the need to introduce additional parameters for each type of association proposed. The results of quantum-mechanical calculations may allow us to be able to determine these parameters a priori, and the work here represents a first step in this process. Acknowledgment The authors would like to thank the National Science Foundation (Grant No. CTS-9123434) and the Department of Energy (Grant No. DE-FG02-93ER13436) for financial support of this research. J.W. would like to thank the National Science Foundation for a PreDoctoral Fellowship.
Nomenclature xn ) true mole fraction of an n-mer. Number of n-mers divided by the total number of molecules present after association fn ) mole fraction of original molecules which have formed n bonds. Number of molecules involved in n-mers divided by the number of molecules present before association ni ) number of i-mers present after association occurs N0 ) total number of molecules before association occurs; the true number of molecules NT ) total number of molecules present after association. Treating pure component as a mixture of oligomers. The apparent number of molecules.
Appendix This appendix describes the calculation of the vapor phase composition for the 3B and 4C association models. For the 3B model (commonly used to describe methanol and water), the fraction of the original molecules which are unbonded is given by
f0 ) XAXBXC ) (XA)2(2XA - 1)
(A1)
where the expression for XA is given in Table 1. The mole fraction of monomers is the number of unbonded molecules divided by the total number of molecules present after association which is the number of monomers plus the total number of chains formed. Due to the assumptions of TPT1, only chains and branched chains may be formed. Every unbranched chain has two chain ends, and every branched chain has an additional chain end for each branch. Thus, we can calculate the number of chains formed by knowing the total number of chain ends and branch points. A segment is a chain end if it has formed only one bond and is a branch point if it has formed three bonds. For the 3B model, the fraction of the original molecules that are chain ends is
f1 ) (1 - XA)XBXC + XA(1 - XB)XC + XAXB(1 - XC) ) 2XA(1 - XA)(3XA - 1)
(A2)
while the fraction that are branch points is
f3 ) (1 - XA)(1 - XB)(1 - XC) ) 2(1 - XA)3 (A3) The total number of chains then is
Nchain )
N0f0 - N0f3 ) N0[1 - (XA)2](2XA - 1) (A4) 2
and the mole fraction of monomers is given by
xSAFT ) 1
N0f0 ) N0f0 + Nchain (XA)2(2XA - 1)
(XA)2(2XA - 1) + [1 - (XA)2](2XA - 1)
) (XA)2 (A5)
Calculation of the mole fraction of dimers is slightly more involved. To form a dimer, a molecule that has formed one bond must attach to a second molecule that has only one bond. We begin by recognizing that the fraction of the original molecules that have formed one bond is f1, as given by eq A2, and we need to determine
4050 Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997
to what these molecules have bonded. In doing so, we want to consider only those molecules that have formed bonds, so we remove the unbonded molecules from the mixture and renormalize the distribution of molecules on a monomer-free basis. The fraction of molecules that have formed one bond on a monomer-free basis is
f1′ )
f1 2XA(1 - XA)(3XA - 1) ) ) 1 - f0 1 - (XA)2(2XA - 1) A
X0 ) XAXBXCXD ) (XA)4
1 - XA + 2(XA)2
(A6)
We now select a molecule in the monomer-free mixture that has formed one bond. In a random process, it will be bonded to a second molecule that has formed one bond with a probability equal to the fraction of molecules in the monomer-free mixture that have formed only one bond. Therefore, the result in eq A6 is the probability that a given molecule that has formed one bond is part of a dimer and the fraction of the original molecules that are in dimer form is the product of f1 and f1′. Therefore
4(XA)2(1 - XA)(3XA - 1)2 1 + XA + 2(XA)2
(A7)
The quantity fdimer counts each molecule in the dimer, so the actual number of dimers formed is fdimer/2. The mole fraction of dimers is then the number of dimers divided by the total number of molecules:
(
A RT (1 - X ) P XA
)
N0fdimer 2(XA)2(1 - XA)(3XA - 1)2 2 ) x2SAFT ) N0f0 + Nchain [1 + XA + 2(XA)2](2XA - 1) (A8) Inserting the results of eqs A5 and A8 into eq 31 allows us to determine the relationship between XA and the quantum-mechanically derived equilibrium constant:
(A12)
Equation 33 is the result of dividing eq A12 by eq A9. In the 4C association model (another model that has been used for water), the mole fraction of unbonded molecules is given by
A
2X (3X - 1)
fdimer ) f1f1′ )
∆AB )
(A13)
where XA is given for the 4C model in Table 1. As was the case for the 3B model, we will take advantage of the assumption that only chains and branched chains are allowed to form. Unbranched chains have two chain ends, chains with branch points contain an additional chain end for each branch, and chains that have doublebranch points, or molecules which have formed 4 bonds, will have two additional chain ends for each doublebranch point. The fractions of the original molecules that are chain ends is given by
f1 ) 4(1 - XA)(XA)3
(A14)
while the fraction of branch points is
f3 ) 4(1 - XA)3(XA)
(A15)
and the fraction of double-branch points is
f4 ) (1 - XA)4
(A16)
The total number of chains may then be calculated from
Nchain )
N0f1 - N0f3 - 2N0f4 ) 2 A A 3 N0(1 - X )[(X ) + (XA)2 + XA - 1] (A17)
and the mole fraction of monomers is given by
SAFT
PKeq x2 2(1 - XA)(3XA - 1)2 ) SAFT 2 ) 1 atm (x ) [1 + XA + 2(XA)2](2XA - 1)(XA)2 1
(A9) Our solution procedure is as follows: (a) Using the experimental vapor pressure and the quantum-mechanically determined equilibrium constant, solve eq A9 to fix the value of XA. (b) Invert the expression for XA for the 3B model given in Table 1 and solve the resulting equation for the quantity F∆AB:
F∆AB )
(1 - XA) XA(2XA - 1)
(A10)
(c) Determine the value of the quantity N0/NT for the 3B model:
N0 x1SAFT ) ) (2XA - 1)-1 NT f0
(A11)
(d) Use the experimental vapor pressure, the experimental temperature, and the value of N0/NT from eq A11 in eq 17 to eliminate the vapor density in eq A10. This allows the determination of the value of the quantity ∆AB at the experimental temperature as follows:
x1SAFT )
N0f0 ) N0f0 + Nchain N0(XA)4
N0(XA)4 + N0(1 - XA)[(XA)3 + (XA)2 + XA - 1] (XA)4 )
2XA - 1
(A18)
Since we will be using this procedure at low temperatures and pressures in the vapor phase, we expect small degrees of association and values of XA ≈ 1, thereby avoiding the neighborhood of the singularity in eq A18 at XA ) 0.5. The fraction of the original molecules present as dimers is calculated in the same manner as for the 3B model:
fdimer ) f1f1′ )
(f1)2 16(XA)6(1 - XA)2 ) ) 1 - f0 1 - (XA)4 16(XA)6(1 - XA) (1 + XA)[1 + (XA)2]
The mole fraction of dimers is
(A19)
(
Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 4051
)
N0fdimer 2 x2SAFT ) ) N0f0 + Nchain 8(XA)6(1 - XA) (1 + XA)[1 + (XA)2](2XA - 1)
(A20)
The relationship between XA and the quantum-mechanically derived equilibrium constant is
PKeq x2SAFT 8(2XA - 1)(1 - XA) ) SAFT 2 ) 1 atm (x ) (1 + XA)[1 + (XA)2](XA)2 1
(A21)
We then follow the same procedure as for the 3B model, with the equivalent of eq A10 being
F∆AB )
(1 - XA) 2(XA)2
(A22)
and the analog of eq A11 given by
N0 x1SAFT ) ) (2XA - 1)-1 NT f0
(A23)
to obtain
∆AB )
A A RT (1 - X )(2X - 1) P 2(XA)2
(A24)
Equation 34 is then derived by dividing eq A24 by eq A21. Literature Cited Alder, B. J.; Young, D. A.; Mark, M. A. Studies in Molecular Dynamics. X: Corrections to the Augmented van der Waals Theory for Square-Well Fluid. J. Chem. Phys. 1972, 56, 3013. Anderko, A. A Simple Equation of State Incorporating Association. Fluid Phase Equilib. 1989a, 45, 39. Anderko, A. Extension of the AEOS Model to Systems Containing Any Number of Associating and Inert Components. Fluid Phase Equilib. 1989b, 50, 21. Barker, J. A.; Henderson, D. J. Perturbation Theory and Equation of State for Fluids. II. A Successful Theory of Liquids. J. Chem. Phys. 1967, 47, 4714. Becke, A. D. Density-functional thermochemistry. I. The effect of the exchange-only gradient correction. J. Chem. Phys. 1992a, 96 (3), 2155. Becke, A. D. Density-functional thermochemistry. II. The effect of the Perdew-Wang generalized-gradient correlation correction. J. Chem. Phys. 1992b, 97 (12), 9173. Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98 (7), 5648. Bertagnolli, H.; Hertz, H. G. Preservation and Loss of Structural Features of Solid Acetic Acid during the Melting Process. Phys. Status Solidi A 1978, 49, 463. Carnahan, N. F.; Starling, K. E. Equation of State for Nonattracting Rigid Spheres. J. Chem. Phys. 1969, 51, 635. Chapman, W. G.; Jackson, G.; Gubbins, K. E. Phase Equilibria of Associating Fluids: Chain Molecules with Multiple Bonding Sites. Mol. Phys. 1988, 65, 1057. 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.
Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. New Reference Equation of State for Associating Liquids. Ind. Eng. Chem. Res. 1990, 29, 1709. Chen, S. S.; Kreglewski, A. Applications of the Augmented van der Waals Theory of Fluids. I: Pure Fluids. Ber. Bunsen-Ges Phys. Chem. 1977, 81, 1049. Economou, I. G.; Donohue, M. D. Chemical, Quasi-Chemical and Perturbation Theories for Associating Fluids. AIChE J. 1991, 37, 1875. Economou, I. G.; Donohue, M. D. Equation of State with Multiple Associating Sites for Water and Water-Hydrocarbon Mixtures. Ind. Eng. Chem. Res. 1992, 31, 2388. Gaussian 92/DFT, Revision G.2: M. J. Frisch, G. W. Trucks, H. B. Schlegel, P. M. W. Gill, B. G. Johnson, M. W. Wong, J. B. Foresman, M. A. Robb, M. Head-Gordon, E. S. Replogle, R. Gomperts, J. L. Andres, K. Raghavachari, J. S. Binkley, C. Gonzalez, R. L. Martin, D. J. Fox, D. J. Defrees, J. Baker, J. J. P. Stewart, and J. A. Pople, Gaussian, Inc., Pittsburgh PA, 1993. Gaussian 94, Revision B.2: M. J. Frisch, G. W. Trucks, H. B. Schlegel, P. M. W. Gill, B. G. Johnson, M. A. Robb, J. R. Cheeseman, T. Keith, G. A. Petersson, J. A. Montgomery, K. Raghavachari, M. A. Al-Laham, V. G. Zakrzewski, J. V. Ortiz, J. B. Foresman, J. Cioslowski, B. B. Stefanov, A. Nanayakkara, M. Challacombe, C. Y. Peng, P. Y. Ayala, W. Chen, M. W. Wong, J. L. Andres, E. S. Replogle, R. Gomperts, R. L. Martin, D. J. Fox, J. S. Binkley, D. J. Defrees, J. Baker, J. P. Stewart , M. Head-Gordon, C. Gonzalez, and J. A. Pople, Gaussian, Inc., Pittsburgh PA, 1995. Huang, S. H.; Radosz, M. Equation of State for Small, Large, Polydisperse, and Associating Molecules. Ind. Eng. Chem. Res. 1990, 29, 2284. Ikonomou, G. D.; Donohue, M. D. Thermodynamics of HydrogenBonded Molecules: The Associated Perturbed Anisotropic Chain Theory. AIChE J. 1986, 32 (10), 1716. Jo¨nsson, P. G. Hydrogen Bond Studies. XLIV. Neutron Diffraction Study of Acetic Acid. Acta Crystallogr. 1971, B27, 893. Kreglewski, A. Equilibrium Properties of Fluids and Fluid Mixtures; The Texas Engineering Experiment Station (TEES) Monograph Series; Texas A&M University Press: College Station, 1984. Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of electron density. Phys. Rev. B 1988, 37, 785. Wei, S.; Shi., Z.; Castleman, A. W., Jr. Mixed Cluster Ions as a Structure Probe: Experimental Evidence for Clathrate Structure of (H2O)20H+ and (H2O)21H+. J. Chem. Phys. 1991, 94, 3268. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. I. Statistical Thermodynamics. J. Stat. Phys. 1984a, 35, 19. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. II. Thermodynamic Perturbation Theory and Integral Equations. J. Stat. Phys. 1984b, 35, 35. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. III. Multiple Attraction Sites. J. Stat. Phys. 1986a, 42, 459. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. IV. Equilibrium Polymerization. J. Stat. Phys. 1986b, 42, 477. Wolbach, J. P.; Sandler, S. I. Thermodynamics of Hydrogen Bonding from Molecular Orbital Theory. 1. Water. AIChE J. 1997a, 43 (6), 1589. Wolbach, J. P.; Sandler, S. I. Thermodynamics of Hydrogen Bonding from Molecular Orbital Theory. 2. Organics. AIChE. J. 1997b, 43 (6), 1597.
Received for review November 8, 1996 Revised manuscript received February 19, 1997 Accepted February 19, 1997X IE9607255
X Abstract published in Advance ACS Abstracts, September 1, 1997.
