Dynamical properties and transport coefficients of Kihara linear fluids

Dynamical properties and transport coefficients of Kihara linear fluids L. G. MacDowell, B. Garzo´n, S. Calero, and S. Lagoa) Departamento de Quı´mica Fı´sica I, Facultad de Ciencias Quı´micas, Universidad Complutense, 28040 Madrid, Spain

~Received 31 July 1996; accepted 10 December 1996! Transport properties of spherical and linear molecules modeled by the Kihara potential are studied by molecular dynamics simulations. Diffusion coefficients, shear viscosities, and thermal conductivities are calculated for a wide range of the fluid region and for several elongations. The corresponding individual and collective correlation functions are discussed along with angular velocity and reorientational correlation functions. Relaxation times and simple models relevant to orientational motion are also studied. The results obtained are discussed in a corresponding states framework, using previous Gibbs ensemble Monte Carlo data for the liquid–vapor equilibria of the models. In this way, the role of elongation can be studied. It is found that in most of the liquid region, the diffusion coefficient is weakly dependent on elongation. On the other hand, both viscosity and thermal conductivity are found to decrease with elongation. The dependence of transport coefficients on density and temperature is also discussed. On testing the Stokes–Einstein relation, it was observed that, unlike previous findings for hard spheres, stick boundary conditions perform just as good as slip boundary conditions for the Lennard-Jones fluid and the low-elongated Kihara fluid. © 1997 American Institute of Physics. @S0021-9606~97!50111-5# I. INTRODUCTION

The study of liquids by computer simulations has proved to be very successful in testing different theoretical approaches1–3 and in obtaining interesting correlations in the absence of a suitable theoretical framework.4 Contrary to conventional experiments, simulation experiments offer the possibility of studying the influence of an individual molecular feature, e.g., the dipole or quadrupole, simply by adding it to an otherwise completely unchanged molecular model. Considering such a singular advantage, we recently performed an exhaustive study of the influence of elongation and multipolar moments on the liquid–vapor equilibrium of simple molecules.4–6 Encouraged by the results, presented here is a corresponding study for transport properties, as there were now good reasons for attempting this problem. Although much progress has been done in the theoretical prediction of thermodynamic properties, both for simple7 and relatively complex fluids,8–10 the theory of transport processes presents great difficulties which still remain to be solved.7 In particular, the prediction of transport properties of simple real fluids in the dense regime remains unsolved, even though some progress has recently been reported.11–14 However, even those approaches using the pair potential as the only input, such as memory function formalism,14 rely on approximations whose significance is not totally understood. It is in this situation where simple correlations between the transport properties of different fluids are particularly useful. Probably the best way for searching for these correlations is by means of a corresponding states approach and such an attempt requires, first of all, a good knowledge of the phase diagram of the different substances to be compared. The development of the Gibbs ensemble Monte Carlo technique15 has permitted the determination of the vapor– a!

Author to whom correspondence should be addressed.

J. Chem. Phys. 106 (11), 15 March 1997

liquid equilibrium in a relatively fast and reliable way. Thus, the phase diagram of the Kihara model is now very well known.4–6 With this background, comparison of the transport properties of the different models in corresponding states warranted a study of their variation with elongation. Previous molecular dynamics studies of simple molecular fluids, using either equilibrium16,17 or nonequilibrium methods18 have been reported, but with minor attention to the effect of elongation, if any. When such was the case, attention was paid to individual transport properties and correlation functions,19 although the lack of knowledge of the phase diagrams could have rendered this problem difficult. Indeed, extensive simulations of the transport properties of hard ellipsoids with different length-to-breadth ratio are now available and very promising theoretical approaches have been developed for this model,20,21 but a systematic study of an anisotropic model with attractive forces on all of these properties had not been attempted previously. We concentrated on molecular models in which dispersive forces were present, so that our conclusions could be eventually applied to real substances of interest to the chemical industry. To imitate the interaction forces of the molecules, the Kihara model potential22 was used. Reorientational relaxation times, rotational and translational diffusion coefficients, shear viscosities and thermal conductivities of linear molecules of different elongations were calculated throughout most of the fluid regions of the models. Relaxation times and transport properties on the dense branch of the liquid–vapor coexistence curve, well inside the liquid pocket and in the supercritical region were evaluated. In order to fulfill the task, equilibrium molecular dynamics simulations were carried out. This method is particularly suitable, as the three transport properties can be evaluated simultaneously, together with the relaxation times. The paper is presented as follows: In Sec. II an account is given of the main equations needed to evaluate transport

0021-9606/97/106(11)/4753/15/$10.00

© 1997 American Institute of Physics

4753

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

MacDowell et al.: Dynamics of Kihara linear fluids

4754

properties. Section III deals with the simulation details, methodology, and model potential. Simulated thermodynamic states and the criteria used to compare transport properties of molecules with different elongation are also given here. Section IV gives an account of the dynamical properties of the potential studied and a simple model for reorientation is tested. The results obtained for the transport properties and collective modes are discussed in Sec. V. The final section provides the main conclusions drawn from our study.

II. THEORETICAL REMARKS

The different transport properties of a fluid can be expressed as integrals of the time correlation functions of the corresponding microscopic fluxes. In order to express the microscopic fluxes of molecular fluids, either an atomic or a molecular reference frame may be chosen.23 In the case of a nonflexible model such as the one used in this work, the fluxes given by Evans16,24 in the molecular reference frame seem the most appropriate. The diffusion coefficient, D, is evaluated using both the Green–Kubo and Einstein relations,7

E^ `

D51/3

0

v~ t ! •v~ 0 ! & dt,

~1a!

lim ^ @ Dr~ t !# & 56Dt,

~1b!

2

t→`

where r and v are the center of mass positions and velocities of the molecule, respectively, Dr(t)5r(t)2r~0! is the displacement of a particle from an arbitrary origin in time t and ^ & stands for a canonical average. The shear viscosity can be obtained from the correlation function of any of the off-diagonal elements of the symmetric part of the pressure tensor, through the following Green– Kubo equation:16

h5

V k BT

E^ `

0

xy J xy p ~ t ! •J p ~ 0 ! & dt,

~2!

where V and T are the volume and temperature of the system, k B is Boltzman’s constant and J xy p is the xy element of the ~symmetrized! pressure tensor. An expression for this tensor in terms of the microscopic variables of a rigid body has been given by Evans,16 1 VJp 5 m

N

(

i51

1 pi pi 2 2

N

N

]u

( ( ri j ] riijj ,

i51 jÞi

~3!

where ri j 5ri 2r j , pi is the momentum of particle i and u i j is the potential between i and j, assumed pairwise additive.22 For spherical particles with central forces, this tensor is necessarily symmetric. On the other hand, this is not generally the case for molecules interacting through noncentral potentials, so one must take care to symmetrize the Jp term. This is readily done by taking the average of the unsymmetrized pressure tensor and its transpose,16 so that

VJ xy p 5

1 m

N

( i51

1 1 2

p xi p iy 2

N

1 2

]u

ij r iyj x ( ]rij jÞi

N

( i51

D

S

1 2

N

]u

ij r xi j y ( ]rij jÞi

~4!

analogous expressions apply for the yz and zx components of the symmetrized pressure tensor. Note that the symmetrization procedure is not needed in the kinetic part. Actually, it turns out that on average, the pressure tensor of an isotropic fluid is always symmetric,25 so that here, this procedure is only used to gain some statistical accuracy and for the sake of consistency. For the same reason, three further components of the pressure tensor have been evaluated, i.e., those obtained by rotating 45° around each of the axis of the reference frame,26 – yy yy J xx 51/2~ J xx p p 2J p ! ,

J ypy – zz 51/2~ J ypy 2J zz p !,

~5!

– xx xx 51/2~ J zz J zz p p 2J p ! .

The thermal conductivity can be obtained from the correlation function of the heat flux vector, Jq ,7 l5

V 3k B T 2

E^ `

0

Jq ~ t ! •Jq ~ 0 ! & dt.

~6!

An expression for the heat flux in terms of rigid body molecular variables has also been given by Evans,24 N

VJq ~ t ! 5

