Structure of hard ellipsoids: The Gaussian overlap model

Journal of Molecular Liquids 122 (2005) 1 – 10 www.elsevier.com/locate/molliq

Structure of hard ellipsoids: The Gaussian overlap model Ram Chandra Singh Department of Physics, Meerut Institute of Engineering and Technology, Baghpat Road-Bypass Crossing, Meerut-250002, India Received 20 January 2004; accepted 10 January 2005 Available online 27 April 2005

Abstract The equilibrium properties of isotropic fluids composed of hard ellipsoids of revolution represented by a Gaussian overlap model are studied using the integral equation theories. The Percus–Yevick and the hypernetted-chain integral equations have been devised for calculating the angular pair-correlation function for fluids of hard ellipsoids of a revolution. This model is computationally simple and shears some similarities with the widely used hard ellipsoids. The methods used involve an expansion of angle-dependent functions appearing in the integral equations in terms of spherical harmonics and the harmonic coefficients are obtained by an iterative algorithm. We have considered all the terms of harmonic coefficients which involve l indices up to less than or equal to 6. The numerical accuracy of the results depends on the number of spherical harmonic coefficients considered for each orientation-dependent functions. Ellipsoids with length-to-width ratios of 2.0 and 3.0 are considered and the harmonic coefficients are compared with the molecular-dynamics results for prolate ellipsoids. We find that both the Percus–Yevick and hypernetted-chain theories are in reasonable agreement with the computer simulation results. The pressure obtained from the virial and compressibility routes of these fluids have also been compared with the computer simulation results which shows that these integral equations are thermodynamically inconsistent. D 2005 Elsevier B.V. All rights reserved. Keywords: Hard ellipsoids; Gaussian overlap model; Isotropic fluids

1. Introduction The repulsive and attractive parts of the interaction potential seem to play different roles in the thermodynamics and structure of a fluid. It is generally accepted that, in atomic fluids, the repulsive part of the potential is mainly responsible for the structure of the fluid and the attractive interactions provide the cohesive energy but have little effect on the structure. Structure of equilibrium fluids is suitably characterized by the pair-distribution functions [1,2]. Among them, the pair-correlation function, g, plays the most important role. In molecular fluids, the intermolecular forces depend on the molecular orientations, vibrational coordinates, etc., in addition to intermolecular separation r. The dependence of the intermolecular forces on the molecular orientations leads to qualitatively new features in the fluid properties, when compared to atomic

E-mail address: rcsingh [email protected]. 0167-7322/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2005.01.003

fluids [3]. The angular correlations between the molecules are manifested in both the structural and thermodynamic properties. For molecular fluids whose molecules possess nonspherical shape, the pair-correlation functions (PCFs) are the lowest order microscopic quantities which, on one hand, contain information about the translational and the orientational structures of the fluid and, on the other, have direct contact with the underlying intermolecular (as well as with intramolecular) interactions [1–3]. Compared to atomic fluids for which solutions to the Ornstein–Zernike (OZ) equation have been obtained for a variety of pair potentials over wide ranges of temperature and density, our knowledge of the correlation functions of the fluid of nonspherical molecules is meager. In recent years, integral equation methods, which generally involve solving the OZ equation with the Percus–Yevick (PY) or hypernetted-chain (HNC) closure approximation, have been widely used to study the PCFs of classical simple fluids [2]. They have also been applied to isotropic fluids of nonspherical particles [4,5] and aniso-

2

R.C. Singh / Journal of Molecular Liquids 122 (2005) 1–10

