5115
J . Phys. Chem. 1990, 94, 5 1 15-5120
A Modified Thermodynamic Perturbation Theory Equation for Molecules with Fused Hard Sphere Cores John M. Walsh and Keith E. Cubbins* School of Chemical Engineering, Cornell University, Ithaca, New York 14853 (Received: October 9, 1989; In Final Form: January 19, 1990)
The Wertheim thermodynamic perturbation theory (TPT), rigorously derived for molecules whose cores are made up of tangent hard spheres, is modified here for molecules in which the spheres overlap. The modified expression is derived by equating the second virial coefficient of TPT with that of scaled particle theory (SPT). By use of this approach, the number of spheres per molecule in the TPT is interpreted as an effective number of spheres. The effective number of spheres per molecule is then calculated from the shape factor of scaled particle theory (a),which depends on molecular geometry alone. The modified TPT is compared with published computer simulation results, and good agreement is found. Thus, the modified TPT gives a useful method of extending the original TPT to molecules with fused hard sphere cores without introducing empirical parameters.
I. Introduction Wertheim has presented a theory for the statistical mechanics of molecules that interact with hard core repulsive plus directional short-range attractive forces.1-8 Such attractive interactions are characteristic of covalent and hydrogen bonds. The theory, referred to as the thermodynamic perturbation theory (TPT), is based on a reformulation of the cluster expansion technique. Two separate density variables are introduced, one for the number of monomeric (nonbonded) molecules and the other for the total number of molecules. It turns out that using these two density variables allows a straightforward topological reduction of the cluster expansion series; such a reduction is intractable using only one density variable (total density). The result is a model which has the physically appealing characteristic that monomers are, in a sense, treated as a separate chemical species from dimers, trimers, and higher order oligomers. Also, the model allows direct calculation of the equilibrium number of monomers, which depends on the strength of bonding interactions between molecules and on the thermodynamic state of the fluid. Calculations using the TPT have shown excellent agreement with computer simulation results for a variety of fluids. Assuming only one bonding site per segment, Wertheim applied the theory to dimerization and obtained an equation of state for an equilibrium mixture of monomers and dimers.’ This application was tested in the zero monomer limit (complete dimerization) using Monte Carlo data obtained by Streett and Tildesley for the homonuclear diatomic tangent hard sphere fluid.9 Jackson et a1.I0 further developed the theory to provide a closed form expression for the equilibrium densities of monomers, dimers, trimers, etc. They found excellent agreement with Monte Carlo simulation results for the equilibrium fraction of monomers. Wertheim,8 and independently Chapman et al.,” used the TPT to derive an equation of state for hard chainlike molecules. Conceptually, chain molecules are formed from a mixture of spheres having one or two bonding sites each. Spheres having one bonding site become chain-end segments; spheres having two ( I ) Wertheim. M. S.J. Chem. Phys. 1983, 78,4619. (2) Wertheim, M. S. J. Chem. Phys. 1983, 78, 4625. (3) Wertheim, M. S. J . Star. Phys. 1984, 35, 19. (4) Wertheim, M. S. J. Stat. Phys. 1984, 35, 35. (5) Wertheim, M. S. J. Star. Phys. 1986, 42, 459. (6) Wertheim, M. S . J. Stat. Phys. 1986, 42, 477. (7) Wertheim, M. S. J. Chem. Phys. 1986, 85, 2929. (8) Wertheim, M . S. J . Chem. Phys. 1987, 87, 7323. (9) Tildesley, D. J.; Streett. W. 9. Mol. Phys. 1980, 4 / , 8 5 . (IO) Jackson, G.;Chapman, W. G.;Gubbins, K. E. Mol. Phys. 1988, 65, I. ( 1 I ) Chapman, W . G.;Jackson, G.; Gubbins, K. E. Mol. Phys. 1988,65, 1057.
bonding sites become internal segments. The equation of state is derived by considering the limit of complete bonding (zero monomer concentration). In the first-order theory, only fully flexible molecules are treated, and there is no distinction between a linear and a branched structure. The second-order theory accounts for the effect of chain conformation and branched structure but requires for its evaluation a fairly detailed knowledge of the triplet correlation function of the hard sphere fluid. The equation of state for hard chain molecules in the first-order TPT is given by pp/p = 1
+ u- 4y - 2y2 + ( u (1 - Y Y
Y 2-y
- 2) (1) 1-y
in which 0 = 1/ k T , P is the pressure, p is the number density ( p = N / V ) , u is the number of spheres per molecule, and y is the packing fraction Cy = V,p, in which V , is the volume of a molecule). For the hard sphere fluid the number of spheres per molecule is one ( u = l ) , and the equation of state of Carnahan and Starling is obtained.12 Equation 1 was tested with the Monte Carlo results of Dickman and HallI3 for fully flexible chainlike molecules of 4, 8, and I6 tangent hard spheres; good results were obtained for u 5 8. In ref 1 i , Chapman et al. also considered dispersion forces and interactions leading to molecular association between the chains. These results have not been tested against simulations. The TPT provides, for the first time, a consistent theoretical equation of state for both associating and nonassociating molecules with spherical or chainlike geometry. The application of this equation to real fluids requires the specification of molecular parameters1’ that characterize the diameter of the spheres ( d ) , the number of spheres per molecule ( u ) , and, in the case of associating molecules, the range or volume of interaction ( K H B ) and the interaction energy (eHB). Although the TPT does not provide a model for dispersion interactions, a perturbation theory could be used which would introduce a dispersion energy parameter (edis). Upon application of the T P T to real fluids, a difficulty arises in determining the parameter u. This parameter cannot be set equal to the number of segmental units that make up a real chainlike molecule because in real molecules the segmental units overlap considerably; such real molecules have bond lengths that are much less than their segment diameters. In the applications that have been made so far,I0J4 u is treated as an adjustable parameter, but (12) Carnahan, N . F.; Starling, K. E. J. Chem. Phys. 1969, 5 1 , 635. (13) Dickman, R.; Hall, C. K. J. Chem. Phys. 1988, 89, 3168. (14) Jackson, G.; Gubbins, K. E. Pure Appl. Chem. 1989, 61, 1021. Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. Fluid Phase Equilib., in press. Chapman, W. G.; Gubbins, K. E.; Jackson, G . ;Radosz, M. Ind. Eng. Chem. Res., in press.
0022-3654/90/2094-5 1 15$02.50/0 0 1990 American Chemical Society
Walsh and Gubbins
5116 The Journal of Physical Chemistry, Vol. 94. No. 12, I990
r\
Y-
o /
/-\
a
b
I
I
I
I
c
HOHON1'C'I.FAR T A N G E N T
Figure 1. Schematic representation of molecules from three classes of fluids.
unless care is taken, this approach can lead to unphysical results. As discussed in greater detail below, it is the aim of this paper to provide a specification for the parameter v for molecules of overlapping segmental units. In Figure 1 we show examples of three models of hard core molecules. Molecules composed of homonuclear tangent hard spheres (Figure la), for which eq 1 was derived, are characterized by the number of spheres ( v ) and by the sphere diameter (d). The molecular volume is calculated as the sum of sphere volumes ( V , = v r d 3 / 6 ) . The bond length is equal to the sphere diameter (I, = d ) and does not appear in the equation of state. A more realistic model is the united atom or fused hard sphere representation; this is shown for butane in Figure I b.I5 In this model, hydrogen atoms in methyl and methylene groups are collapsed upon their parent carbon atom. As shown, the butane molecule is modeled as four fused hard sphere segments. In general, fused hard sphere molecules are characterized by the number of spheres, the diameter of each sphere, the bond length between adjacent spheres, and the bond angle between each group of three adjoining spheres. Conceptually, the opposite extreme of the tangent hard sphere geometry is the convex body geometry. An example is the spherocylinder shown in Figure IC. A molecule is a convex body if any line whose ends are inside the body is also entirely inside the body a t every point along its length. By this criterion, a molecule of tangent hard spheres is an extreme example of a nonconvex body. By the same criterion, homonuclear hard sphere molecules are also nonconvex. The properties of convex molecules can be modeled by scaled particle theory (SPT).I6 The SPT equation, in terms of the compressibility factor, is given by
01
30
I
0 1
1
= RmSm/(3J")
(3)
in which R, is the mean radius of curvature, S, is the surface area, and V, is the volume of the molecule. The equations for calculating R,, S,, V,, and hence a in terms of geometrical functionals for convex bodies are given by Kihara." A straightforward numerical procedure for calculating these quantities is discussed below. As reviewed in refs 18 and 19, many modifications and extensions have significantly advanced the understanding of the liquid structure, virial coefficients, and equation of state of convex and nonconvex hard body fluids. Two extensions of S P T concern us here and are discussed further below. The first is a definition of the shape factor for fused hard sphere (nonconvex) molecules. The second is the hard convex body equation of state (more ( 1 5 ) Hsu, C. S.; Pratt, L. R.; Chandler, D. J . Chem. Phys. 1978.68.4213. ( 1 6) Gibbons, R. M. Mol. Phys. 1969, 17, 8 I . (17) Kihara, T. Intermolecular Forces. In Physical Chemistry, An Aduanced Treatise; Eyring. H.,Ed.; Academic Press: New York, 1980: Vol. S. p 663. (18) Gray, C. G.; Gubbins, K . E. Theory of Molecular Fluids; Clarendon Press: Oxford, in press; Vol. 2, Chapter 6. (19) Boublik, T.;Nezbeda, I . Collect. Czech. Chem. Commun. 1986, 51. 2301.
0 :
Figure 2. Wertheim's TPT and the improved scaled particle theory (also known as the hard convex body equation of state) versus literature-reported Monte Carlo results for flexible chainlike homonuclear tangent sphere molecules. The compressibility factor ((PIP)is shown as a function of packing fraction (J = V,p) for three chain lengths. The number of spheres per chain ( u ) is shown in the figure. The Monte Carlo results are from ref 13.
recently referred to as the improved scaled particle theory, ISPT).'9,20 lMost fused hard sphere molecules are nonconvex bodies for which the mean radius of curvature, and hence the shape factor, is not rigorously defined. However, Boublik and Nezbeda2' have shown that the shape factor for a nonconvex body can be calculated from that of an analogous convex body using the following approximation scheme. For the calculation of the mean radius of curvature, the real (nonconvex) body is replaced by an "equivalent" convex body which is generated by figuratively wrapping the molecule in a taut membrane. The taut membrane gives a convex surface over which the mean radius of curvature is defined. For example, the equivalent convex body for a diatomic homonuclear tangent sphere molecule (with sphere diameter d ) would be a prolate spherocylinder (of length d, excluding the hemispherical caps). Surface area and volume are defined for both convex and nonconvex geometry; therefore, the actual nonconvex geometry is used (without a membrane) in calculating these quantities. The shape factor for nonconvex bodies is then given by = Rc&j"/(3Vm)
in which the parameter a , the shape factor, characterizes the geometry of a convex body and is given by
I
0 3 0 4 pocking fraction (y=V&)
0:
(4)
in which Rcbdesignates the mean radius of curvature for the analogous convex body, and S, and V, are respectively the surface area ana volume of the actual nonconvex body. Boublik has shown that, within the SPT, a for linear fused hard sphere molecules is given rigorously by eq 4.22 For other fused hard sphere molecules, the second virial coefficient has been calculated numerically by using the exact expression and compared with predictions using eq 2 with a from eq 4.2' In all cases, eq 4 was found to give second virial coefficients accurate to within a few percent. Equation 4 has also been used to calculate a values for the ISPT which has been used as the nonconvex hard body reference term in several perturbation t h e ~ r i e s . ~ ~ - ~ O The SPT and its modifications (ISPT) provide an accurate account of convex and nonconvex hard body fluids for a: values up to about 2. but their accuracy falls off considerably with higher (20) Boublik, T. J . Chem. Phys. 1975, 63, 4084. (21) Boublik, T.; Nezbeda, I . Chem. Phys. Lett. 1977, 46, 315. (22) Boublik, T. Mol. Phys. 1981, 44, 1369. (23) Nezbeda, 1.; Labik. S. Mol. Phys. 1982, 47, 1087. (24) Kohler, F.; Quirke, N.; Perram, J. W . J. Chem. Phys. 1979. 7 / , 4128. (25) Tildesley, D. J. Mol. Phys. 1980, 41, 341. (26) Fischer, J. J . Chem. Phys. 1980, 7 2 , 5371. (27) Lombardero, M.; Abascal, J. L. F.; Lago, S . Mol. Phys. 1981, 42, 999. (28) Enciso, E.; Lombardero, M. Mol. Phys. 1981, 44. 725, (29) Lustig, R. Mol. Phys. 1986, 59, 173. (30) Bohn. M.; Fischer. J . : Kohler. F. Fluid Phase Equilih. 1986, 31, 233.
The Journal of Physical Chemistry, Vol. 94, No. 12, 1990 5117
Molecules with Fused Hard Sphere Cores TABLE I: Second Virial Coefficients from the Original TPT Homonuclear Tangent Hard Sphere Molecules Y BIINVm B 1
4.0
2 3
5.444“
~N vm ~ 4.0
~
I
5.5 7.0
6.853*
“Exact analytical result from ref 31. *Numerical result from ref 22.
elongati~ns.’’*~~ This is in contrast to the TPT which gives good results even for relatively elongated molecules (u = 8 for which a = 4.5). Results for the two approaches for flexible hard chain molecules are shown in Figure 2 . In this paper we present a modified TPT for the equation of state of molecules with fused hard sphere cores. The main motivation for this work stems from the accuracy of the original TPT equation for the relatively elongated chainlike tangent sphere molecules of Figure 2 and the possibility then of devising an equation of comparable accuracy for molecules with more realistic overlapping hard sphere cores. The basic idea is to interpolate between the geometry of tangent hard sphere molecules and that of convex molecules. We impose the restriction that the modified TPT must give the same result as the original TPT for the tangent hard sphere molecules for which the original TPT was derived. The modified TPT is derived in section 11, and in section 111 it is evaluated by using published computer simulation results. In applications of the Wertheim TPT to real fluids, this modified treatment provides a means for estimating v from a knowledge of molecular geometry. 11. Theory In general, the second virial coefficient is a function of temperature and a functional of the intermolecular pair potential. In principle, it can be calculated exactly from the well-known expression
in which the indicated bond is given by the Mayer f function. For the case of hard molecules, the intermolecular potential is harshly repulsive and the second virial coefficient reduces to a functional of molecular geometry alone. For our purposes, the second virial coefficient provides a means of characterizing the geometry of hard molecules and interpolating between the geometry of tangent hard sphere molecules and that of hard convex body molecules. According to the TPT, the second virial coefficient for tangent hard sphere molecules is given by
-B2- - -1( 3 u NV,
2
+ 5)
As shown in Table I, second virial coefficients predicted by eq 6 are in good agreement with exact analytical and accurate numerical results. According to SPT, the second virial coefficient for hard convex molecules is given by B , / ( N V , ) = 1 + 3a (7) Some time before the SPT was derived for convex bodies, lsihara showed that eq 7 is exact for convex molecules of arbitrary shape.3’ T o interpolate, the TPT is rewritten by replacing the number of spheres per molecule (u) with an effective number of spheres (u,). The second virial coefficient from T P T is set equal to that from SPT, which gives u, = 2a - 1 (8) This expression relates the effective number of spheres (u,) in the modified T P T to molecular geometry through the shape factor ( a ) . Thus, the problem of determining u, for molecules of fused hard spheres is reduced to the problem of determining a . (31) Isihara. A. J . Chem. Phys. 1951, 19, 397.
For many molecular geometries, the calculation of RcbrS,, V,, and hence a can be done analytically with relative ease. However, for some polyatomic fused hard sphere molecules, the calculation of S, and V , involves surface and volume integrals over regions of complicated geometry where three or more spheres overlap. Lustig has derived analytical expressions for the surface area and volume of these region^,^^,^^ but for polyatomic molecules, we found a numerical procedure (described below) to be more expedient. In implementing the numerical procedure though, it was quite useful to check our results using Lustig’s analytical expressions. Alejandre et al. have presented a numerical procedure for . ’ ~molecular evaluating a for molecules of arbitrary ~ h a p e . ~ ~The volume and surface area are calculated by using a Monte Carlo integration scheme, while the mean radius of curvature is calculated by direct Simpson’s rule integration. Specifically, the surface area is calculated in the following way. A random point is generated which lies on the surface of a particular sphere; then the position of this point is checked to determine whether it lies on an exposed part of the molecular surface or whether it lies inside an adjacent, overlapping sphere. This process is repeated for several million points on each sphere. The surface area is then calculated by (9)
in which N E kis the number of points which are found to lie on the exposed surface of sphere k , NTk is the total number of random points generated on the surface of sphere k , the quantity Sk is the surface area of sphere k (Sk= ndk2),and the summation is taken over all spheres. The molecular volume is calculated in a similar way by figuratively placing the molecule in a box, generating points within the box, and then checking the position of each point to determine whether it lies inside the molecular volume. The molecular volume is then
in which N1 is the number of points inside the molecule, NT is the total number of points, and V,, is the volume of the box surrounding the molecule. Before evaluating the modified T P T (eq 8) for various fused sphere molecules, it is necessary to verify that it gives the correct limiting result when applied to the tangent hard sphere molecules for which the original TPT was derived. For these molecules the components of the shape factor ( a ) are given by Rcb = y4(u
+ 1)d;
S, =
UT&;
V, = u(n/6)d3 (11)
From these expressions and eq 8, a and u, are given by
in which the desired result for u, is obtained (u, = u). As in the original TPT, the packing fraction is calculated as y = u ( r / 6 ) d 3 p . Thus, the modified TPT reduces to the original TPT when applied to tangent hard sphere molecules. In the next section, the modified TPT is evaluated for several fused hard sphere molecules. 111. Results and Discussion
In this section, the modified TPT is applied to the prediction of the compressibility factor ( p P / p ) for several molecular geom(32) Lustig, R. Mol. Phys. 1985, 55, 305. (33) Lustig, R. Mol. Phys. 1986, 59, 195. (34) Alejandre, J.; Chapela, G. A. Mol. Phys. 1987, 61, 11 19. (35) Alejandre, J.; Martinez-Casas, S.E.; Chapela, G . A. Mol. Phys. 1988, 65, 1185. (36) Jolly, D.;Freasier, B. C.; Bearman, R. J. Chem. Phys. Leu. 1977, 46,
75. (37) Nezbeda, I.; Boublik, T. Czech. J. Phys. 1977, B27, 1071. (38) Boublik, T.Czech. J. Phys. 1983, 833, 121. (39) Nezbeda, 1.; Reddy, M. R.; Smith, W. R. Mol. Phys. 1985, 55,447.
5118
Walsh and Gubbins
The Journal of Physical Chemistry, Vol. 94, No. 12, 1990
TABLE 11: Diatomic Heteronuclear Hard Molecules; Comparison of Monte Carlo (MC)" with Modified TPT (MTPT)* geometryC
d21dl
Ll4
~
I I II I1 Ill Ill IV IV
I .5 1.5 I .5 1.5 1.5 I .5 1.8 1.8 1.2 I .2
0.5
1.905 1.905 2.094 2.094 2.237 2.237 3.370 3.370 1.248 1.248
V V
0.5 0.75 0.75 1 .o
I .o
0.9 0.9 0.6 0.6
4
3
a
ye
Y
BPIPMC
PplPMrm
1.0334 1.0334 1.1113 1.1113 1.2312 1.2312 I .0932 1.0932 1.1217 1.1217
1.0668 1.0668 1.2226 1.2226 1.4624 1.4624 1.1864 I . 1864 1.2434 1.2434
0.209 44 0.366 52 0.209 44 0.366 52 0.209 44 0.366 52 0.209 44 0.366 52 0.209 44 0.366 52
2.66 6.04 2.75 6.45 2.93 7.13 2.68 6.34 2.78 6.51
2.57 5.92 2.71 6.42 2.91 7.19 2.67 6.31 2.72 6.49
Equation 1 with u replaced by u, = 2a - 1. The a values were obtained from ref 37 and verified by using the a Monte Carlo results of ref 36. numerical procedure described in the text. 'Schematic representations of molecular geometries are shown in Figure 3.
TABLE 111: Linear Homonuclear Hard Molecules; Comparison of Molecular Dynamics (MD)' with Modified TPT (MTPT)* geometryC V J d 3
VI VI VI VI1 VI1 VI1 VI1 Vlll Vlll
1.13 1.13 1.13 1.43 1.43 1.43 1.73 1.73 1.73
a
Y.
1.19 1.19 1.19 1.31 1.31 1.31 1.45 1.45 1.45
1.38 1.38 1.38 1.62 1.62 1.62 1.90 1.90 1.90
y 0.339 0.396 0.452 0.358 0.429 0.501 0.433 0.5 19 0.554
pPIpMD pP/pMrm 5.93 8.30 12.02 7.35 11.77 18.67 13.57 24.90 33.24
5.87 8.28 11.91 7.27 11.50 18.85 13.29 24.65 32.20
I'd
etries. Results for the modified TPT are calculated via eq 1, with u replaced by u, = 2a - 1 (eq 8). The results are presented here in roughly the order of increasing molecular complexity. Figure 3 shows schematically the molecular geometries considered. In general, computer simulation results for the compressibility factor are believed to be accurate to within a few percent. In all cases considered here, the density of the fluid is reported in terms of the packing fraction 0,= V , p ) . For reference purposes we have estimated the hard core packing fraction of liquid propane at its normal boiling point to be 0.410 (dimensionless). This was calculated from the experimentally measured molar volume of 13.2 mol L-'and with a value of the hard core diameter of 0.372 nm for each of three spheres in the united atom model. For homonuclear diatomic molecules, the modified TPT is compared with published Monte Carlo results in Figure 4. Several different bond lengths are shown over a range of density. The components of the shape factor ( a ) are given by the analytic expressions
'[
=- 1 2
3
+ - d s,
= xd2(1
+ L*); v,,,= -( 2 + 3L* - L*')d3 (1 3 ) 12 A
and the shape factor is given by ( I + L*)(2 a=
2
+ L*)
+ 3L* - L*3
(14)
in which L* is the reduced bond length (L* = L / d ) , and &b is the mean radius of curvature of the analogous hard convex body, as previously discussed. The geometry of these molecules, and
v; 1.
a
a Molecular dynamics results of ref 34. *Equation 1 with Y replaced by uC = 2 a - 1 . The a values were obtained from ref 34 and verified by using the numcrical procedure described in the text. 'All sphere diameters are equal, and all reduced bond lengths are equal ( L l d = 0.3979). Schematic representations of molecular geometries are shown in Figure 3.
RCb
ViI
VI
I