( i51 2

1 2

S

p 2i

1 1 1 v i •I• v i 1 2m 2 2 N

N

( ( ri j i51 iÞ j

S

N

uij ( jÞi

D D

pi m

pi ] u i j • 2 v i •Ti j , m ] ri j

~7!

where vi is the angular velocity of particle i and Ti j is the torque that j exerts on i. This equation is also valid for atomic fluids by simply removing those terms related with the rotational motion. Usually, the properties of different fluids are most conveniently compared in the context of a law of corresponding states. In this respect, the Green–Kubo formula have proved to be very helpful, as their generality allows for such an approach in the particular case of transport properties.27 Thus, if one has a potential of the following general form: u i j ~ r! 5 e f ~ r/ s ! ,

~8!

27

then it can be shown that in the classical limit, the following reduced transport properties are universal functions of the reduced temperature T * 5k B T/ e and numerical density n * 5n s 3 , D * ~ T * ,n * ! 5D ~ T,n !~ m/ e s 2 ! 1/2,

~9!

h * ~ T * ,n * ! 5 h ~ T,n ! s 2 / ~ m e ! 1/2,

~10!

l * ~ T * ,n * ! 5l ~ T,n ! s 2 ~ m/ e k 2B ! 1/2.

~11!

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

MacDowell et al.: Dynamics of Kihara linear fluids

Deviations from this universal law appear as soon as the pair potential does not obey the form given in Eq. ~8!, as in our case. However, we found it convenient to keep the same reduced variables for comparison.

III. SIMULATION METHODOLOGY, MODEL POTENTIAL, AND SIMULATED THERMODYNAMIC STATES

The Kihara model fluid used in this work consists of infinitely thin rods of length L which interact according to the following pair potential: u i j ~ r ! 54 e @~ s / r ! 122 ~ s / r ! 6 # ,

~12!

where as usual, e and s are strength and size parameters and r is the shortest distance between rods. Note that r depends on positions and mutual orientations, so that in fact, r5r(r i j ,V i ,V j ), but we have suppressed the explicit dependence for the sake of clarity. The shortest distance between rods with a given orientation is evaluated using the algorithm of Ref. 28. Further details of this potential can be obtained from Refs. 4 and 6. Molecular dynamics ~MD! simulations were carried out in the microcanonical ensemble ~NVE! using the leapfrog algorithm.1 The rotational dynamics was implemented using the ‘‘orthogonal’’ constraint method proposed by Fincham.29 The forces required to run the MD algorithm are computed from the following equation:30 fi j 52

]uij m , ]r i j

~13!

where mi j is a unit vector in the direction of the shortest distance between two rods. It can be seen that in this potential, the force is exerted along the line joining the closest points between the molecular cores. To our knowledge, only one previous MD simulation has been carried out for the Kihara potential,30 and we refer to that paper for details. The moment of inertia of the rods is arbitrarily set equal to that of homonuclear diatomics with punctual masses. All the simulations were started from an a –N 2 lattice of 256 molecules and equilibrated for periods 5000 time steps long. The system was driven to the required temperature by means of the isokinetic thermostat,31 which was turned off 500 time steps before the beginning of the production run. The production length varied between 100 000 and 200 000 time steps, although it was usually 120 000 time steps long. The length of the time step was Dt *52•1023 reduced units @t * 5t( e /m s 2 ) 1/2#, except for those simulations with L *50 ~the Lennard-Jones potential!, where Dt * was set to either 3.95•1023 or 4.95•1023 reduced units. The Kihara potential was truncated at the maximum cut-off distance compatible with the minimum image convention, i.e., B/22L *, with B being the length of the simulation box. This gave a cutoff distance of 3.3 s in the worst case. Dynamical variables involved in the correlation functions were correlated every 20, 25 time steps during the simulation, so that collective correlation functions were averaged about 5000 times. A

4755

considerably greater statistical accuracy was obtained for the individual correlation functions, as they were averaged over the 256 particles. Transport coefficients were obtained by numerical integration of the corresponding correlation functions, according to the Green–Kubo equations of Sec. II. The diffusion coefficient was also evaluated using the Einstein relation, as it has been shown that it provides a better convergence than its Green–Kubo analog.32 However, Einstein relations for the collective properties present additional problems that make them unsuitable for practical purposes.33 The uncertainty in the diffusion coefficient has been obtained from the prescription made by Zwanzig and Ailawadi.34 As to the viscosity coefficient, all six elements of the pressure tensor were used in its evaluation to improve statistical accuracy. Thus, the error in h was estimated from the root mean square deviation of these six elements, considered independently. The same procedure was used to evaluate the uncertainty in the thermal conductivity. In this work, attention is paid to linear apolar Kihara molecules, whose phase diagram is known from Ref. 5. In that work, it was shown that the liquid–vapor equilibria of molecules with different elongation presented a slight departure from the principle of corresponding states, so that the coexistence curves were broadened as elongation increased. This actually means that corresponding states can no longer be defined through the usual reduction with the critical density and temperature, since this would lead to compare transport properties of spherical molecules on the coexistence curve, with those of elongated molecules on the metastable region of the phase diagram. Rather, thermodynamic states on the coexistence curve with equal T to T c ratio are compared. With this idea in mind, the following thermodynamic states were simulated for molecules with L *50.0 ~the Lennard-Jones fluid!, 0.3, 0.6, and 0.8. ~1! States on the liquid–vapor coexistence curve. Three different states for each elongation were studied in this region. One state near the critical point, with T/T c 50.90; a second one with T/T c 50.65, and finally one near the triple point, such that T/T c 50.55. These points were named A1, B2, and T3, respectively. Several other thermodynamic states were then simulated using these three initial states as a reference. ~2! States inside the liquid pocket. In this region three different states were studied. Two along the T/T c 50.90 isotherm, with densities equal to those of points B2 and T3, and one along the T/T c 50.65 isotherm, with a density equal to that of point T3. These states were named A2, A3, and B3, respectively. ~3! States in the gas region. Four other states were studied, with T/T c 51.5. One of them in the semidilute regime, with n/n c 50.60 and three others with densities equal to those of points A1, B2, and T3. These states were named C0, C1, C2, and C3. The ten states studied are represented on the phase dia-

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

MacDowell et al.: Dynamics of Kihara linear fluids

4756

TABLE II. Simulated thermodynamic states as defined by their temperature and density in units of the corresponding critical densities and temperatures. For a schematic representation of these points, see Fig. 1. L *50.0 States .T/T c A1 B2 T3 A2 A3 B3 C0 C1 C2 C3

FIG. 1. Schematic representation of the thermodynamic states studied on the phase diagram of the Lennard-Jones fluid. For the elongated models, this representation is only roughly approximate, since the different phase diagrams do not follow exactly the principle of corresponding states. Note that each point is labeled so that equal letters represent states along the same isotherm, while equal numbers represent states along the same isochore.

gram of a Lennard-Jones fluid in Fig. 1. The phase diagram of the other models will look somewhat different, but the figure still shows roughly the relative location of the points studied. Summarizing, each state is represented by a letter indicating the T/T c ratio, followed by a number which indicates the n/n c ratio. Thus, points A1, A2, and A3 lay along the same isotherm, while points C3, A3, B3, and T3 lay along the same isochore. In Table I, the critical parameters of the different Kihara models as determined in Ref. 5 and the Lennard-Jones critical parameters as determined from Ref. 35 are shown. As the NVE ensemble was used, fluctuations of the order of 1/N 1/2 with respect to the desired temperature are found.1 The thermodynamic states actually simulated are given in reduced units in Table II. IV. DYNAMICAL PROPERTIES

A glance at the thermodynamical properties, obtained as a by-product of this work, may give us some hints as to the dependence of transport properties with elongation. In Table III the results for the reduced configurational energies and pressures are shown. Both configurational energies and pressures are seen to decrease gradually with elongation, when compared within the correspondence criteria of the preceding TABLE I. Critical parameters of the different Kihara models studied.a T c* 5 kT c / e , n c* 5 n c s 3 , and P c* 5 P c s 3 / e are the reduced critical temperature, 3 density, and pressure. y c 5 ( p /6)(1 1 2L * )n c* is the critical packing fraction.

a

L*

T c*

n c*

P c*

yc

0.0 0.3 0.6 0.8

1.310 1.114 1.000 0.952

0.314 0.219 0.161 0.140

0.130 0.073 0.051 0.038

0.162 0.166 0.160 0.161

References 5 and 35.

0.90 0.65 0.55 0.90 0.90 0.65 1.50 1.50 1.50 1.50

L *50.3

L *50.6

L *50.8

T/T c

n/n c

T/T c

n/n c

T/T c