tropic fluids of perfectly aligned molecules [6,7]. Zhong and Petschek [8] have developed an approach to extend Percus– Yevick approximation to anisotropic nematic fluids with partial spontaneous orientational order. PY and HNC theories for isotropic fluids with orientational dependent interactions have also been introduced by Pynn [9], Fries and Patey [4] and Singh and coworkers [10]. Several authors have calculated direct correlation functions (DCFs) in isotropic liquid crystals using theories of this kind [4,5,9,11–16]. The resulting PCFs were generally in reasonable agreement with available simulation data. An essential ingredient of the theory of fluid is the direct correlation function, which allows the determination of the mechanical stability (the compressibility) and the orientational stability (the Kerr constant) of the fluid, as well as the PCF to which it is linked by the OZ equation [2]. The DCF of isotropic phases has now been well studied theoretically [9,16–22] and by simulation [23], few studies have been devoted to the DCF of nematic mesophases [24]. Recently Phuong et al. [25] have proposed a method to calculate the DCF in uniaxial nematic fluids without approximations from computer simulations. Density-functional theories in general provide a suitable framework for studying all the various liquid crystal phases of nonspherical particles in particular [26–28]. The DCF plays a key role in density-functional theories, since it is the second functional derivative of the excess free energy with respect to the local density [2,27]. The DCFs of the isotropic phase are used as input informations in the densityfunctional theory. The freezing transitions in complex fluids can be predicted reasonably well with the density-functional method if the values of the DCFs in the isotropic phase are accurately known. Much effort has therefore been devoted to the investigation of DCFs in molecular fluids. The principal purpose of this paper is to compare integral equation theories with molecular-dynamics (MD) calculations for fluids of hard ellipsoids of revolution (HER). The motivation for studying these systems stems from the fact that they capture some of the basic features of the real ordered phases. For example, a system of HER is found [29] to exhibit four distinct phases, viz., isotropic fluid, nematic fluid, plastic solid and ordered solid. Systems of ellipsoids serve as model systems which exhibit generic properties of nematic liquid crystals. They are studied for the general insight they can give into the physics of liquid crystals, and also because simulations of ellipsoids are computationally much cheaper than simulations of a comparable number of more realistic molecules. However, fluids of ellipsoids are also fascinating from a much more general point of view: They are simple systems where purely entropic effects generate an amazing variety of different structures, and lead to unusual longrange static and dynamic properties. The paper is arranged as follows. In the next section, we briefly describe the hard ellipsoids of revolution model. The

procedure for solving the OZ equation using PY and HNC closure relations and a method in which the PCFs are expanded in a series of spherical harmonics are discussed in Section 3. In Section 4, we discuss the results for the PCFs of hard ellipsoids of revolution and make a thorough evaluation of the PY and HNC theories. We also compare our results with those of the computer simulations wherever the latter are available.

2. The HER model One of the simplest (convex) repulsive interaction methods is the hard Gaussian overlap model, which can be considered as an approximation for the excluded volume of HER. The Gaussian overlap model was originally proposed by Berne and Pechukas [30] as an attempt to mimic the short-ranged repulsive interactions between elongated molecules. Each molecule is represented by a Gaussian distribution, and the repulsive interaction between a pair of molecules is assumed to be defined in terms of the overlap of their Gaussians. The potential energy of the interaction of a pair of HER is represented as uð1; 2Þ ¼ uðr12 ; X1 ; X2 Þ
l for r12 bDðrˆ 12 ; X1 ; X2 Þ ¼ 0 for r12 zDðrˆ 12 ; X1 ; X2 Þ

ð1Þ

where rˆ12 is a unit vector along r 12 = r 2  r 1. D(rˆ12, X 1, X 2) is the distance of closest approach of the two molecules in contact with orientations X 1 and X 2 in a direction given by rˆ12. For HER, D(rˆ12, X 1, X 2) is taken to be that given by the Gaussian overlap model of Berne and Pechukas [30], i.e. Dðrˆ 12 ; X1 ; X2 Þ " #1=2 ðrˆ 12d eˆ 1 Þ2 þ ðrˆ 12d eˆ 2 Þ2  2vðrˆ 12d eˆ 1 Þðrˆ 12d eˆ 2 Þðeˆ 1d eˆ 2 Þ ; ¼ d0 1 v 1  v2 ðeˆ 1 eˆ 2 Þ2

ð2Þ eˆ1 and eˆ2 are unit vectors along the symmetry axes of two interacting hard ellipsoids, d 0 =2b, the minor axis of the ellipsoid and v is a shape (anisotropy) parameter given by v¼

x20  1 ; x20 þ 1

where x 0 (= 2a/2b) is the length (major axis)-to-width (minor axis) ratio of these ellipsoids. The molecular volume and are, respectively, given  the packing fraction  as v ¼ p6 d03 x0 and g ¼ p6 qf x0 where q f* =q fd 03. Eq. (2) has a closed form but is not an exact representation of the ellipsoid overlap. This model includes as limiting cases the hard-spheres, hard-platelets, and hard needles systems. All these systems are of physical interest because they

R.C. Singh / Journal of Molecular Liquids 122 (2005) 1–10

represent primitive models for real fluids, solids, plastics, and liquid crystals.

3. Integral equations for correlation functions The single-particle density distribution q(1) is defined as X  qð1Þ ¼ qðr; XÞ ¼ dðr  ri ÞdðX  Xi Þ ; i¼1

