J. Phys. Chem. 1983, 87, 1972-1978
1972
Perturbation Theories of Diatomic Fluids N. Qulrke' chemistry Department, Royal Holloway College, Egham, Surrey, TW20 OEX, England
and D. Tlldesley Chemistry Department, University of Southampton, Southampton, SO9 SNH, England (Received: September 20, 1982; In Final Form: November 4, 1982)
The structural and thermodynamic predictions of a number of different perturbation theories for diatomic liquids are compared with simulation results. For the diatomic Lennard-Jones (DLJ) liquid, zero-order perturbation theories based on nonspherical reference systems give a more accurate site-site distribution function, g&), than those on spherical reference systems. The spherical perturbation theories discussed do not improve at first order, but related expansions show the principal features in g&). For the DLJ liquids there is little to be gained by working in a site frame (SF)rather than a center frame (CF) for the zero-order f expansion but for repulsive DLJ liquids and hard-dumbell (HD) fluids the site frame is better. The reasons for this are discussed in the text.
1. Introduction
The diatomic liquids discussed in this paper are composed of homonuclear interaction site molecules. The total intermolecular potential between a pair of molecules is the sum of four interaction site potentials u,(r). We consider three different liquids: the diatomic Lennard-Jones (DW) where u,(r) = 4(r-12 - r*) (all distances are reduced by Q and energies by E where these parameters refer to the interaction site potential, uss(r)),a soft repulsive model where u,(r) = 4(r-12- rd + 0.25) for r < 2lI6 and u,(r) = 0 for r 1 2lI6, and the hard dumbell (HD) where u,(r) = m for r < 1 and u,(r) = 0 for r L 1. These are basic model molecular liquids with high symmetry and no multipolar interactions. Despite its simplicity, the DLJ potential provides a useful effective pairwise model for liquids N2,1 C12,2Br2,2and even C02.2 There have been extensive molecular dynamics simulations of the DLJ f l ~ i d l and -~ both str~cturall-~ and thermodynamic datal-" are available for a variety of elongations at liquid densities and temperatures. These simulations, although routine, are expensive (typically 1 h of CDC 7600 time to determine equilibrium properties). This would make the calculation of, for example, a coexistence curve for a molecular liquid a major undertaking. Furthermore, properties such as the light-scattering factor g25and sections through the total distribution function Gcc(R,Q1,Q2) are still difficult to calculate by simulatione6 For this reason a successful perturbation theory for simple molecular liquids, i.e., one which produced all the desired properties in a short time (- 1min of CDC 7600 time) would be useful. Such a description of these simple molecular fluids will provide a solid basis for theories of polar liquids and molecular liquid mixtures. In this paper we present new results for the structural and thermodynamic properties of D U , repulsive DW, and HD fluids obtained from a variety of perturbation expansions. The results are compared with simulation data at the same state points. The perturbation theories are based on spherical and nonspherical reference potentials and the expansions are generally truncated at zero and first order. In section 3 we introduce a novel expansion, in which part of a second-order term is included. 'Address correspondence to this author at the Department of Chemistry, University of Maine at Orono, Orono, ME 04469. 0022-365418312087- 1972$0 1.5010
Perturbation theories for molecular fluids are here divided into two kinds, depending on whether the expansion is made in a coordinate frame defined with respect to the center of force separation R or the site-site separation r. (In principle there are an infinite number of possible coordinate frames.) Our aim is to examine a number of current perturbation expansions and to test their ability to predict thermodynamic and structural properties. We are particularly interested in discriminating between site-frame (SF) and center-frame (CF) expansions for spherical reference systems and in comparing the structural predictions of zero-order SF expansions based on spherical and nonspherical reference systems. Throughout this paper we use the normal' system of reduced units for thermodynamic properties. Variables which describe the relative orientation of molecules in the center frame are {R,Ql,Qzjand in the site frame {r,01,02). The aim of a perturbation theory should be to describe the full molecular distribution function Gc,(R,Ql,Q2)= Gss(r,w1,w2),since these functions provide a complete description of the diatomic fluid. However, Gcc(R,Ql,Q2) is difficult to obtain from a simulation and we shall often calculate the angle-averaged functions in the site frame g&) and the center frame g,(R). The use of the site frame has an advantage in that internal energy, pressure, mean squared torque, and mean squared force depend on only a few of the spherical harmonic coefficients of G&-,01y2) whereas in the center frame they depend on an infinite sum of spherical harmonics of Gcc(R,Q1,Q2), which does not converge rapidly for large values of the site separation LE In the site frame, the configurational energy per particle U of a diatomic fluid is given by
(1) P. S. Y. Cheung and J. G . Powles, Mol. Phys., 30, 921 (1975). (2) K. Singer, A. Taylor, and J. V. L. Singer, Mol. Phys., 33, 1757 (1977). (3) W. B. Streett, and D. J. Tildesley, Proc. R . SOC. London, Ser. A , 353. 239 (1977). (4) N. Quirke and G . Jacucci, Mol. Phys., 45, 823 (1982).
(5)R. W. Impey, P. A. Madden, and D. J. Tildesley, Mol. Phys., 44,
1319 (1981). (6) P. Cummings, I. Nezbeda, W. R. Smith, and G. Morris, Mol. Phys.,
in press. (7) N. Quirke and D. J. Tildesley, Mol. Phys., 45, 811 (1982). (8) D. J. Tildesley, W. B. Streett, and D. S. Wilson, Chem. Phys., 36, 63 (1979).
0 1983 American Chemical Society
The Journal of Physical Chemistry, Voi. 87, No. 11, 1983 1973
Perturbation Theories of Diatomic Fluids
+
and the pressure P by
G,,(R,Q1,Q2)= (e-5v(R~n1~n2)/ (e-@v(R~nl~n2))cF]{go(R)
P=
47
5
1=0
+
(YIo(QI)Yoo(Q~) Yoo(QJ
Ylo(Q2))TW (9)
( I even)
The equation is written in reduced units and e W ( r )is the second spherical harmonic coefficient in the expansion of Gss(r,w1,w2).The calculation of the mean square force are discussed in ref 7. and torque In sections 2 and 3 we present an outline of the perturbation expansions used in this paper. Sections 5 and 6 present structural and thermodynamic predictions, and in section 7 we draw some conclusions about the use of the various perturbation expansions. 2. Expansions about Spherical Reference Potentials In this section we consider expansionsoploabout a spherical reference potential, VO1lJ2
where Tl(R)13involves integrations over the reference three-body distribution function, which can be approximated by the Kirkwood superposition approximation. The f S F is similar but contains odd-1 harmonic terms in the summation due to the lower symmetry of G,(r,wl,w2). The indirect correlation function Y,(R,Q1,Q2) is obtained from G,(R,Q1,Q2) by dividing out the low density limit of G. An alternative route is to divide the full G(R,Ql,Q2)by its zero-order result (which is exact at low densities and known to give a good representation of the structure at fluid densities13),producing the function Aw(R,Ql,Q2); then G,(R,Qi,Q2) = (go(R)e-@u(R,nlpn2) / ( e-@U(RVnlsn2) ) cF]Acc(R,Q 1,Q,) ( 10) and A, represents the modification of G, at high density. Expanding A,, to second order produces
A,, = 1 + Alee
pV0(R) = -In (e-@v(R*nl*nz))cF(center frame) (3)
pvo(r)= -In
(e-@"(r~wl~oz))SF (site frame)
(4)
+ A*,,
(11)
where
(...) denotes an unweighted average over all relative orientations with respect to R(CF) or r(SF). Since the development of these perturbation theories is the same in either frame, we continue with VOW. We can define a Mayer function, f (U(R,Ql,Q2);X)
f(U,,(R,Q,,Q,);~) = f(Vo(R)) + ~Lf(~C,(R,Ql,Q2)) - f(Vo(R))l (5) which is appropriate to the reference system for h = 0 and the full molecular potential for X = 1. Equilibrium p r o p erties of the fluid, X = 1,can be expanded in a Taylor series around those of the fluid X = 0. If we choose either eq 3 or eq 4 for Vo, the structure of the full molecular fluid is exact in the low density limit (the second virial coefficients of the real and reference fluids are the same). To obtain the structure of a molecular fluid from this type of theory, it is essential to expand Y,,(R,Ql,Q2) = G,,(R,Q2,,Q2)exp(+PU(R,Q1,Q2))rather than G(R,Q1,Q2) it~e1f.I~This expansion, to zero order, gives GcC(R,Qi,Q2)= go(R) exp(-P U(R,Q1, Q J ) / ( exp (- U(R$1
$2))
) CF (6)
and a similar expression in the site frame. In the center frame we call this expansion fCF and = (Gcc(R,Ql,Q2))CF = gcc(R)
(7)
the distribution of molecular centers. In the site frame the theory is denoted by fSF and gdr) = (Gss(r,u1,u2))~~ = g8&9
(8)
the site-site distribution function.12 The f expansions can be continued to higher order. In the center frame f C F predicts13 ~
(9)J. W. Perram and L. R. White, Mol. Phys., 24,1133(1972);28,527 (1974).See also. 'abstracts", Van der Waals Centenary Conference, 1973, pp 151-2;L. R. White, Doctoral Thesis, Australian National University, 1974. (10)W. R. Smith, Can. J. Phys., 62,2022 (1974). (11)D.Cook and J. S. Rowlinson, Proc. R. Soc. London, Ser. A , 219, 405 (1953). (12)I. Nezbeda and W. R. Smith, Chem. Phys. Lett., 81,79 (1981). (13)N.Quirke, J. W. Perram, and G. Jacucci, Mol. Phys., 39, 1311 (1980).
(13) To first order this expansion is identical with f c F but the second-order term contains both first- and second-order derivatives of Y,, with respect to A. (This expansion we call the upper-case F expansion to second-order, FCF.) If one neglects the second-order derivative of Y,,, AZccis simple to evaluate and becomes A2,*. Then the expansion FCF*,defined as Gcc(R,Qi,Q2)= go(R)(e-Bv,(R.n,,n2)/(e-8v,(R,n1,n2))CF1(1+ AIcc + Azo*) (14) provides an improved estimate of the angular correlation function (see section 5). It is interesting to note that the [1,2] Pade approximant to the Taylor series expansion of
A,,
1 + Ala" ** A,, = 1- A2,
when approximated by using A2,,* becomes G,,(R,fii,Q2)/ ( G c c ( ~ , ~ ~ , ~= ~ ) ) C F G,,'(R,Qi,Q2)/(Gcc1(R,QiQ2) )CF (16)
the reduced flCF expansion proposed by Smith and Nezbeda.14 (The superscript 1 indicates the first-order expansion results for G,,.) 3. Expansions about Nonspherical Reference Potentials The two expansions considered under this heading connect the properties of the DLJ fluids to those of (14)T.W. Melnyk, W. R. Smith, and I. Nezbeda, Mol. Phys., in press.
1974
The Journal of Physical Chemistry, Vol. 87, No. 11, 1983
hard-dumbell fluids whose properties are nknownn. They are both generalizations to molecular liquids of the Weeks, Chandler, and Andersen (WCA) e ~ p a n s i o n for ' ~ atomic liquids. For molecular liquids there is a greater flexibility in the division of the total potential into attractive and repulsive parts. In the site frame we divide each of the interaction site potentials u,(r) into a purely repulsive part ui,(r) = 4(r-12 r+ 0.25) for r < 21/6 and uia(r) = 0 for r > 2lI6 and an attractive part u&(r)= u,(r) - uis(r)for all r. Then with the bond length fixed we find hard-sphere potentials which are thermodynamically equivalent to ui,(r) by minimizing the blip function. This procedure defines an effective hard-core model. The free energy of the reference system can be obtained from a number of accurate equations of state16J7and, in addition, the HD g,(r) are available from the RISM integral equation. The details of the expansion WCA(SF) are given in ref 18. It is zero order in structure and first order in thermodynamics. In the center frame, sections through the full molecular potential U(R,Q1,Q2) are considered at each Q, and Q2. This radial potential is used to defiie a purely repulsive system, CP(R,Ql,Q2), by cutting the potential at its minimum and raising it by the well energy. This amounts to a WCA division of the potential at each orientation Ql, a,. The center-center correlation function of the repulsive reference fluid can be obtained from the corresponding harddumbell structure through a blip function expansion or directly from the f expansion for CT'(R,Ql,Q2). The details of this WCA(CF) expansion can be found in ref 19-21. Again it is zero order in structure and first order in free energy. Both expansions considered in this system take the structure of the full molecular fluid to be that of the corresponding soft-repulsive fluid and are high-temperature expansions.
Quirke and Tildesley
+ +
4. Numerical Techniques In the f expansions, the reference system potentials are of the form given by eq 3 and 4. In the case of f C F , go(R) is obtained by solving the Percus-Yevik (PY) approximation using Baxter's formulation.22 At the high, liquidlike densities considered, the PY results for go@) agree with the simulated go(R) to within 3%.13 To obtain g,(r) in the f C F expansion we use13 gSs(d= ( (go(R)exp[-PU(R,Ql,Qz)l/ (exp[-PU(R,Ql,Q2)1 )CFJ)SF
(17)
where the argument of the average is just the full Gee(R,Ql,Q2) at this order of perturbation theory and the angle average is carried out in the site frame. The triple integral in eq 17 is performed by using a (16)3Gauss-Legendre integration. In the case of fSF there has been some debate as to how to calculate go(r) (note to zero-order g,(r) = g d r ) ) . Smith and Nezbeda have shown that a variety of integral equations give quite different results for the HD system.23 To (15)D.Chandler, Annu. Reu. Phys. Chem., 29,441 (1978). (16)D. J. Tildesley and W. B. Streett, Mol. Phys., 41, 85 (1980). (17)T. Boublik and I. Nezbeda, Chem. Phys. Lett., 46,315 (1977). (18)D.J. Tildesley, Mol. Phys., 41,341 (1980). (19)F. Kohler, N. Quirke, and J. W. Perram, J. Chem. Phys., 71,4128 (1979). ~--.-,
(20)J. Fisher, J. Chem. Phys., 72,5371 (1980). (21)F. Kohler and N. Quirke in "Molecular Based Study of Fluids", J. M. Haile and G. A. Mansoori, Eds., American Chemical Society, Washington, DC, ACS Adu. Chem. Ser., in press. (22)F. Kohler, J. W. Perram, and L. R. White, Chem. Phys. Lett., 32, 42 (1975).
Figure 1. She-site distrlbution function 9 4 )for a fluid of homonuclear diatomic molecules, I = 0.5471, p = 0.524, T = 2.26. a: (0) md simulation, (---) f C F ,(-) f,. b: (0) md simulation, (-) WCA(SF).
avoid any ambiguity we have simulated the reference potential, Vo(r) (eq 4). We use a molecular dynamics (md) simulation of 256 atoms with a table look-up for the force (the grid size in the table is 0.01 and we use a second-order interpolation). One important difficulty with these potentials is that the reference system (Vo(r)at p and r ) may be in a two-phase gas/liquid region. The simulated g&) under these circumstances is characterized by unusually high peaks and by a failure to oscillate about unity.24 We have avoided the state points where this difficulty occurs (see section 7). The thermodynamic properties for the fCF and f S F fluid are calculated by using the appropriate equations from section 2 at the required order of the theory. The WCA(SF) nonspherical expansion requires fewer than 10 evaluations of the blip function for minimization, and all the integrals in the theory are one-dimensional. The structural prediction of WCA(SF) is gaa(r)and, to calculate pressures and mean squared torques, g&(r) harmonic coefficients are constructed from the SSAP appr~ximation.~ Alternatively the free energy at a number of state points can be differentiated with respect to p or T. 5. The Structure In comparisons of the structural predictions of theories there are a number of points which should be taken into consideration. (i) A good theory w ill produce the characteristic features of a particular site-site distribution function, particularly the rounded cusp at 1 + u (which becomes more pronounced with increasing I). (ii) The theory should produce the correct peak heights in the distribution functions. An overestimate in the first peak of g,(r) generally produces an internal energy which is too negative. In addition, g,(r) should follow the simulation result closely on the leading side of the first peak. Theoretical curves which are too high in the region 0.8-1.1~ give pressures which are significantly too high. In our view a theory can be said to produce good structure if it copes adequately with the two points mentioned above. In this section we shall examine the structural predictions of a number of theories: DLJ, repulsive DLJ, and HD fluids. In the main we concentrate on the site-site (23)W. R. Smith and I. Nezbeda, Chem. Phys. Lett., 82, 96 (1981). (24)G.Jacucci and N. Quirke, Nuouo Cimento, 58, 317 (1980).
The Journal of Physical Chemistty, Vol. 87, No. 11, 1983
Perturbation Theories of Diatomic Fluids
OOd'
"
"
10
1
15
"
'
(a)
1
1
1
'
"
20
I
'
'
"
25
'
'
1075
J
r
30
2O r
t:A
08
1
,
1.0
1
,
1.5
(b,
2.0
1
1
1
#
#
2.5
1
1
1
1
1
3.0
,
10
,
,
I
1
2
la1
I 14
I
I
I
16
18
,
]
PO
r
1
r
Figure 2. She-she distribution function g,(r) for a fluM of homonuclear diatomic molecules, I = 0.3292, p = 0.66, T = 1.79. a: (0) md simulation, (- - -) f,, (-) f,. b: (0)md simulation, (-) WCA(SF).
distribution function, gJr), but we shall also consider the distribution of molecular centers g,,(R). Diatomic Lennard-Jones Liquid. The simulated sitesite distribution functions for DLJ liquid are shown in Figure l a and 2a. Figure l a , 1 = 0.5471, corresponds to a bromine elongation and Figure 2a, 1 = 0.3292 corresponds I:d I , , , , , J to a nitrogen elongation. The results of the two zero-order 08 10 1 2 ibl 1 4 16 18 20 r spherical reference system perturbation expansions, fCF and Flgure 3. (a) Comparison of g,(r) from f , and WCA for a repulsive fSF, are included for comparison. DW fluid, I = 0.793, p = 0.422, T = 1.345 (CO,): (0)md simulation,* In Figure 1fCF produces a g,(r) which is too low on the (---) WCA, (-) f,. (b) Comparison of f,, WCA(SF), and a repulsive leading side of the first peak (pressure too low) and which DW simulation for I = 0.5471, T = 2.18, p = 0.524. Symbols as in part a. shows no sign of a cusp at r = 1.55. The range of the correlation function is too short. f s produces ~ a g,(r) which is the full DLJ fluid.2) The WCA theory is in excellent is too high on the leading side of the first peak and with agreement with the simulation; fSF is poor. In Figure 3b a high first maximum. It produces no cusp and is too for a bromine elongation, 1 = 0.5471, the perturbation short-ranged. The same inadequacies in fCF and fSF are theory based on the spherical reference system fails to apparent in Figure 2 at the lower elongation. It has been reproduce the structure between 1.1and 2.0. suggested12 that the rather unusual features in the f S F Although fSF is poor for this model, it is better than the reference potential at 1 + u should produce a cusp in the spherical expansion based on the centers frame. In Figure site-site distribution function. For DLJ potentials this 4a, fCF is used to calculate g,(r) from the total distribution does not appear to be the case. Where we were able to function G,,(R,O,,O,); the resulting site distribution produce the fSF results and the fcF result at the same state function is extremely poor. However, the center distripoints both theories were equally poor and failed to meet bution g,,(R) obtained from fCF is quite reasonable (see the criteria given at the beginning of this section. In Figure 4b). As we might expect, the f S F expansion based Figures l b and 2b the same state points are examined by on the site-site distance r gives reasonable gs,(r)'s while using the WCA(SF) nonspherical reference system perf C F based on R gives better gcc(R)'s. turbation theory. This theory correctly determines the This study of the repulsive DLJ fluids confirms the range of the correlation functions and produces cusplike conclusions reached for D U fluids, but now the fCF results features at 1 u. However, it is a high-temperature apfor g&) are worse than the fSF. proximation and therefore tends to overestimate the Hard-Dumbell Fluids. In this section we only consider leading edge of the first peak. Further, the theory relies the f expansions since the nonspherical perturbation upon the accuracy of the RISM equation. As can be seen theories take the hard-dumbell structure as their starting in Figure 2b, these two faults can combine to produce both point. The functions g&) from f s F have been published an overestimate of the leading edge of the first peak and by Smith and N e ~ b e d a . They ~ ~ are in reasonable agreean underestimate of the height of the cusp. There is clearly ment with the simulation data for homonuclear diatomics. room for improvement in the WCA(SF) but at the states At low densities they are better than the RISM integral which we examined it gave an overall more accurate picture equation g-(r). The extension of f S F to heteronuclear and of g,(r) than either fSF or fCF. ~~ polyatomic dumbells has not proved so s u c c e s ~ f u l .We Repulsive Diatomic Lennard-Jones Fluids. Our expehave tried to use fCF to calculate g,,(r). The results are rience in simulating this fluid is that the structure obtained is very close to that obtained for the full DLJ p ~ t e n t i a l . ~ shown in Figure 5a for a dumbell elongation 1 = 0.4 at a density p = 0.6. The results for r C 1 + 1 (inside the cusp) In Figure 3a we show a comparison of fSF and WCA(SF) for a highly elongated COz molecule. (We note that the (25) W. R. Smith and I. Nezbeda, Mol. Phys., in press. simulated distribution function included for comparison 1
+
1
,
1978
The Journal of F'oysical Chemistry, Vol. 87, No. 11, 1983
2o
r
Quirke and Tildesley
TABLE I: Thermodynamic Predictions for DLJ Liquids of the f ' g F , f Q C F , f l C F , F'cF*, and WCA(SF) Expansions"
c
1 = 0.3292, p = 0.66, T = 1.791
U
source 0
0
f
-16.9 -17.9 -16.6 F'cF* -17.7 WCA(SF) - 1 7 . 1 SIM -17.64 'SF
f'CF flCF
1
1 = 0.5471, p = 0.524, T = 2.26
1 = 0.793, = 0.422, T = 1.345
p
P
U
P
U
P
i1.4 -3.0 -3.7 -1.5 i3.1 -0.06
-12.9 -13.0 -11.5 -12.9 -12.4 -12.50
+3.2 -0.1 -0.8 +1.2 +3.6 +3.31
-10.7 -10.8 -9.2 -10.8 -10.5 -10.54
-2.1 -1.8 -1.8 -1.2 -0.0
+0.14
" Equations 1 and 2 integrated between r = 0 . 8 and 3.0 with appropriate long-range corrections. g,, repulsive DLJ used with f o S F ( I = 0.793);see text.
I
TABLE 11: Thermodynamic Predictions for a Repulsive DLJ Fluid from the f ' S F and f ' c Expansions ~ p =
source f'SF f 'CF
WCA(SF) S IM
2'oc
I
,
1 = 0.5471, 0.5239, T = 2.18
U
P
2.4 1.7 2.4 2.48
12.6 9.3 12.4 12.69
is indicative of a poorly convergent perturbation expansion which contains a considerable cancellation of the higher order terms. The F C F expansion to second order as approximated by eq 14 improves the situation. The height of g8,(r)at the first maximum is still too low but there is a cusp at r = 1.55. The results obtained from the [1,2] Pade to F C F eq 16, are very close to F C F * . Similar results for the carbon dioxide elongation I* = 0.793 are shown in Figure 6b. 6. Thermodynamic Results
0
0 0 0
14
Figwe 4. (a)Comparison of g,(r) obtained from f,, with simulation, for a repulsive DW, I = 0.5471, p = 0.524, T = 2.18: (0)md sim~ l a t i o n(-), ~ f,. (b) g,(R) from the same theory at the same state point. Symbols as in part a.
are poor and attempts to improve the perturbation theory by using the simulated g,,(R) rather than the reference system g,,(R) obtained from the Percus-Yevick equation do not improve the agreement. f S F appears to be a better theory of g,, (hard dumbell) than f c F . Again, Figure 4b shows that f C F gives a reasonable center distribution function g,,(R). High-Order f Theories. The perturbation theories discussed in sections 2 and 3 are all zero order in structure. I t is possible to extend the f expansions to higher order at the expense of numerical simplicity. In this section we consider the DLJ fluid and work with the f C F expansion (note that ~ C Fand f S F give comparable results for this system). In Figure 6a the simulation resulb are shown for DLJ fluid, 1 = 0.5471, p = 0.524, T = 1.23. The results to zero order are as discussed above. Surprisingly the results to first order in the f C F expansion are much poorer. This
In this section we look at the thermodynamic predictions of the various expansions, obtained both by integration of G,,, G, and directly from an expansion of the Helmholtz free energy. The accuracy of the results obtained from G , or G, provides another indication of the quality of the predicted structure. Table I gives the zero-order results for the internal energy and pressure of DLJ liquids obtained from the f , F, and WCA(SF) expansions. The WCA(SF) pressures are obtained by using the WCA(SF) g,,(r) plus the site superposition approximation SSA2.' The good result for the f S F pressure at I* = 0.5471 appears to be due to cancellation of errors in g,, (see Figure la). It is clear that neither of the f expansions is satisfactory whereas WCA(SF) is good except for the lowest elongation where the structural prediction is clearly poor. Table I1 gives the zero-order results for the repulsive DLJ liquid. There the f O s F expansion is an improvement on fOCF (see section 7). The higher order results for D U liquids in Table I show that the fmt-order terms are not an improvement but that F C F * produces some improvement in the thermodynamics as well as the structure. Note that these predictions can be further improved by the use of g,(r) obtained from a simulation of a molecular liquid in eq 1 and 2, leaving only the higher order coefficients to be predicted. In fact, Nezbeda and Smith%have obtained quantitative results for the pressure of repulsive DLJ liquids using this route. However, a strong motivation for developing perturbation theories is the desire to avoid the necessity for a simulation of the full molecular liquid. (26) I. Nezbeda and W.
R. Smith, J. Chem. Phys., 75, 4060 (1981).
The Journal of phvslcal Chemistry, Vol. 87, No. 11, 1983 1977
Perturbation Theories of Diatomic Fluids
2.07
0
1.5-
-
0
0 0 0
1.0L I
_
d -
I
0.5-
i
la)
1.2
4
1.6
1.4
,
r
0 0
08
10
12
(b)
16
14
18
r
20
Figure 6. Site-site dlsttibution function g,(f) for a fluM of homonuclear diatomic molecules. a: I = 0.5471, T = 1.23, p = 0.524; (-) focF, (A)f',, (A)F 2 w * , (0) md simulation. b: I = 0.793, T = 1.345, p = 0.422. Symbols as In part a.
TABLE 111: Helmholtz Free Energy for a Model of Nitrogen from the WCA(SF)'*and WCA(CF)'9*20a A,
T 3.0 3.0 2.0 1.55
p
0.7 0.6 0.7 0.7
A,
WCA(SF) WCA(CF) -3.11 -3.57 -6.02 -8.70
-3.06 -3.53 -6.02 -8.75
ACb -3.07 -3.53 -6.02 -8.71
a At T = 1.55, f 2 C F 1 4 (not F*cF)predicts -8.96; it is consistently too negative. The simulation data are from ref 4. 1 = 0.3292.
0.01-J
1.0
1.2
(bj
1.4
I
1.6
I
R
,
Figure 5. (a) Comparisonof g,(r) for a haMumbell system with the prediction of the f , perturbatlon theory, I = 0.4, p = 0.6: (0) Monte (b) g,(R) for the Carlo (mc) simulation,le (-) fcF(PY), (m) f,.&im). same system: (0) mc simuiationle (-) f,.
An alternative route to the thermodynamic properties which can be relatively insensitive to the underlying detailed structure2' is provided by direct expansion of the free energy. Table I11 shows that both WCA expansions (27) F. Kohler, W. Marius, N. Quirke, J. W. Perram, C. Hoheisel, and H. Breitenfelder-Manske, Mol. Phye., 38, 2057 (1979).
provide accurate predictions of the free energy, whereas f C F is too negative at high den~itie8.l~ Part of this may be due to the use of the Kirkwood SSA.28 It has been shown elsewhere that the free energy expansion leads to excellent and internal energies.18 The disadvantage of this procedure is that it is necessary to fit curves to the free energy data at several state points, in order to differentiate the free energy with respect to volume or temperature. 7. Conclusions
The principal thrust of this paper has been to compare a number of different perturbation theories for homonuclear diatomic fluids with simulation data for a variety of molecular elongations at a number of characteristic liquid (28) F. Kohler, private communication.
1978
The Journal of Physical Chemistry, Vol. 87, No. 11, 1983
Quirke and Tildesley
state point gives a fluid in the two-phase region and cannot be used to predict gsa(r). Note also that, for the f C F expansion, the predicted g,,(R) at high density is much improved by using the repulsive DLJ For any structural property other than g,,(R), the f C F expansion has severe difficulties in treating short-range potentials such as repulsive DLJ or HD. This is because the range of angular ordering is given by the range of the molecular Boltzmann factor (see eq 6). For repulsive and hard diatomic liquids this means that Gcc(R,Ql,Q2) = g,,(R) for R > 1 + 2 l l s or R > 1 + 1, respectively. The use of eq 17 to produce g&) at small r involves integrating over a wide range of R including those where the theory predicts no angular ordering but simulations show it clearly exists. This error produces the very bad estimates of g,(r) shown in Figures 3b and 4b. The f C F expansion is not improved by going to first order, but a number of related expansions, using essentially the same information, FCF and the [ 1,2] Pade of F C F (the reduced P C F expansion14),improve the results for g@(r). These higher order f type expansions are at present certainly no better than zero-order WCA(SF) expansions for g&). One motivation for continuing work on f expansions for DLJ liquids is that it is relatively easy Flgure 7. Reference fluid in the f , expansion for the state I = 0.793, to introduce multipolar interactions21 and thereby treat T = 1.345, p = 0.422 (A)reference potential and (0)reference force a wider class of molecular liquids. For dipolar liquids for the DLJ model; (A)reference potential and (0) reference force for Gcc(R,fll,Q2)is of great interest and the f expansions, in the repulsive DW model. contrast to the WCA theories, attempt to predict the whole function. Note however that good preliminary results for densities and temperatures. Most of the perturbation the thermodynamic properties of polar liquids have been theories discussed are zero order in structure and first or obtained from WCA(CF).21 second order in free energy. At zero order in structure the Although the WCA(SF) is very successful in predicting WCA(SF) expansion based on a nonspherical reference g&), it could be improved in two ways. Firstly, the hard system gives a better description of g&) than either f C F diatomic g,(r) from fSF could be used as input to WCAor f s ~The . f s F expansion gives a better description of g J r ) (SF). In some cases they are better than the RISM ges(r) than f C F for repulsive DLJ fluids and for HD fluids. which are used at present. Secondly, the high-temperature However, in the case of DLJ and repulsive DLJ fluids f S F approximation could be replaced by ORPA or EXP apdoes not reproduce the cusplike feature a t 1 + 1. proximation~.~ f S F and f C F expansions produce a g,(r) distribution We have concentrated on g&) in this paper, not the full function whose range is too short. This indicates that the G,(r,wl,wz). This is because g,(r) can be related directly packing fraction of the reference fluid is too small comto neutron and X-ray scattering and also gives the internal pared to the full system. energy and mean squared force of the fluid. The additional There is an interesting inconsistency in both the f S F and two harmonic coefficients required in the site frame to give f C F expansions. The site-site and center-center distributhe mean squared torque and pressure can often be aption functions produced by the DLJ and repulsive DLJ proximated from gea(r).7 fluids from molecular dynamics simulations are ~ i m i l a r . ~ Finally, we feel that a practical perturbation theory However, in the f s F expansion the reference potentials for should not rely on simulation data as input, since even the a repulsive and full DLJ potential are quite different and simulation of atomic systems takes far longer than the do not produce the same structure. In Figure 7 we have solution of an integral equation for a spherical potential. drawn the modulus of the force E&) = IVVo(r)j where Vo(r)is defined by eq 14 for the DLJ and repulsive DLJ (29) J. Fischer and N. Quirke, Mol. Phys., 38, 1703 (1979). model of COz. The g,(r) from Vo(r)for the repulsive DLJ (30) H. C. Andersen, D. Chandler, and J. D. Weeks, Adu. Chem. Phys., is shown in Figure 3a whereas Vo(r)for DLJ at the same 34, 105 (1976).