n/n c

T/T c

n/n c

0.8840; 0.6356; 0.5508; 0.8884; 0.8597; 0.6563; 1.483; 1.523; 1.515; 1.496;

1.78 2.48 2.69 2.48 2.69 2.69 0.60 1.78 2.48 2.69

0.9200; 0.6719; 0.5580; 0.8900; 0.8587; 0.6613; 1.544; 1.502; 1.480; 1.410;

1.90 2.50 2.68 2.50 2.68 2.68 0.60 1.90 2.50 2.68

0.9128; 0.6545; 0.5390; 0.8809; 0.8783; 0.6383; 1.502; 1.527; 1.458; 1.501;

1.96 2.63 2.82 2.63 2.82 2.82 0.60 1.96 2.63 2.82

0.896; 0.637; 0.542; 0.876; 0.889; 0.648; 1.51; 1.484; 1.446; 1.500;

1.94 2.61 2.82 2.61 2.82 2.82 0.60 1.94 2.61 2.82

section. The decrease of configurational energy means that as elongation is increased, the system takes on a more repulsive character. It is therefore expected that diffusion will increase30 and viscosity will decrease as elongation is increased. This latter statement is supported by the fact that the trace of the pressure tensor, i.e., the pressure itself, also decreases. We shall now analyze the different individual correlation functions computed during our simulations.

A. Translational dynamics

A simple way of describing the translational dynamics of a fluid is by means of the velocity autocorrelation function, defined by the following expression:

w v ~ t ! 5 ^ v~ t ! •v~ 0 ! & ,

~14!

where ^ & means ensemble averaging. In Fig. 2 the normalized velocity correlation function of states C3 and T3 is shown for the four elongations studied. Both states C3 and T3 show roughly a similar behavior. There is an initial sudden decay, followed by a much smoother one, and a minimum is clearly visible. As the elongation decreases, this minimum is more evident and becomes negative for the less elongated models. For the longest molecule ~L *50.8!, the development of an inflexion point which becomes a negative minimum for the T3 state can be seen. In this case, the formation of a double negative minimum is apparent. The nonexponential decay of wv (t) is a consequence of correlated collisions where the structure of the system plays a major role. In particular, it is well known that the negative minimum of spherical fluids at high densities is related to caging effects followed by backscattering of the molecules.2,36 Indeed, it has been shown that a very simple structural model in which particles diffuse along tunnels made up of their own neighbors can account for this negative minimum.37 On the other hand, the double minimum structure for the more elongated molecules shows evidence for the occurrence of two different correlated processes.

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

MacDowell et al.: Dynamics of Kihara linear fluids

4757

TABLE III. Reduced configurational energy, U * 5U/ e , and pressure P * 5 P s 3 / e as obtained from the MD calculations of this work. L *50.0

L *50.3

L *50.6

L *50.8

States

.T/T c

2U *

P*

2U *

P*

2U *

P*

2U *

P*

A1 B2 T3 A2 A3 B3 C0 C1 C2 C3

0.90 0.65 0.55 0.90 0.90 0.65 1.50 1.50 1.50 1.50

3.860; 5.555; 6.086; 5.276; 5.663; 5.933; 1.239; 3.510; 4.683; 4.932;

20.007 20.051 0.141 1.456 2.549 0.994 0.301 1.409 4.644 6.566

3.693; 5.104; 5.601; 4.893; 5.242; 5.473; 1.099; 3.394; 4.391; 4.677;

0.065 0.084 0.044 0.996 1.662 0.606 0.185 1.046 3.133 4.151

3.535; 5.005; 5.500; 4.782; 5.106; 5.378; 1.010; 3.227; 4.314; 4.491;

0.027 0.024 20.038 0.748 1.374 0.399 0.122 0.774 2.433 3.493

3.485; 4.977; 5.496; 4.751; 5.093; 5.364; 0.993; 3.188; 4.290; 4.489;

0.004 20.028 0.010 0.645 1.304 0.422 0.103 0.626 2.086 3.182

In order to illustrate this point, it is common to split the velocity autocorrelation function into two components, wvi (t) and w'v (t), so as to study separately the velocity parallel and perpendicular to a unit vector, e, along the molecular axis,38

w v ~ t ! 5 ^ @ v~ t ! •e~ 0 !#@ v~ 0 ! •e~ 0 !# & , i

~15a!

w'v ~ t ! 51/2@ w v ~ t ! 2 w v ~ t !# . i

~15b!

w 'v (t)

are shown for In Fig. 3, normalized wv (t), wv (t), and state T3 and elongations L *50.3 and L *50.8. The figure clearly shows the differences in the motion of the parallel and perpendicular modes of diffusion, with the latter decaying much faster at the short time. This figure seems to indicate that the decreasing depth of the minimum as elongation increases is not only due to the smaller numerical density of the elongated molecules,19 but also to a relative enhancement of the motion parallel to the molecular axis. Indeed, Fig. 3 shows how the negative minimum disappears in fv while remaining in the perpendicular component going from L *50.3 to L *50.8. Further support for the latter hypothesis comes from the observation that at the onset of nematic phases in linear molecules, the parallel component of diffusion may even increase with density.39,40 In fact, for a highly idealized model of infinitely thin hard needles, it has been shown that the diffusion parallel to the molecular axis diverges as density increases.41 A striking feature of Fig. 2, where the normalized velocity autocorrelation functions are shown, is that the initial velocity decay up to almost 0.1 time units is almost identical for all the elongations ~this was also observed in most of the other thermodynamic states, not reported here!. The short time decay of the velocity correlation function is given by the mean squared force per particle, ^F2&, so that7 i

F

w v ~ t * ! 5 ^ v 2 & 12

G

^ F* 2 & t * 2 1O ~ t * 4 ! . ^ v* 2 & 2!

~16!

Since the normalized velocity correlation function, w v (t * )/^v2&, decays similarly for all the elongations, it was expected that the ratio ^F*2&/^v*2& should also be quite similar. Further calculations showed that this expansion was capable of predicting w v (t * ) for up to 0.2 reduced time units, so that the ratios ^F*2&/^v*2& are not identical, but close

FIG. 2. ~a! Normalized velocity autocorrelation function for the C3 state and the four elongations. ~b! Same as ~a!, but for the T3 state.

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

4758

MacDowell et al.: Dynamics of Kihara linear fluids

FIG. 3. ~a! Normalized velocity autocorrelation function for state point T3 and L *50.3, resolved into parallel and perpendicular components to the molecular axis. ~b! Same as ~a!, but for L *50.8.

enough to show differences in the velocity decay only after 0.1 reduced time units. For example, for the T3 states, the following ^F*2&/^v*2& were obtained: 274 ~L *50!, 240 ~L *50.3!, 228 ~L *50.6!, 229 ~L *50.8!, and the differences are probably within the simulation error except for L *50. Observation of w v (t) for the other thermodynamic states not reported here support our conclusions. Furthermore, a clear exponential behavior of w v (t) for the C0 state and all the elongations was observed, as expected. Similar behavior was found for states A1 and C1, while all the others shared similar characteristics as those found in Fig. 2.

FIG. 4. ~a! First order reorientational function for state A1 and L *50.3, L *50.6, and L *50.8. ~b! Same as ~a! but for the C3 state. ~c! Same as ~a! but for the T3 state.

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

MacDowell et al.: Dynamics of Kihara linear fluids

4759

B. Rotational dynamics

The rotational dynamics of the system has been monitored by means of the angular velocity correlation function, wv(t) and the first and second order reorientational correlation functions, w l51 (t) and w l52 (t),

w v~ t ! 5 ^ v ~ t ! • v ~ 0 ! & ,

~17!

w l51 ~ t ! 5 ^ P l51 @ e~ t ! •e~ 0 !# & 5 ^ e~ t ! •e~ 0 ! & 5 ^ cos Q ~ t ! & ,

~18a!

w l52 ~ t ! 5 ^ P l52 @ e~ t ! •e~ 0 !# & 5 ^ 3 @ cos Q 2 ~ t ! 21 # & ,

~18b!