where ri and X i give the position and the orientation of ith molecule, the angular bracket represents the ensemble average, and the d the Dirac delta function, is constant independent of position and orientation for an isotropic fluid. It therefore contains no information about the structure of the system. The structural information of an isotropic fluid is contained in the two-particle density distribution q(1,2) that gives the probability of finding simultaneously a molecule in a volume element dr1dX 1 centered at (r1, X 1) and a second molecule in a volume element dr2dX 2 centered at (r2, X 2). q(1,2) is defined as qð1; 2Þ ¼ qðr1 ; X1 ; r2 ; X2 Þ X      dðr1  ri ÞdðX1  Xi Þd r2  rj d X2  Xj : ¼ ipj

The pair-correlation function g(1,2) is related to q(1,2) by the relation gð1; 2Þ ¼

qð1; 2Þ : qð1Þqð2Þ

3

can be expressed as integral over g(1,2) or its spherical harmonic coefficients. The values of PCFs as a function of intermolecular separation and orientations at a given temperature and pressure are found either by computer simulations or by semianalytic approximate methods [3]. In the latter approach, one solves the Ornstein–Zernike equation, cð1; 2Þ ¼ hð1; 2Þ  cð1; 2Þ Z ¼ qf cð1; 3Þ½cð2; 3Þ þ cð2; 3Þ dx3 ;

ð3Þ

where q f is the number density of the fluid, h(1,2) = g(1,2)  1 and c(1,2) are, respectively, the total and direct PCFs, with suitable closure relations such as PY integral equation, HNC equation, mean spherical approximation (MSA), etc. Approximations are introduced in the theory through these closure relations. In Eq. (3), i =x i indicate both the location r i of the centre of the ith molecule and its relative orientation X i , described by the Euler angles h, u, and w. c(1,2) is intrinsically a shorter-ranged function than h(1,2), acting as a kernel in the above relation. Although experiments and simulations do not provide a direct route to this function, it is of equal importance to h(1,2) in the statistical mechanics of liquids, and there has been a dramatic growth in the interest in c(1,2) in recent years. This is because of the rapid development of densityfunctional theories of fluids [26–28]: c(1,2) is simply related to the excess free energy F ex of the fluid by functional differentiation with respect to the local density q(1) =q(r 1,X 1): cð1; 2Þ ¼ d2 ð  F ex =kB T Þ=dqð1Þdqð2Þ;

Since for the isotropic fluid q(1) = q(2) = q f = hNi/ V where hNi is the average number of molecules in volume V, q2f gðr; X1 ; X2 Þ ¼ qðr; X1 ; X2 Þ; where r =(r2  r1). In the isotropic phase, q(1,2) depends only on the distance |r2  r1| = r, the orientation of molecules with respect to each other and on the direction of vector r(rˆ =r/r is a unit vector along r). Structure in atomic and molecular fluids is described by the paircorrelation function g(r 12, X 1,X 2) =g(1,2), where r1 and r2 are the centre of mass coordinates of particles 1 and 2, and X 1 and X 2 are unit vectors defining the orientations. Thus, the pair-correlation function g(1,2) enables one to compute the thermodynamic properties and in addition describes the fluid structure and thus is of fundamental importance to the theory of equilibrium properties of molecular fluids. Because of the large number of variables involved, the complete determination of full g(1,2) is not a simple task, and not yet accessible experimentally to the best of my knowledge. One usually, therefore, either considers spherical harmonic coefficients of g(1, 2) or value of g(1,2) as a function of intermolecular separation r 12 for some fixed angular orientations. The various equilibrium properties

where k B is Boltzmann’s constant and T the temperature. This expression leads to a variety of theories of both homogeneous and inhomogeneous fluids, mostly based on assumptions regarding c(1,2) in the system of interest, or in some reference system used as the basis of a perturbation treatment. The PY and HNC integral equation theories are given by the OZ equation coupled with the closure relation cPY ð1; 2Þ ¼ f ð1; 2Þ½1 þ cð1; 2Þ ;

ð4aÞ

and cHNC ð1; 2Þ ¼ hð1; 2Þ  ln½1 þ hð1; 2Þ  buð1; 2Þ;

ð4bÞ

respectively. Here f(1,2)= exp[bu(1,2)]  1, b = (k BT )1, u(1,2) is a pair potential energy of interaction. We find it convenient to separate the overlap and nonoverlap regions and choose
 1  cð1; 2Þ for overlap region ð aÞ cð1; 2Þ ¼ : cCL ð1; 2Þ for non overlap region ðbÞ ð5Þ

4

R.C. Singh / Journal of Molecular Liquids 122 (2005) 1–10