where P l (x) is the lth order Legendre polynomial and Q(t) is the angle formed between e~0! and e(t). w l51 (t) is shown in Fig. 4 for states A1, C3, and T3 and the three elongations different from zero, L *50.3, 0.6, and 0.8. This figure clearly illustrates how molecules gradually lose their freedom to rotate as elongation and density are increased and temperature is decreased. On doing so, the negative minimum of w l51 (t) characteristic of the inertial regime is gradually lost. @Note the two inflexion points for L *50.3 in Fig. 4~b! reminiscent of this regime.# For the T3 state, rotation of the molecules is strongly hindered, so that it takes a long time before they can escape from the cage surrounding them and lose the memory of their initial orientation. It is expected that the diffusion model of Debye should work on these states, where relaxation must be the consequence of small amplitude reorientations. Although the simplicity of Debye’s diffusion model makes it very appealing, it is well known that it fails to predict correctly the short time reorientational behavior. The search of simple models capable of describing both the inertial short time behavior and the long term diffusion is thus of great interest. We have tested the Gaussian cage model of Lynden-Bell and Steele,42 which is appropriate for dense liquids with relatively high torques. These authors applied a cumulant expansion to w l (t) up to second order and found the following equation: ln w l ~ t ! 52l ~ l11 ! /2

E

1

0

~ t2 t ! w v ~ t ! d t .

~19!

They proceeded by assuming that the molecules in the fluid performed small amplitude librations with an average frequency V and a spread given by D. This assumption results in a simple expression for wv(t),

w v~ t ! 5

2k B T cos~ Vt ! exp~ 2D 2 t 2 /2! . I

~20!

On replacing Eq. ~20! in Eq. ~19!, one gets a simple expression for w l (t) as a function of V and D, which are taken as adjustable parameters. We have tested this model on state point T3. The result fits to both w l51 (t), w l52 (t) and the reduced wv(t) shown in Fig. 5. The model does a good job in fitting w l51 (t) and w l52 (t) up to about t *52 for all the elongations. On the contrary, it is uncapable of describing the ‘‘fanning’’ out of

FIG. 5. Test of the Gaussian cage model for state T3 and L *50.3, L *50.6, and L *50.8. ~a! Logarithmic plot of the first order reorientational function. The symbols are molecular dynamics results, while the lines are a fit of the Gaussian cage model to the data. ~b! As in ~a! but for the second order reorientational function. ~c! Fit of the model to the normalized angular velocity correlation function. Symbols, molecular dynamics results, line, best fit.

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

MacDowell et al.: Dynamics of Kihara linear fluids

4760

TABLE IV. Test of the Gaussian cage model for reorientation at state point T3. The symbols are defined as in the text.

^T2&* 2

L* 0.3

0.6

0.8

^T2&*

l

V*

D*

T* ~fit!

T*2 ~MD!

1 2 wv 1 2 wv 1 2 wv

8.59 10.53 0.004 8.16 8.69 0.006 8.02 8.18 0.001

5.22 4.47 11.8 5.80 5.25 12.93 5.52 5.02 13.05

12.2 14.7 15.3 37.8 38.8 60.5 62.2 60.5 108.5

15.9

TABLE V. Reduced relaxation times for reorientational motion, t 1* , t 2* , and rotational diffusion, t v* , at several state points and elongations. Test of the first and second order Hubbard relations are also shown, together with the ratio of one to the other. These relaxation times were measured at slightly different temperatures from those of Table II. We therefore show the corresponding temperatures in the third column. State L *

69.8

C0

125.9

C2

C3

A1

w l52 (t) at long times, which clearly requires higher order terms in the cumulant expansion. It is also apparent that the fit improves as elongation increases. Figure 5~c! shows that wv(t) is fairly well described for L *50.6 and L *50.8, whose fit is in fact almost identical for both elongations, while it fails to describe it for L *50.3. This, in fact, is not surprising, since the average torque on a molecule greatly diminishes as elongation is decreased. A look at Table IV, where the parameters of the fits are shown, provides a more severe test of this model. It is seen that the parameters are quite similar for the fit to w l51 (t) and w l52 (t) but they are dramatically different for wv(t). This shows that even if wv(t) can be fit reasonably by Eq. ~20!, the high torque regime in which Eq. ~19! is valid is not reached. It can be shown that the mean square torque on a molecule is related to the parameters of the model, V and D.42 The resulting estimates are also shown in Table IV, together with those obtained from the simulations. As has been previously observed,42 the mean square torque is systematically underestimated, even though the fit to wv(t) gives better results. The corresponding relaxation times of the correlation functions studied are also important properties which allow for a test of different reorientational models. In Table V we show tv , t1 , and t2 , defined by the following expressions:

t v 51/^ v 2 & t l5

Ew `

0

Ew `

0

l ~ t ! dt.

v ~ t ! dt,

~21! ~22!

As expected, for a given state point tv always decreases as elongation is increased, while the contrary is true for t1 and t2 . Other less obvious conclusions may be obtained from this table. It is observed that for a given elongation, an increase in temperature is followed by a decrease in the rotational diffusion, tv , while the orientational relaxation, given by t1 and t2 , decreases. The only way to reconcile this apparently paradoxical observation is that the decrease in tv is compensated by the increase in the average angular speed, thus allowing for a faster relaxation of orientations.

A2

A3

B2

B3

T3

0.3 0.6 0.8 0.3 0.6 0.8 0.3 0.6 0.8 0.3 0.6 0.8 0.3 0.6 0.8 0.3 0.6 0.8 0.3 0.6 0.8 0.3 0.6 0.8 0.3 0.6 0.8

T/T c 1.491 1.492 1.52 1.541 1.461 1.47 1.511 1.373 1.37 0.9003 0.8769 0.884 0.9031 0.9064 0.907 0.8846 0.9025 0.876 0.6405 0.6640 0.650 0.6612 0.6588 0.660 0.5507 0.5330 0.548

t v* •10 t *1 •10 t *2 •10 I/k B T t v t 1 I/k B T t v t 2 t1/t2 12.5 10.5 10.2 1.39 1.02 1.03 1.12 0.86 0.83 2.76 2.27 2.34 1.45 1.11 1.14 1.16 0.90 0.87 1.50 1.16 1.18 1.18 0.91 0.86 1.19 0.93 0.87

0.46 1.20 1.77 1.12 3.78 6.26 1.25 4.08 7.93 1.15 3.48 5.46 1.60 5.34 9.23 1.82 6.53 12.14 2.04 7.02 12.28 2.34 8.79 15.96 2.63 10.30 18.92

1.84 2.25 2.62 0.64 1.60 2.50 0.61 1.03 2.74 0.94 1.92 2.78 0.91 2.39 3.69 0.98 2.67 4.69 1.16 3.01 4.89 1.22 3.55 6.02 1.41 4.29 7.43

0.23 0.48 0.57 0.8 1.6 1.8 0.1 1.9 1.9 0.7 1.3 1.5 1.0 1.7 1.8 1.1 1.7 1.8 1.0 1.7 1.8 1.1 1.7 1.8 1.2 1.8 1.8

0.1 0.3 0.4 1.5 3.8 4.5 2.0 7.7 5.6 0.9 2.4 2.9 1.7 3.7 4.3 2.0 4.2 4.8 1.8 3.8 4.5 2.1 4.2 5.0 2.2 4.2 4.8

0.2 0.5 0.7 1.8 2.4 2.5 2.0 4.0 2.9 1.2 1.8 2.0 1.8 2.2 2.5 1.9 2.4 2.6 1.8 2.3 2.5 1.9 2.5 2.6 1.9 2.4 2.6

The Hubbard relations allow us to determine whether the reorientation in a liquid is proceeded by small angular displacements of the molecule or not.43 These relations read as follows:43 I 5l ~ l11 ! . k BT t vt l

~23!

If l51, then Hubbard’s relation predicts a ratio of I/k B T t v t l 53, while if l52, I/k B T t v t l 56. As expected from Debye’s diffusion, it also predicts that t1/t253. In Table V we test Hubbard’s relations for several states. As expected, the regime where this law is obeyed is approached sooner as elongation or density are increased. It is also observed that at high densities, the rotational regime is only slightly dependent on temperature, since the calculated ratios remain almost unchanged. V. TRANSPORT COEFFICIENTS AND COLLECTIVE CORRELATION FUNCTIONS

In this section, transport properties as obtained from the molecular dynamics calculations of this work are discussed. A comparison of our results for the Lennard-Jones system with those of other authors is shown in Table VI. The agreement is satisfactory for all the transport properties. It should be noted that an exact agreement is not expected since it is

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

MacDowell et al.: Dynamics of Kihara linear fluids TABLE VI. Results of the transport properties of the Lennard-Jones compared with other work. Note that the compared results do not always belong to exactly the same state point. All the results were obtained with systems of 256 particles but for those of footnotes a, c, and e, in which it is claimed that the thermodynamic limit was reached. States

D *•102

h*

l*

T3

3.1 3.01b 3.35d

3.3 3.19a 2.98b 3.01d 3.34e 3.30g 2.0 1.96d 1.86f 0.67 0.66d 2.8 2.65f 1.8 1.64f

6.8 6.96d 6.76c

B2

A1 B3 A2

5.7 5.66d 5.89f 21.1 17.7c 4.0 4.02f 8.6 8.35f

5.5 5.59d 5.46f 2.6 2.56d 7.3 7.14f 5.9 5.99f

a

Reference 54. b Reference 53. c Reference 58. d Reference 59. e Reference 44. f Reference 52. g Reference 57.

not always possible to compare with data at exactly the same thermodynamic states. Furthermore, it has been suggested that transport properties, or at least the shear viscosity, may be dependent on the cutoff distance of the potential.44 Whereas this distance has been traditionally taken as 2.5s, in this work it was set to half the box length, so that it was at least 3.3s. This means that our results must suffer from a smaller systematic error than is usually the case. Note also that the value of e needed to describe a real fluid with the Kihara model is in general much bigger than that with a two-centers Lennard-Jones,45 so that a direct comparison in reduced units cannot be done.

4761

TABLE VII. ~a! Reduced diffusion coefficient, D * 5D(m/ e s 2 ) 1/2, as a function of L * for the different corresponding states of this work. Subcritical temperatures. Numbers in parentheses give a measure of the uncertainty in the last decimal digit. ~b! As in ~a! but for the supercritical temperatures. ~a! States Elongation

A1 D*

A2 D *•102

A3 D *•102

B2 D *•102

B3 D *•102

T3 D *•102

L *50.0 L *50.3 L *50.6 L *50.8

0.211~8! 0.192~5! 0.191~4! 0.189~4!

8.6~1! 8.2~4! 7.8~3! 7.9~3!

5.7~1! 5.7~4! 5.9~3! 6.1~3!

5.7~4! 5.9~3! 5.7~3! 5.4~2!

4.0~4! 4.0~2! 4.0~2! 4.0~2!

3.1~4! 3.2~2! 3.2~2! 3.2~2!

~b! States Elongation L *50.0 L *50.3 L *50.6 L *50.8

C0 D*

C1 D*

C2 D*

C3 D*

1.31~1! 1.40~2! 1.47~1! 1.46~1!

0.322~7! 0.285~6! 0.298~6! 0.283~6!

0.15~1! 0.137~6! 0.133~6! 0.135~8!

0.114~5! 0.107~6! 0.108~6! 0.107~2!

L * 50.6 D * ~ T * ! 50.079 007T * 20.010 481,

~24c!

L * 50.8 D * ~ T * ! 50.082 282T * 20.009 843.

~24d!

It is far more illuminating to fit D * for all the elongations as a function of T/T c . The resulting linear function along the C3 – A3 – B3 – T3 isochore, together with the D *s, is shown in Fig. 6. It is seen that D * can be predicted quite well in a corresponding states manner by the following linear function: all elongations

D * ~ T/T c ! 50.083 123T * 20.013 889. ~25!

This supports what is already apparent from direct observation of Tables VII~a! and VII~b!, i.e., that D * is quite independent of elongation for the highest densities and lower temperatures, at least for the elongations considered in this

A. Self-diffusion

The reduced diffusion coefficients for the different states and elongations are shown in Tables VII~a! and VII~b! as obtained via the Einstein relation. The use of the Green– Kubo relation always gave the same result within the experimental error except for the C0 states, in which w v (t) had not decayed to zero as far as 3 reduced time units. In the dense fluid, D is seen to increase with temperature and decrease with increasing density. The temperature dependence at constant density in a range as long as 0.95 times T c , is found to be linear for all the elongations. The same sort of dependence has been found by other authors for the Lennard-Jones46 and the two centers Lennard-Jones.19 Thus, for the C3 – A3 – B3 – T3 isochore, the following fits were obtained: L * 50.0 D * ~ T * ! 50.067 149T * 20.017 856,

~24a!

L * 50.3 D * ~ T * ! 50.079 438T * 20.018 190,

~24b!

FIG. 6. Plot of the reduced diffusion coefficient against T/T c along a given isochore for the four elongations. ~Top! B2 – A2 – C2 isochore. ~Bottom! T3 – B3 – A3 – C3 isochore.

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

MacDowell et al.: Dynamics of Kihara linear fluids

4762

TABLE VIII. Reduced diffusion coefficients as obtained from the Stokes–Einstein relationship. l stands for slip and t for stick boundary conditions. States

C0

C2

A1

S

T3

Elongation

D l*

D t*

D l*

D t*

D l*

D* t

2 D* l •10

2 D* t •10

2 D* l •10

2 D* t •10

L *50.0 L *50.3 L *50.6 L *50.8

1.0 0.97 1.1 1.5

0.68 0.64 0.72 0.98

0.20 0.15 0.17 0.13

0.13 0.10 0.11 0.086

0.28 0.25 0.24 0.28

0.18 0.17 0.15 0.14

6.6 7.0 7.0 6.7

4.4 4.7 4.7 4.5

3.5 4.2 3.8 3.5

2.3 2.8 2.6 2.3

work. We note, however, that if D * is fit as a function of T/T c for each elongation separately, the linear coefficients are 0.088, 0.089, 0.079, and 0.078 from the lowest to the highest elongation. On the other hand, for the C2, A2, and B2 states, with liquidlike density, D * seems to decrease with L * ~although a corresponding states fit for the C2 – A2 – B2 isochore still gives reasonable results, as shown in Fig. 6!. This is contrary to the expectations made from considerations of the dependence of configuration energy with elongation. Indeed, in Sec. IV it was conjectured that the decrease of U * and thus the increasing repulsive nature of the systems would cause an enhancement of D *. It must be concluded that this effect is canceled by the hindrance of translation, due to end to end collisions. Observation of w v (t) for state C0 ~not shown here!, shows an almost exponential decay, so that one might expect that a one mode modified Enskog theory47 can predict the diffusion coefficient reasonably well. On fitting these functions to a single exponential, the following inverse reduced relaxation times were found: 1.48 for L *50; 1.20 for L *50.03; 1.03 for L *50.6, and 0.99 for L *50.8. These should be compared fairly well with the predictions of Enskog’s theory,7 i.e., D * /T *, which give 1.48, 1.19, 1.03, and 0.99, respectively. On the other hand, the self-diffusion of a tagged spherical particle of radius R can be given in terms of the viscosity of the medium, h, by the following hydrodynamical equation:48 D5

B2

D

113 h / b k BT , 6 p h R 112 h / b

~26!

where b is a parameter which measures the friction between the tagged particle and the medium. If b is very small, this equation is equivalent to the Stokes–Einstein relation for slip boundary conditions, while if b goes to infinity, the stick boundary conditions are recovered. Although this equation is soundly founded only for large particles in a continuous medium, it can be shown that it can also apply to molecular processes by introducing reasonable microscopic approximations.48 In order to apply this equation to nonspherical particles, a geometrical correction factor should be included. However, for these elongations it is expected that this factor be very small, so that an effective radius can be simply obtained by equating the volume of the spherocylinder to that of a sphere. In Table VIII, the Stokes–Einstein