Since in the overlap region g(1,2) =1 +h(1,2) is zero, Eq. (5a) is exact. For c CL(1,2), one uses a closure relation given by Eq. (4a,b) or any other form. In order to reduce and solve Eq. (3) with a closure relation, it is convenient to expand all angle-dependent functions in spherical harmonics either in space-fixed (SF) frame or in a body-fixed (BF) frame. For example, the expansion of direct PCF in a BF-frame is written as X cð1; 2Þ ¼ cðr12 ; X1 ; X2 Þ ¼ cl1 l2 m ðr12 ÞYl1 m ðX1 ÞYl2 mP ðX2 Þ;

and

where m = m. To obtain the iterative solution of Eqs. (3) (4a,b) (5), it is necessary to expand the closure relations in rotational invariants. The expansion for HER can be easily performed following the methods given in our earlier papers [14,15]. For the PY and HNC approximations, one obtains the final form of the direct PCFs as [15] X 1V l2V m V cl1 l2 m ðr12 Þ ¼ 4pdl000  cl1Vl2V m Vðr12 ÞcCLV l1V l2V mVðr12 Þ

For hard Gaussian overlap model, c lHNC (r 12) reduces to 1l 2m NOV c lHNCV (r ) where bu (r ) becomes zero. 12 l 1l 2m 12 1l 2m In any numerical calculations, we can handle only a finite number of the spherical harmonic coefficients for each orientation-dependent function. The accuracy of the results depends on this number. As the anisotropy in the shape of molecules (or in interactions) and the value of fluid density q f increases, more harmonics are needed to get proper convergence. We have found that the series get converged if we truncate the series at the value of l indices equal to 6 for molecules with x 0V3.0 [15]. Though it is desirable to include higher-order harmonics, i.e., for l N 6, it will increase the computational time manifold. One can use the data of the harmonics of the DCFs for freezing transitions where only low-order harmonics are generally involved. The only effect the higher-order harmonics appear to have on these low-order harmonics is to modify the finer structure of the harmonics at small values of r whose contributions to the structural parameters (to be defined below) are negligible. The iterative numerical solution can be carried out in a manner described elsewhere [14,15]. An appropriate grid width Dr =0.01 is chosen in configuration space and the various functions are tabulated on M = 1024 grid points, the step size in Fourier space being Dk ¼ MpDr . All onedimensional integrals can be conveniently calculated using the trapezoidal rule. The correlation function expansions

l1 l2 m

b

t

l1Vl2Vm V

V CL  All11Vll22Vm m ðr12 Þ þ cl1 l2 m ðr12 Þ;

ð6Þ

where All11Vll22Vmm V ðr12 Þ ¼

1=2

1 X ð2l1 þ 1Þð2l2 þ 1Þð2l1Vþ 1Þð2l2Vþ 1Þ 4p LLVM ð 2L þ 1Þð 2LVþ 1Þ    Cg ð l1 l1VL; 000ÞCg ð l2 l2VL V; 000ÞCg l1 l1VL; P m mVM P   VM fLLVM ðr12 Þ: ð7Þ  Cg l2 l 2VL V; mm P

Here C g(l 1l 2l;m 1m 2m) are the Clebsch–Gordon coefficients, c l 1lCL (r 12) represents either the HNC or PY 2m approximation, cHNC l1 l2 m ðr12 Þ

X 1 ð2l1Vþ 1Þð2l1Wþ 1Þð2l2Vþ 1Þð2l2Wþ1Þ 1=2 ¼ 4p ð2l1 þ1Þð2l2 þ 1Þ l V l Wl V l W

cPY l1 l2 m ðr12 Þ

X 1 ð2l1Vþ 1Þð2l1Wþ 1Þð2l2Vþ 1Þð2l2Wþ 1Þ 1=2 ¼ 4p ð2l1 þ 1Þð2l2 þ 1Þ l1Vl1Wl2Vl2W X Cg ðl1Vl1Wl1 ; 000ÞCg ðl2Vl2Wl2 ; 000Þ Cg ðl1Vl1Wl1 ; mVmWmÞ mVmW

h i l1Wl2WmW f m Vm Wm ð r Þ 4pd þ c ð r Þ : Cg l2Vl2Wl2 ; P l V l V mV 12 12 l Wl WmW 1 2 000 1 2 P P 



ð10Þ

1 1 2 2

 Cg ðl1Vl1Wl1 ; 000ÞCg ðl2Vl2Wl2 ; 000Þ X   Cg ðl1Vl1Wl1 ; mVmWmÞCg l2Vl2Wl2 ; P m Vm  P Wm P mVmW Z l

B h hl1Wl2WmWðr12 VÞ VÞ  cl1Vl2V mVðr12 Br12 V r12 i V Þ dr12 V  buNOV þ buNOV l1Vl2VmV ðr12 l1 l2 m ðr12 Þ; 

ð8Þ