equation is tested for a few states using the viscosity coefficients obtained in this work ~see next section!. Comparing with the simulation results of Tables VII~a! and VII~b!, it is seen that the Stokes–Einstein equation shows roughly the dependence of D * with elongation described before. Furthermore, stick boundary conditions are seen to predict the diffusion coefficient just as well as slip boundary conditions, although the latter seem to improve as elongation is increased. However, previous evidence has usually supported the idea that the slip boundary conditions are more appropriate to describe microscopic behavior.2,19,49–51 This idea probably traces back to the pioneering work of Alder et al., on hard spheres.2 Actually, on investigating this problem further, we observed that there exist many molecular dynamics calculations giving better diffusion coefficients using stick rather than slip boundary conditions, at least in the liquid density region. We refer to results by Borgelt et al.52 and Heyes,53 which are probably the two most extensive calculations of diffusion coefficients and viscosities for the Lennard-Jones fluid. In order to assess which of the boundary conditions gives better results, the diffusion coefficients given by the Stokes–Einstein relationship were computed and the mean square relative deviation from the molecular dynamics data was evaluated ~i.e., we calculate the mean of @~D stick/slip/D MD21!#2. Borgelt et al. give shear viscosities and diffusion coefficients for three reduced densities, 0.781, 0.842, and 0.884 and reduced temperatures ranging roughly from 0.66 to 2.7. For each of the three isochores, the following deviations are obtained when using stick boundary conditions, 0.025, 0.022, 0.018, respectively. On the other hand, slip boundary conditions give the following deviations, 0.092, 0.094, 0.15. It is clear that in this case, the deviations obtained using stick boundary conditions are much smaller than when using slip boundary conditions. The work of Heyes gives transport coefficients along most of the phase diagram of the Lennard-Jones, but those states with reduced viscosities less than unity are not studied here, since it has been argued that the nonequilibrium method used by Heyes is somewhat unreliable when dealing with low viscosity states.54 Those state points whose viscosity appears in parentheses53 are not used either. The deviations obtained using the 34 remaining states are 0.046 and 0.11 for stick and slip boundary conditions, respectively. Again, the predictions given by the stick boundary conditions are much better, despite the fact that in his paper,

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

MacDowell et al.: Dynamics of Kihara linear fluids TABLE IX. ~a! Reduced shear viscosity, h*5hs2/(m e ) 1/2, as a function of L * for the different corresponding states of this work. Subcritical temperatures. Numbers in parentheses give a measure of the uncertainty in the last decimal digit. ~b! As in ~a! but for supercritical temperatures. ~a! States Elongation L *50.0 L *50.3 L *50.6 L *50.8 ~b! States Elongation L *50.0 L *50.3 L *50.6 L *50.8

A1 h*

A2 h*

A3 h*

B2 h*

B3 h*

T3 h*

0.67~4! 0.58~2! 0.49~1! 0.49~2!

1.8~2! 1.3~2! 1.2~1! 1.3~3!

2.4~1! 2.0~4! 1.6~3! 2.0~2!

2.0~1! 1.5~3! 1.2~2! 1.1~2!

2.8~1! 2.2~2! 1.8~2! 1.6~2!

3.3~2! 2.1~2! 1.8~4! 1.8~2!

C0 h*

C1 h*

C2 h*

C3 h*

0.30~3! 0.25~3! 0.18~4! 0.12~2!

0.77~3! 0.66~3! 0.56~5! 0.41~2!

1.6~2! 1.5~3! 1.1~2! 1.3~5!

2.2~2! 1.7~3! 1.5~3! 1.7~4!

Heyes uses the slip boundary conditions in order to evaluate the hydrodynamic radius of the different states. It was noted in the results by Heyes, however, that for reduced densities close to one, the slip boundary conditions gave almost quantitative predictions for the diffusion coefficient. Indeed, Eyring48 showed a long time ago that this is to be expected when the local structure of the liquid is almost solid like. B. Shear viscosity and the stress correlation function

The results obtained for the shear viscosity are given in Tables IX~a! and IX~b!. Inspection of these tables shows that the viscosity decreases monotonously as elongation is increased, as expected from considerations on thermodynamic properties, put through in Sec. IV. This decrease is particularly pronounced on going from L *50.0 to L *50.3. For the liquid density region, viscosity is seen to increase steeply with density, while the temperature dependence is not very clear, although it seems to show a weak decrease as temperature is raised. According to a law by de Guzman Carrancio,55 viscosity follows an Arrhenius sort of behavior, with a negative viscous flow activation energy. The failure to show a clear trend in the temperature dependence illustrates the difficulties encountered in the calculation of shear viscosity by molecular dynamics. The difficulties arise from two different sources. First of all, the collective nature of this property, which makes it necessary to perform very long runs in order to gain the same statistical accuracy as that of individual properties such as selfdiffusion. Secondly and most important, the relatively long time tails which appear at high densities. These long time tails are very hard to reproduce in independent runs and give rise to most of the statistical uncertainty in the values of h.56 Furthermore, the impossibility of integrating to infinity introduces a systematic error which is impossible to evaluate. In this work, it was found necessary to integrate the stress autocorrelation function up to rather long times, in an attempt to consider the difference between the long time tails of molecules with different elongation. For this reason, the trans-

4763

port coefficients are given with a somewhat greater statistical error than is usually the case. In Fig. 7 the normalized stress autocorrelation functions of states C0, C3, and T3 are shown for all the elongations. Figures 7~b! and 7~c! show two of the features which are usually seen in the stress correlation functions of rigid anisotropic molecules,17 i.e., a very sharp decay which is steeper for the longest molecules and a long time tail which gives the stress autocorrelation function its peculiar ‘‘chair’’ form. The long time tail becomes more important as elongation increases. It is also apparent from these figures that the value of the tail is so small for reduced time units of about 0.5 that it becomes very hard to distinguish it from statistical noise. It is usually said that low density states present no problem in the evaluation of h, since they lack a long time tail. However, in Fig. 7~a! it can be seen that although this seems to be true, the decay is so slow that they must be evaluated for times as long as 3 reduced units. A very interesting feature observed in the three figures is the gradual development of a minimum followed by a maximum as the thermodynamic states become more condensed and elongation increases @see Figs. 7~b! and 7~c!#. This behavior is similar to that shown by the velocity correlation function due to the appearance of two different modes of translation. Obviously, this comparison cannot be taken too far, since the latter are individual, while the former are collective properties. In a recent paper, an Enskog-type theory for the viscosity of hard ellipsoids was developed where for the first time an account was made of the two collective modes that appear as a consequence of the anisotropy of linear molecules.21 The results of this theory were somehow disappointing because it predicted viscosities very similar to those of a single mode Enskog theory.20 However, Fig. 7~a! shows that even at densities where the Enskog theory should give results in reasonable agreement with simulation, there seems to be a peculiar behavior characteristic of the elongated molecules, which should be related to the anisotropic interactions. It is also noted that the short time decay of the stress correlation function is very similar for the three linear models studied, but differs much more for the Lennard-Jones. As was the case for wv ~1!, the initial decay of the stress correlation function is also determined by a well defined equilibrium property, namely the mean squared momentum flux time derivative. Apparently, the greater number of degrees of freedom in the linear molecules allows for a faster initial relaxation, which is provoked as soon as the molecules are slightly nonspherical and depends weakly on elongation. The shear rigidity modulus, h` , characterizes the infinite frequency response of a system to shear. Any consistent attempt to calculate shear viscosity from theory must properly evaluate this property which is obtained from equilibrium considerations solely. In a molecular dynamics simulation, it is obtained through the following relation:

h `5

V 2 ^ @ J ag p # &, k BT

~27!

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

4764

MacDowell et al.: Dynamics of Kihara linear fluids TABLE X. ~a! Total (T) and potential ( P) part of the reduced shear rigidity modulus, h `* 5 h ` s 3 / e , as a function of L * for the different corresponding states of this work. Subcritical temperatures. The errors are found in the last decimal digit. ~b! As in ~a! but for the supercritical temperatures. ~a! L *

h `*

A1

A2

A3

B2

B3

T3

0.0

T P T P T P T P

8.87 8.22 8.82 8.40 7.66 7.37 7.26 7.03

22.6 21.7 18.4 17.9 16.7 16.3 15.9 15.6

27.5 26.6 22.0 21.4 20.9 20.5 21.4 21.1

19.2 18.5 17.3 16.9 15.3 15.0 13.9 13.7

25.8 25.1 21.3 20.9 18.6 18.3 18.4 18.2

24.0 23.4 19.6 19.2 17.4 17.2 17.3 17.1

0.3 0.6 0.8

~b! L *

h `*

C0

C1

C2

C3

0.0

T P T P T P T P

1.29 0.92 1.03 0.80 0.78 0.64 0.76 0.64

11.7 10.6 10.6 9.9 9.41 8.93 8.70 8.32

28.8 27.3 22.8 21.9 21.1 20.5 20.1 19.6

35.6 33.9 27.0 26.1 27.1 26.4 26.8 26.6

0.3 0.6 0.8

where a and g stand for any of the three Cartesian coordinates. From Eqs. ~2! and ~3!, it can be easily shown that h` can be split into a pure kinetic and a pure potential term. The kinetic term may be trivially calculated and shown to be equal to r k B T. In Tables X~a! and X~b!, we show the total h` , together with the potential contribution solely, obtained as h`2r k B T. It is apparent that the potential part of h` increases with increasing temperature and density, while it consistently decreases with elongation when compared with our correspondence criteria. The same dependence with elongation, density, and temperature is found for the kinetic term. For a long time, the value of the shear viscosity at the pseudo-triple point of the Lennard-Jones was the subject of much controversy. Very thorough work is now available in which it is suggested that this controversy was mainly due to a strong system size effect.44,54 According to these papers, the shear viscosity at the triple point for systems smaller than about 500 particles will be underestimated. However, the value obtained in this work, along with that of a recent paper by Stassen and Steele57 is higher than those of Refs. 44 and 54 ~in the former case, only after removing the extended mode coupling correction to the final result of 3.345! in which it is claimed that the thermodynamic limit was reached. In the case of Ref. 54, the discrepancy is certainly due to different estimations of the correlation time, since the shear rigidity modulus obtained in this work ~23.99 reduced units of pressure! is in good agreement with that reported i.e., 24.02 for G `,2 of that paper.

FIG. 7. ~a! Normalized stress correlation function for state C0 and the four elongations. ~b! As in ~a! but for the C3 state. ~c! As in ~a! but for the T3 state.

C. Thermal conductivity and the heat flux correlation function

Tables XI~a! and XI~b! show the results obtained for the thermal conductivity for the different states and molecular elongations.

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

MacDowell et al.: Dynamics of Kihara linear fluids

4765

TABLE XI. ~a! Reduced thermal conductivity, l*5ls2(m/ e k 2B ) 1/2, as a function of L * for the different corresponding states. Subcritical temperatures. Numbers in parentheses give a measure of the uncertainty in the last decimal digit. ~b! As in ~a! but for the supercritical temperatures. ~a! States Elongation L *50.0 L *50.3 L *50.6 L *50.8 ~b! States Elongation L *50.0 L *50.3 L *50.6 L *50.8

A1 l*

A2 l*

A3 l*

B2 l*

B3 l*

T3 l*

2.6~1! 2.5~2! 2.2~3! 2.1~2!

5.9~2! 5.0~1! 4.8~3! 4.0~3!

6.3~5! 5.4~1! 6.0~6! 5.4~5!

5.5~5! 4.9~5! 4.5~5! 3.5~2!

7.3~1! 6.2~4! 5.3~5! 5.5~5!

6.8~1! 6.15~7! 4.6~7! 4.6~2!

C0 l*

C1 l*

C2 l*

C3 l*

1.12~2! 1.04~6! 0.89~6! 0.78~2!

3.03~5! 3.2~1! 3.0~3! 2.4~2!

6.2~6! 6.5~7! 5.2~3! 4.7~7!

8.3~2! 6.4~9! 6.1~5! 6.4~9!

As observed for the shear viscosity, l* decreases with increasing elongation when compared with the correspondence criteria of this work. Once more, this decrease seems to be more pronounced when going from the Lennard-Jones to the L *50.3 model than in the rest of the cases. As expected, l* increases with density and temperature, although the temperature dependence as observed from the results of this work shows some irregularities which are probably due to the uncertainty in its determination. The normalized heat flux correlation functions of states C0 and T3 for all the elongations studied in this work are shown in Fig. 8. Figure 8~a! shows the same behavior as that observed for the stress correlation function at state C0. Although there seems to be no long time tail, the decay is so slow that it becomes necessary to integrate up to at least 3 reduced time units. Note also that if a modified Enskog theory47 was to predict l* for this state, one would expect this correlation function to decay as a single exponential,2 while Fig. 8~a! shows that this is not true even for the Lennard-Jones fluid. This argument also holds for h*, since Fig. 7~a! does not show evidence of a single relaxation time either. Figure 8~b! shows a steep initial decay of the heat flux correlation function but no evidence of a long time tail, at least of a size comparable with that of the stress correlation function. The steep decay varies quite drastically as soon as the molecules are nonspherical but from then on it does not seem to vary significantly with elongation. As was observed for the stress correlation function, the initial decay is always very similar for the elongated molecules and much steeper than that of the Lennard-Jones. Again, it seems that the greater number of available degrees of freedom allow a faster randomization of the energy fluctuations, so that the correlation times decrease with increasing elongation. On the other hand, in this case the correlation functions have a much simpler structure. The L *50.3 model shows a small shoulder at intermediate times, but this feature disappears completely for the more elongated molecules in the liquid state. In fact, only in the low density state does this function seem to show a complex structure.

FIG. 8. ~a! Normalized heat flux autocorrelation function for state C0 and the four elongations. ~b! As in ~a! but for state T3.

The fact that there is a discontinuity when one considers the thermal conductivity as a function of elongation is clearly illustrated by considering the thermal rigidity modulus, l` . This property is obtained from simulation by the following relation: l `5

V ^ J 2& . 3k B T 2 q

~28!

The results obtained for l` are shown in Table XII in reduced units. It is apparent that l` decreases on going from L *50.3 to L *50.8, while it clearly increases on going from L *50.0 to L *50.3. The relaxation mechanisms instigated in the elongated molecules are capable of canceling this sudden

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

MacDowell et al.: Dynamics of Kihara linear fluids

4766

TABLE XII. ~a! Reduced thermal rigidity modulus, l `* 5 l ` m s 3 /k B e , as a function of L * for the different corresponding states. Supercritical temperatures. The errors are found in the last decimal digit. ~b! As in ~a! but for the supercritical temperatures. ~a! States Elongation

A1 l `*

A2 l* `

A3 l `*

B2 l `*

B3 l `*

T3 l* `

L *50.0 L *50.3 L *50.6 L *50.8

23.7 37.6 32.8 29.8

61.4 79.5 69.6 64.4

73.9 92.3 89.1 89.4

52.9 77.6 63.4 58.8

73.4 94.5 79.5 80.1

64.6 88.6 70.8 71.2

~b! States Elongation

C0 l `*

C1 l `*

C2 l `*

C3 l `*

L *50.0 L *50.3 L *50.6 L *50.8

3.31 4.19 3.28 3.13

32.6 45.1 40.5 36.2

80.4 102.9 89.5 82.4

104.0 116.0 124.1 120.3