cHNCV l1 l2 m ðr12 Þ

X 1 ð2l1Vþ 1Þð2l1Wþ 1Þð2l2Vþ 1Þð2l2Wþ 1Þ 1=2 ¼ 4p ð2l1 þ 1Þð2l2 þ 1Þ l Vl Wl Vl W 1 1 2 2

 Cg ðl1Vl1Wl1 ; 000ÞCg ðl2Vl2Wl2 ; 000Þ X   Cg ðl1Vl1Wl1 ; mVmWmÞCg l2Vl2Wl2 ; P m Vm Wm  P P 

mVmW Z l r12

hl1Wl2WmWðr12 VÞ

 B  V ;  cl1Vl2V m Vðr12 V Þ dr12 Br12 V

ð9Þ

Fig. 1. Pair-correlation function of the centre of mass g(r) for x 0 =2.0, g =0.3702. The solid and dashed curves are, respectively, HNC and PY results. The solid circles are MD results of Talbot et al. (Ref. [12]).

R.C. Singh / Journal of Molecular Liquids 122 (2005) 1–10

5

Fig. 2. The spherical harmonic coefficient g 200 (r)/4k in BF frame for x 0 =2.0, g =0.3702. The curves are the same as in Fig. 1.

Fig. 4. The spherical harmonic coefficient g 240 (r) / 4k in BF frame for x 0 =2.0, g =0.3702. The curves are the same as in Fig. 1.

include all terms for which l 1, l 2, mV6 involving a total of 30 independent projections. The iterative cycle is started by making initial guesses for the coefficients c l 1l 2m (r 12). The iteration proceeds as follows:

large differences between successive iterates and hence prevent divergence.

cl1 l2 m ðr12 ÞðiÞ Y cl1 l2 m ðr12 Þ Y cl1 l2 m ðr12 ÞðnewÞ ;

4.1. Pair-correlation functions

where the superscript (i) denotes the ith iteration. After both steps c l 1l 2m (r 12)(i) and c l 1l 2m (r 12)(new) are compared, and if they agree with the desired accuracy (here we take the maximum difference to be less than 103), the iteration is terminated. If they do not satisfy this criterion then the (i + 1)th approximation is obtained by mixing the bithQ and the bnewQ values according to the equation

where a is the mixing parameter adjusted in the range 0bab1. The function of the mixing parameter is to avoid

The HNC and PY results for the values of paircorrelation function g(r) = 1 +h 000(r) / 4k of the BF-frame are compared for x 0 = 2.0 at g =0.3702 in Fig. 1 with the MD results of Talbot et al. [12]. The spherical harmonic coefficients g 200(r)/4k, g 220(r)/4k, g 240(r) / 4k, g 400(r) / 4k, g 440(r)/ 4k, g 221(r) / 4k of the HNC and PY equations are also compared with the MD results for x 0 = 2.0 at g =0.3702 in Figs. 2–7. It can be seen from Figs. 1–7 that at g = 0.3702 for x 0 = 2.0, both the PY and HNC theories are in reasonably good agreement with the MD calculations. The ellipsoids overlap model of Talbot et al. [12], however, differs from the Gaussian overlap model of Eq.

Fig. 3. The spherical harmonic coefficient g 220 (r)/4k in BF frame for x 0 =2.0, g =0.3702. The curves are the same as in Fig. 1.

Fig. 5. The spherical harmonic coefficient g 400 (r) / 4k in BF frame for x 0 =2.0, g =0.3702. The curves are the same as in Fig. 1.

OZ

4. Results and discussion

PY

cl1 l2 m ðr12 Þðiþ1Þ ¼ acl1 l2 m ðr12 ÞðiÞ þ ð1  aÞcl1 l2 m ðr12 ÞðnewÞ ;

6

R.C. Singh / Journal of Molecular Liquids 122 (2005) 1–10

Fig. 6. The spherical harmonic coefficient g 440 (r)/4k in BF frame for x 0 =2.0, g =0.3702. The curves are the same as in Fig. 1.

Fig. 8. Pair-correlation function of the centre of mass g(r) for x 0 =3.0, g =0.3702. The curves are the same as in Fig. 1.

(1) used by us. It is therefore reasonable to expect that the difference in the values of these two results may partly be due to the difference in the overlap model. For x 0 = 3.0, the HNC and PY results for g(r) =1 +h 000(r)/ 4k are plotted at g =0.3702 in Fig. 8. In Fig. 8, we also plot the MD result of Talbot et al. [12]. The other spherical harmonic coefficients of HNC and PY equations for x 0 = 3.0 at g = 0.3702 are compared in Figs. 9–14. Comparing Figs. 1–7 for x 0 = 2.0 and Figs. 8–14 for x 0 = 3.0 at g = 0.3702, we see that for similar packing fractions, the curves retain the same general shape but the structural features become broader as x 0 increases. For x 0 = 3.0 (Figs. 8–14), the agreement between the integral equation theories and MD result is similar to that found for x 0 =2.0 (Figs. 1–7) at g =0.3702. It is seen from Figs. 1 and 8 that g(r) exhibits a pronounced maximum at a scaled separation r* (= r/d 0) just greater than c1.25 corresponding to adjacent ellipsoids lying parallel to each other. There is a second, weaker maximum at r* c2.5 and this may be associated with a shell of next nearest neighbours, again with symmetry axes

parallel. The structure in the g(r) is lost for larger separations as expected for an isotropic liquid and becomes less pronounced as x 0 is increased. This is most likely due to decreasing tendency of ellipsoids to form parallel configurations. This effect is also apparent in the g 220(r) harmonic coefficients. It can be seen from Figs. 1–14 that for prolate ellipsoids, both the HNC and PY theories are in generally good agreement with MD calculations. Although the HNC theory does appear to have slight edge in overall accuracy than PY, it is obvious that the different approximations predict very different results as the density is increased. The HNC and PY equations for HER have also been solved by Perera, Kusalik and Patey [4]. Both the HNC and PY theories were found to give results for PCFs which are in reasonable agreement with the MD results for prolate molecules. In the figures plotted here, we do not show the results of Perera, Kusalik and Patey [4] due to the difference in the two ellipsoid overlap models. Phuong et al. [25] have also calculated the DCF in PY approximation for fluids of uniaxial ellipsoids by computer simulation. In their PY

Fig. 7. The spherical harmonic coefficient g 221 (r)/4k in BF frame for x 0 =2.0, g =0.3702. The curves are the same as in Fig. 1.

Fig. 9. The spherical harmonic coefficient g 200 (r) / 4k in BF frame for x 0 =3.0, g =0.3702. The curves are the same as in Fig. 1.

R.C. Singh / Journal of Molecular Liquids 122 (2005) 1–10

Fig. 10. The spherical harmonic coefficient g 220 (r) / 4k in BF frame for x 0 =3.0, g =0.3702. The curves are the same as in Fig. 1.

calculations, the spherical harmonics expansion is extended up to l = 4. 4.2. The equation of state Once the PCFs are known, they can be used to calculate pressure from virial and compressibility equations in a straightforward way [3]. In an exact calculation both routes must give the same result, but this is not true for approximate theories such as HNC and PY approximations. The compressibility equation of state can be expressed in the form

7

Fig. 12. The spherical harmonic coefficient g 400 (r)/4k in BF frame for x 0 =3.0, g =0.3702. The curves are the same as in Fig. 1.

Expanding Eq. (11) in spherical harmonics in SF-frame and after simplification, we get Z bPC 1 q ð0Þ ¼1 ð12Þ cˆ ½qV dqV; q 0 00 q where ð0Þ cˆ 00

qf ffi ¼ pffiffiffiffiffi 4p

Z

2 r12 c000 ðr12 Þdr12 :

ð13Þ

In general for isotropic fluids, the virial equation of state can be expressed in the form

ð11Þ

Z Z bPV 2p 3 Buðr12 ; X1 ; X2 Þ bq dX1 dX2 dr12 r12 ¼1 3 Br12 q  gðr12 ; X1 ; X2 Þ; ð14Þ

Fig. 11. The spherical harmonic coefficient g 240 (r)/4k in BF frame for x 0 =3.0, g =0.3702. The curves are the same as in Fig. 1.

Fig. 13. The spherical harmonic coefficient g 440 (r)/4k in BF frame for x 0 =3.0, g =0.3702. The curves are the same as in Fig. 1.

  Z B bPC ¼ 1  qf dX1 dX2 dr12 cðr12 ; X1 ; X2 Þ: Bq

8

R.C. Singh / Journal of Molecular Liquids 122 (2005) 1–10

Fig. 14. The spherical harmonic coefficient g 221 (r)/4k in BF frame for x 0 =3.0, g =0.3702. The curves are the same as in Fig. 1.

where P V represents the virial pressure. Expansion of Eq. (14) in spherical harmonics yields X X bPV 2p 1 q ¼1þ ð2l1Vþ 1Þð2l1Wþ 1Þ 3 ð4pÞ2 l Vl Wl Vl W mVmW q 1 1 2 2 1=2  ð2l2Vþ 1Þð2l2Wþ 1Þ Cg ðl1Vl1W0; 000Þ  Cg ðl2Vl2W0; 000ÞCg ðl1Vl1W0; mVmW0Þ    Cg l2Vl2W0; P m Vm W0 P Z Bfl1Vl2VmV ðr12 Þ 3  yl1Wl2WmWðr12 Þr12 dr12 ; Br12