increase, resulting in an overall decrease of l with elongation. As for the density and temperature dependence, the same sort of behavior as h` is observed. VI. CONCLUSIONS

priate but it was shown here that a study of results of other authors supports this conclusion. It is expected, however, that slip boundary conditions will do better for more elongated molecules, as has been suggested by the work of Tang and Evans.51 Individual orientational relaxation has been studied by means of a simple model capable of predicting both the short time inertial regime and the long time diffusion. It is observed that the model is quite flexible in fitting the tested correlation functions but care must be taken when interpreting the physical meaning of the resulting parameters. It is observed that at high densities, rotational diffusion decreases with temperature, while the contrary behavior is seen for the orientational relaxation. ACKNOWLEDGMENTS

The Spanish DGICYT ~Direccio´n General de Investiga´ cion Cientı´fica y Te´cnica! is gratefully acknowledged for financial support ~Project PB94-0285!. One of us ~L.G.M! wishes to thank a predoctoral grant from Universidad Complutense de Madrid. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids ~Clarendon, Oxford, 1987!. 2 B. J. Alder, D. M. Gass, and T. E. Wainwright, J. Chem. Phys. 53, 3813 ~1970!. 3 C. Vega and S. Lago, Chem. Phys. Lett. 185, 516 ~1991!. 4 B. Garzo´n, S. Lago, C. Vega, E. de Miguel, and L. F. Rull, J. Chem. Phys. 101, 4166 ~1994!. 5 C. Vega, S. Lago, E. de Miguel, and L. F. Rull, J. Phys. Chem. 96, 7431 ~1992!. 6 B. Garzo´n, S. Lago, C. Vega, and L. F. Rull, J. Chem. Phys. 102, 7204 ~1995!. 7 J. P. Hansen and I. R. McDonald, Theory of Simple Liquids ~Academic, London, 1986!. 8 C. G. Gray and K. E. Gubbins, Theory of Molecular Liquids ~Clarendon, Oxford, 1984!. 9 C. Vega, S. Lago, and P. Padilla, J. Phys. Chem. 96, 1900 ~1992!. 10 C. Vega, S. Lago, R. Pospisil, S. Labik, and A. Malijevsky, J. Phys. Chem. 96, 1895 ~1992!. 11 K. R. Harris, Mol. Phys. 77, 1153 ~1992!. 12 J. Amoros, Mol. Phys. 83, 771 ~1994!. 13 R. Castillo and J. Orozco, Mol. Phys. 93, 343 ~1993!; J. Orozco and R. Castillo, J. Chem. Phys. 99, 1300 ~1993!. 14 R. K. Sharma, K. Tankeshwar, and K. N. Pathak, J. Phys. Condensed Matter 7, 537 ~1995!. 15 A. Z. Panagiatopoulos, Mol. Phys. 61, 813 ~1987!. 16 D. J. Evans and W. B. Street, Mol. Phys. 36, 161 ~1978!. 17 C. Hoheisel, J. Chem. Phys. 89, 3195 ~1988!; C. Hoheisel and H. Luo, Nuovo Cimento D 12, 499 ~1990!; H. Luo and C. Hoheisel, Phys. Rev. A 44, 6421 ~1991!. 18 D. J. Evans, Mol. Phys. 42, 1355 ~1981!; M. P. Allen and D. Kivelson, ibid. 44, 945 ~1981!; M. P . Allen and G. Marechal, ibid. 57, 7 ~1986!. 19 K. Singer, J. V. L Singer, and A. J. Taylor, Mol. Phys. 37, 1239 ~1979!. 20 P. Bereolos, J. Talbot, M. P. Allen, and G. T. Evans, J. Chem. Phys. 99, 6087 ~1993!. 21 S. Tang, G. T. Evans, C. P. Mason, and M. P. Allen, J. Chem. Phys. 102, 3794 ~1995!. 22 T. Kihara, J. Phys. Soc. Jpn. 16, 289 ~1951!. 23 G. Mare´chal and J. P. Ryckaert, Chem. Phys. Lett. 101, 548 ~1983!; M. P. Allen, Mol. Phys. 52, 705 ~1984!. 24 D. J. Evans and S. Murad, Mol. Phys. 68, 1219 ~1989!. Note that Ref. 14 already contains an expression for the heat flux. However, due to a misprint, the sign is changed in the translational potential term. 25 D. J. Evans, Mol. Phys. 32, 1171 ~1976!. 26 B. L. Holian and D. J. Evans, J. Chem. Phys. 78, 5147 ~1983!. 1

The results of this work show that both individual and collective time correlation functions of linear molecules exhibit characteristics which are peculiar to them and not found in simple spherical fluids. However, it was observed that for the normalized velocity correlation function, the short time behavior of molecules ranging from L *50 to L *50.8 is very similar when compared in corresponding states. The short time behavior of the normalized collective correlation functions of the linear molecules is also quite similar, though significantly different from that of the Lennard-Jones. The overall effect is that the corresponding short time integrals ~those appropriate for h* and l*! decrease with elongation. In the case of h*, this initial decrease may be somewhat compensated by the long time correlations which are seen to become more important as elongation increases. The heat flux correlation function does not show evidence of a long time tail for any elongation, so that this compensating effect is not expected. Overall, both h* and l* are found to decrease with elongation when compared with the correspondence criteria of this work. As to the translational motion, it was observed that the correlation times increase with elongation. This is due to the decrease of density and to the relatively long persistence of the motion parallel to the molecular axis. On the other hand, the normalizing factor of the Green–Kubo integral appropriate to D * ~i.e., T *! decreases accordingly so that D * is found to be almost independent of elongation in most of the liquid region. This was shown by means of a corresponding states plot of D * as a function of T/T c , where a linear dependence was observed. The Stokes–Einstein relation was tested and it was found that stick boundary conditions gave relatively good predictions of D * for the liquid densities. Traditionally it had been thought that slip boundary conditions were more appro-

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html

MacDowell et al.: Dynamics of Kihara linear fluids E. Helfand and S. A. Rice, J. Chem. Phys. 32, 1642 ~1960!. C. Vega and S. Lago, Comput. Chem. 18, 55 ~1994!. 29 The algorithm named ORT by D. Fincham, Mol. Simul. 11, 79 ~1993! and also described in Ref. 1. 30 C. Vega and S. Lago, J. Chem. Phys. 93, 8171 ~1990!. 31 D. J. Evans, W. G. Hoover, B. H. Failor, B. Moran, and A. J. C. Ladd, Phys. Rev. A 28, 1016 ~1983!. The implementation of this thermostat to the leap frog algorithm was done following Ref. 1 and D. Brown and J. H. R. Clarke, Mol. Phys. 51, 1243 ~1984! for the translational motion and D. Fincham, N. Quirke, and D. J. Tildesley, J. Chem. Phys. 84, 4535 ~1986! for the rotational motion. 32 D. J. Tildesley and P. A. Madden, Mol. Phys. 48, 129 ~1983!. 33 M. P. Allen and A. J. Masters, Mol. Phys. 79, 435 ~1993!; M. P. Allen, D. Brown, and A. J. Masters, Phys. Rev. E 49, 2488 ~1994!. 34 R. Zwanzig and N. K. Ailawadi, Phys. Rev. 182, 280 ~1969!. 35 A. Lotfi, O. Vrabec, and J. Fischer, Mol. Phys. 76, 1319 ~1992!. 36 B. J. Alder and T. E. Wainwright, Phys. Rev. Lett. 18, 988 ~1967!. 37 R. Zwanzig and M. Bishop, J. Chem. Phys. 60, 295 ~1974!. 38 J. Barojas, D. Levesque, and B. Quentrec, Phys. Rev. A 7, 1092 ~1973!. 39 M. P. Allen, Phys. Rev. Lett. 65, 2881 ~1990!. 40 E. de Miguel, L. F. Rull, and K. E. Gubbins, Phys. Rev. A 45, 3813 ~1992!. 41 D. Frenkel and J. F. Maguire, Mol. Phys. 49, 503 ~1983!. 42 R. M. Lynden-Bell and W. A. Steele, J. Phys. Chem. 88, 6514 ~1984!. 27 28

4767

43

A. I. Burshtein and S. I. Temkin, Spectroscopy of Molecular Rotation in Gases and Liquids ~Cambridge University Press, Cambridge, 1994!. 44 J. J. Erpenbeck, Phys. Rev. A 38, 6255 ~1988!. 45 T. Boublik, J. Chem. Phys. 87, 1751 ~1987!. 46 D. Levesque and L. Verlet, Phys. Rev. A 2, 2514 ~1970!. 47 H. J. M. Hanley, R. D. McCarty, and E. G. D. Cohen, Physica 60, 322 ~1972!. 48 H. Eyring, D. Henderson, B. J. Stover, and E. M . Eyring, Statistical Mechanics and Dynamics, 2nd ed. ~Wiley, New York, 1972!. 49 R. Zwanzig, Phys. Rev. A 2, 2005 ~1970!. 50 J. A. Montgomery, Jr. and B. J. Berne, J. Chem. Phys. 67, 4580 ~1977!. 51 S. Tang and G. T. Evans, Mol. Phys. 80, 1443 ~1993!. 52 P. Borgelt, C. Hoheisel, and G. Stell, Phys. Rev. A 42, 789 ~1990!. 53 D. M. Heyes, J. Chem. Soc. Faraday Trans. 74, 1741 ~1983!. 54 M. Schoen and C. Hoheisel, Mol. Phys. 56, 653 ~1985!. 55 W. J. Moore, Physical Chemistry ~Prentice-Hall, Englewood Cliffs, 1963!. 56 See, for example, G. Mare´chal, J. P. Ryckaert, and A. Bellemans, Mol. Phys. 61, 33 ~1987! or C. J. Mundy, J. I. Siepmann, and M. L. Klein, J. Chem. Phys. 102, 3376 ~1995! for an illustration of this effect in the particular case of n-alkanes. 57 H. Stassen and W. A. Steele, J. Chem. Phys. 102, 932 ~1995!. 58 R. Vogelsang, C. Hoheisel, and G. Ciccotti, J. Chem. Phys. 86, 6371 ~1987!. 59 M. F. Pas and B. J. Zwolinsky, Mol. Phys. 73, 471 ~1991!.

J. Chem. Phys., Vol. 106, No. 11, 15 March 1997

Downloaded¬26¬Feb¬2001¬to¬145.18.129.48.¬Redistribution¬subject¬to¬AIP¬copyright,¬see¬http://ojps.aip.org/jcpo/jcpcpyrts.html