Fig. 15. Pressure as a function of fluid density for x 0 =2.0. The solid circles represent the Monte Carlo result of Mulder and Frenkel (Ref. [29]), the dash-dotted and the dashed curves are, respectively, PY and HNC results. Here, V and C denote the virial and compressibility results, respectively.

ð15Þ

where g ðr12 ; X1 ; X2 Þ ¼ ½1 þ f ðr12 ; X1 ; X2 Þ yðr12 ; X1 ; X2 Þ:

plays an important role. Here P l is a Legendre polynomial of degree l and angles refer to a space-fixed z-axis. Note that (0) (0) cˆ 00 is related to the isothermal compressibility and cˆ 22 and the higher-order coefficient to the freezing parameters. The

ð16Þ

In Figs. 15 and 16, we plot the virial and compressibility pressures as a function of density for x 0 = 2.0 and 3.0, respectively. In Figs. 15 and 16, we also plot the computer simulation results of Mulder and Frenkel [29]. One may, however, note that the simulation results given in these figures are for the exact HER overlap which is not identical to the Gaussian overlap model used by us. It may therefore be possible that the results found for the Gaussian overlap model may not equally match with the simulation results for all the harmonic coefficients for the PCFs. 4.3. The structural parameters In a theory of freezing of molecular fluids into nematic phase the structural parameter cˆ l (0) defined as [28] 1l 2 ð0Þ

cˆ l1 l2 ¼ ð2l1 þ 1Þð2l2 þ 1Þqf

Z

 Pl1 ðcosh1 ÞPl2 ðcosh2 Þ;

dr12 dX1 dX2 cðr12 ; X1 ; X2 Þ Fig. 16. Pressure as a function of fluid density for x 0 =3.0. The curves are the same as in Fig. 15.

R.C. Singh / Journal of Molecular Liquids 122 (2005) 1–10

Fig. 17. The structural parameter cˆ LL(0) for HER at x 0 =3.0. The solid and the dashed curves are, respectively, PY and HNC results. (0) (0) quantities cˆ 22 and cˆ 44 are found to be very sensitive to the approximations involved in a given integral equation theory. In Fig. 17, we compare the values of the structural (0) (0) parameters cˆ 22 and cˆ 44 found from these theories for x 0 = 3.0. The PY theory is believed to underestimate the angle-dependent part of the PCFs while the HNC theory overestimates them. This is seen from the general stability condition for isotropic fluid which is derived from the Kerr constant [31] and is written as

1 ð0Þ 1  cˆ 22 V0: 5 It can be seen from Fig. 17 that both the PY and HNC values increase with q f*. However, the HNC results deviate rapidly from the low density linear behaviour and increase steeply in the vicinity of the phase transition. These steep increase can in fact be related to the growth of long-range orientational correlations. It is possible that the PY theory would eventually give a transition if numerical solutions could be obtained at higher values of q f*. It is likely, however, that such a transition would lie above the freezing density. The HNC theory predicts that the isotropic fluid of HER will be unstable for g z 0.435 which is close to the isotropic–nematic transition packing fraction [29]. For PY theory, it is found that the isotropic phase remains stable even at very high densities. A similar behaviour has been found by Perera and coworkers [32]. In the present paper we make another contribution to the theoretical description of understanding the equilibrium properties of a fluid of hard ellipsoids of revolution. Solutions of the HNC and PY integral equations are of considerable interest in the study of the structure and thermodynamic properties of fluids, since they are usually accurate and require much less computational effort than approaches based on Monte Carlo (MC)

9

and MD methods. Integral equations can be applied economically to a much broader range of problems than MC and MD methods. For fluids of axially symmetric nonspherical molecules, the angle-dependent PCFs are expanded in products of spherical harmonics either in a BF- or SF-coordinate frame. In a numerical calculation, one has to truncate the expansion series and has to solve the OZ equation with only a finite number of the harmonic coefficients. The convergence of the series depends on the values of x 0 and q f. We have examined the convergence of the expansion series and found that for x 0 V 3.0, one gets good results at all fluid densities if one considers all the terms of harmonic coefficients which involve l indices up to less than or equal to 6. Molecular fluids are much more difficult to treat theoretically and understanding the behaviour of anisotropic fluid systems is a continuing theoretical challenge.

Acknowledgements I am grateful to Y. Singh and J. Ram for many stimulating suggestions during the progress of this work and also for providing computational facilities. I thank the Director and Management of the Institute for providing the required facilities at MIET. The work was supported by the Department of Science and Technology (India) through a project grant.

References [1] J.A. Barker, D. Henderson, Rev. Mod. Phys. 48 (1976) 587. [2] J.P. Hansen, I.R. McDonald, Theory of Simple Fluids, Academic, London, 1976. [3] C.G. Gray, K.E. Gubbins, Theory of Molecular Fluids, Clarendon, Oxford, 1984. [4] P.H. Fries, G.N. Patey, J. Chem. Phys. 82 (1985) 429; P.H. Fries, G.N. Patey, J. Chem. Phys. 85 (1986) 7307; A. Perera, P.G. Kusalik, G.N. Patey, J. Chem. Phys. 87 (1987) 1295; A. Perera, P.G. Kusalik, G.N. Patey, J. Chem. Phys. 87 (1988) 5969E. [5] A. Perera, G.N. Patey, J. Chem. Phys. 89 (1988) 5861; A. Perera, G.N. Patey, J. Chem. Phys. 91 (1989) 3045. [6] J.M. Caillol, J.J. Weis, J.Chem. Phys. 90 (1989) 7403. [7] J.M. Caillol, J.J. Weis, G.N. Patey, Phys. Rev., A 38 (1988) 4772; J.M. Caillol, J.J. Weis, J. Chem. Phys. 92 (1990) 3197. [8] H. Zhong, R.G. Petschek, Phys. Rev., E 51 (1995) 2263; Phys. Rev., E 53 (1996) 4944. [9] R. Pynn, Solid State Commun. 14 (1974) 29. [10] R.C. Singh, J. Ram, Y. Singh, Phys. Rev., E 54 (1996) 977. [11] A. Perera, P.G. Kusalik, G.N. Patey, Mol. Phys. 60 (1987) 77. [12] J. Talbot, A. Perera, G.N. Patey, Mol. Phys. 70 (1990) 285. [13] R. Pospisil, A. Malijevsky, W.R. Smith, Mol. Phys. 79 (1993) 1011. [14] J. Ram, Y. Singh, Phys. Rev., A 44 (1991) 3718. [15] J. Ram, R.C. Singh, Y. Singh, Phys. Rev., E 49 (1994) 5117; S. Gupta, J. Ram, R.C. Singh, Physica, A 278 (2000) 447. [16] M. Letz, A. Latz, Phys. Rev., E 60 (1999) 5865. [17] A. Wulf, J. Chem. Phys. 67 (1997) 2254. [18] J.D. Parsons, Phys. Rev., A 19 (1979) 1125.

10

R.C. Singh / Journal of Molecular Liquids 122 (2005) 1–10

[19] S.D. Lee, J. Chem. Phys. 87 (1987) 4972. [20] M. Baus, J.L. Colot, X.G. Wu, H. Xu, Phys. Rev. Lett. 59 (1987) 2184. [21] J.F. Marko, Phys. Rev., A 39 (1989) 2050. [22] A. Chamoux, A. Perera, J. Chem. Phys. 104 (1995) 1493. [23] M.P. Allen, C.P. Mason, E. de Miguel, J. Stelzer, Phys. Rev., E 52 (1995) R25. [24] P.G. de Gennes, J. Prost, The Physics of Liquid Crystals, Oxford University Press, Oxford, 1995. [25] N.H. Phuong, G. Germano, F. Schmid, J. Chem. Phys. 115 (2001) 7227; N.H. Phuong, F. Schmid, J. Chem. Phys. 119 (2003) 1214;

[26] [27] [28] [29] [30] [31] [32]

N.H. Phuong, G. Germano, F. Schmid, Comput. Phys. Commun. 147 (2002) 350. M.P. Allen, G.T. Evans, D. Frenkel, B.M. Mulder, Adv. Chem. Phys. LXXXVI (1993) 1. R. Evans, in: D. Henderson (Ed.), Fundamentals of Inhomogeneous Fluids, Dekker, New York, 1992. Y. Singh, Phys. Rep. 207 (1991) 351. D. Frenkel, B. Mulder, J.P. McTague, Phys. Rev. Lett. 52 (1984) 287; B. Mulder, D. Frenkel, Mol. Phys. 55 (1985) 1193. B.J. Berne, P. Pechukas, J. Chem. Phys. 56 (1972) 4213. J.C. Filippini, Y. Poggi, J. Phys. Lett. 35 (1974) 99. A. Perera, G.N. Patey, J.J. Weis, J. Chem. Phys. 89 (1988) 6941.