J. Phys. Chem. 1987, 91, 21-30 0 FIR,)
2
1
-
Lnm
Figure 1. Emission spectrum (1) and excitation spectrum (2) of the luminescence at 4.2 K of Bi3Re08. @A gives the spectral radiant power per constant wavelength interval and qr the relative quantum output, both in arbitrary units. Curve 1 was recorded for 325-nm excitation. Curve 2 was recorded for 560-nm emission. The diffuse reflectance spectrum (3) was measured at 300 K and plotted as the Kubelka-Munk function F(R.).
see ref 3). Luminescence of this type is characterized by the unexpected combination of a very large Stokes shift (>2 eV) and a high quenching temperature (approximately room temperature). Some illustrative data are given in Table I. This, however, is not observed for the luminescence of Bi3Re08. The Stokes shift is smaller than 2 eV and the thermal quenching temperature is much lower than room temperature. Therefore, this possibility is rejected. We are left then with a third possibility, viz. Bi3+-Re7+ charge-transfer transitions. This is in analogy with the lu-
21
minescence of bismuth tungstates and molybdates. This type of luminescence is characterized by broad band emission with low quenching temperatures (of the order of 100 K) and Stokes shifts of about 1.4 eV.4 Examples are given in Table I. These values are of the same order of magnitude as observed for the emission of Bi3Re08. Actually, the charge-transfer transition within the perrhenate group itself is situated at about 280 nm (in KRe04) with tetrahedral Reo4- (ref 2) as well as in EU3ReOs with octahedral Re065-(ref 7)). This value is not much different from that in the molybdate or tungstate group.2 Just as the optical absorption edge of bismuth molybdates and tungstates has shifted to longer wavelengths compared to that of molybdates and tungstates without s2-configurationcations,2 so the absorption edge of Bi3Re08is at longer wavelengths than in KRe04 or Eu3Re08. Following earlier models2 this shift is ascribed to a low-energy absorption transition characterized by charge transfer from the 6s orbital of Bi3+to the empty 5d orbital of Re7+. The emission consists of the reversed transition. If this assignment is correct, the luminescence of Bi3Re08 presents the first example of an emission transition in which the Re7+ ion plays a role. Up till now perrhenates appeared to be nonluminescent.2 In view of all the data on luminescent and nonluminescent tetroxo ions with a closed-shell configuration, the absence of luminescence in the perrhenate group has to be ascribed to a very large reorganization of the excited state, making crossover to the ground state a probable p r o c e ~ s The . ~ ~present ~ ~ ~ results suggest that this reorganization is smaller in Bi3Re08. This is confirmed by comparable observations for the luminescence of the bismuth tungstates and molybdates and has been ascribed to the antibonding nature of the 6s2 electrons of the Bi3+ ion.2
FEATURE ARTICLE Description of the Molecular Trajectories in Simple Liquids Cheruvu S. Murthyt Department of Chemistry, The University of Texas at Austin, Austin, Texas 78712
and Konrad Singer* Department of Chemistry, Royal Holloway and Bedford New College, Egham, Surrey TW2O OEX, England (Received: March 1 1 , 1986; In Final Form: August 14, 1986)
In order to gain insight into the nature of molecular motions and the mechanism of diffusion, attempts have been made to use molecular dynamics simulation data primarily in two ways: (1) the analysis of the velocity autocorrelation function, and (2) the division of trajectories into “steps” according to different criteria. This paper gives a selective summary of the literature and presents data on the analysis of translational and rotational trajectories in terms of the quasi-oscillations and librations of the Cartesian components of the translational and angular velocity vectors, and in terms of right-angle turns of these vectors. The curvature of the local potential field and energy fluctuations are discussed briefly.
I. Introduction The simulation of classical many-particle systems by Monte Carlo and molecular dynamics (MD) simulation has 1 4 to great advances in the understanding of liquids and solids during the past 25 years. The derivation of all physically significant mean values t Present address: Chemical Dynamics Corporation, 9560 Pennsylvania Avenue, Upper Marlboro, MD 20772.
0022-3654/87/2091-0021$01.50/0
and correlations in space and time from a given interaction potential in many-particle systems has made it possible to test the adequacy of pair potentials and theories, and to assist in the interpretation of experiments, s-x”nmes supplementing them by otherwise inaccessible data, such as structural and time correlations in polyatomic liquids, interfacial properties, etc. Although the relaxation of fluctuations described by time autocorrelation functions (AF) leads to important insights into the 0 1987 American Chemical Society
22
The Journal of Physical Chemistry, Vol. 91, No. 1, 1987
TABLE I: Potential Parameters, State Points, Computational Details“ system elk UIA P* MI) Ar(I1) Ar(II1)
3.405 3.405 3.405 2.825 3.332 3.332
119.8 119.8 119.8 52.8 178.3 178.3
F2 C12U) C12(W
0.835 0.731 0.626 0.608 0.539 0.497
Murthy and Singer
I*
T*
(0.0) (0.0) (0.0) 0.505 0.63 0.63
0.810 0.977 1.135 1.06 1 .oo 1.35
Atls 2 x 10-14 2 x 1044 2 x 1044 0.5 x 10-14 1.0 x 10-14 1.0 x 10-14
time span/ps 50
50 50 8 16 16
‘For argon: I-center Lennard-Jones (LJ) potential y, = 4 ~ [ ( u / r i > )-’ ~( ~ / q ~ )For ~ ]F2, . CI,: 2-center LJ potential, Le., a site-site LJ potential between each atom in molecule i with each atom in molecule j. p* = NAv u3/V,; V, = molar volume; I* = l / u ; I = bond length; T* = kT/e; At = time step in the MD simultations; time spanlength of simulations which were carried out with 256 particles.
dynamics of liquids, M D simulation yields a vast amount of information which is discarded by the essential averaging: the complete dynamical history of the system consisting typically of 100-1000 particles over a time span of 5-100 ps is reduced to lo3 numbers. The problem of obtaining from this wealth of information an answer to the question “how do molecules in a liquid move?” is far from simple. In the gas phase and in hard-sphere (HS)fluids, collisions provide a natural criterion for the subdivision of trajectories which is not available in “real” liquids; the complete absence of regular features in the translational and rotational motion of dense liquids is illustrated in Figure la,b. The motivation to search for an intuitively meaningful description of trajectories is threefold: (1) scientific curiosity, (2) to test the physical basis of models which have been proposed, and (3) to assist in the search for realistic new models. The search for a detailed and useful description of trajectories has proceeded mainly along two lines. One is the analysis of the velocity A F (VAF) according to various criteria; the other, which we explore, is the subdivision of trajectories into physically meaningful sections and their statistical analysis. The plan of the paper is as follows. Section I1 is a selective summary of papers in which the VAF is analyzed or resolved in various ways. In section 111, information is obtained from the subdivision of the particle trajectories according to different criteria: for the hard-sphere fluid from the distribution of free paths between co1lisions;l for monatomic and diatomic Lennard-Jones (LJ) liquids from the space and time intervals produced (1) by the zeros of the projection of the velocity vector v onto a fixed axis, and (2) by changes of the direction of v by 90’. Sections IV and V deal with the rotational motion of diatomics simulated by 2-center LJ potentials; they are analogous to sections I1 and 111 but with “rotational” substituted for “translational.” Section VI discusses the curvature of the potential field acting on a particle and, briefly, fluctuations of the kinetic and potential energies. A summary of conclusions is given in section VII. Our work deals with monatomic and diatomic liquids simulated with 1- and 2-center LJ potentials. Some of the material on the former has been reported previously.2 The systems are specified in Table I .
-
11. The Velocity Autocorrelation Function A . General Features. The diffusion constant of an isotropic
fluid is determined by D = ( L ; ~ ) ~ ~ dt C ~ ~ ( C )
(2.1)
where C , is the normalized VAF: (2.2) where v is the velocity vector of the center of mass (COM) and ( ) denotes ensemble or time averages (the latter over the time origins t = 0). The VAF cannot be determined experimentally; its importance lies not only in eq 2.1 but also in the light which it sheds on details of the molecular motion. The variation of the VAF with density of the HS fluid is shown in Figure 2a. The same qualitative features apply to fluids with C,(t) = (d)-’(v(O).v(t))
= ( u , 2 ) - ~ ( u x ( 0 )u,(t))
(1) Einwohner, T.; Alder, B. J. J . Chem. Phys. 1969, 49, 1458. (2) Murthy, C. S.; Singer, K. Proe. R . Soe. A 1983, 389, 299.
b
b
Figure 1. (a) Projections of typical trajectories in LJ argon onto the xy plane. Duration: 4 ps; distance between dots: 1 ps. T* = 0.81; p* = 0.835 from ref 2. (b) Stereoscopic pictures of trajectories generated by the rotation of the bond vector of a linear molecular (F2) on the surface of a sphere. Duration: 5 ps; T* = 1.06, p* = 0.608, I* = 0.505.
simple continuous interaction potentials, eg., the repulsive “soft sphere” r-I2, and 1- and 2-center LJ potentials. At high density (near the triple point in real liquids), the curve passes through a shallow negative minimum and then rises slowly to the asymptotic zero when the correlation with the initial velocity has disappeared (Figure 2b). With decreasing density (and increasing temperature) the negative lobe disappears and the VAF decays monotonically and more rapidly to zero. The passage through zero to negative values at t = T~ shows that the majority of trajectories have turned through more t h a n
The Journal of Physical Chemistry, Vol. 91, No. I, 1987 23
Feature Article 0.08
1.o
0.04 0.5 Q0
e
LL
a >
0.1
\
-0.04 0
1.5
I
0.6
1.0
=hulkBT
-0.1
/ \
0:2
i
0
I
\
I
1
0
,'
'-' 10
20
30
0
1.o
0.5
1.5
TIME I ps
S
7.5
c
.03 5.0 N
.02
oa \
5?
-
3
U
>
N
.01
'-
a
Q
V
2.5
0
5
10
15
TIME I ps
Figure 2. (a) VAF of the HS fluid at different densities. E denotes the exponential decay predicted by the Enskog theory; positive and negative deviations from this occur at low and high density; from ref 3. (b) VAF from the first MD simulation of liquid Ar. The continuous curve is the average over 64 time origins. Open circles and crosses represent maximum departures from the mean. Full circles show the exponential decay calculated with the Langevin equation. The inset shows F T s of the MD- and the Langevin-VAF; from ref 5 . (c) A portion of the VAFs ,C, CWl1,and CWL(full line, dashes, and dot-dashes) for 2-center LJ Nz;from ref 17. (d) Mean squared displacements produced by vi1 and v l at different state points in simulations of CS2. At 192 K 10 ps elapse before the slopes are equal! From ref 18. 90° after
T~ because of back-scattering from near-neighbors. Because of the spread of reversal times, the negative lobe is shallow (usually ((1); it disappears at lower density because of the greater spread of reversal times and the enhancement of forward as opposed to back-scattering produced by a larger and looser cage of near neighbors. Figure 2a also shows the VAF predicted by the Enskog theory, which is based on the assumption that successive collisions are uncorrelated. While back-scattering accounts for the negative deviations a t high density, the explanation of the positive deviations or faster diffusion at low densities is more subtle. It has been accounted for by inspection of the MD trajectories and by hydrodynamic calculation-by a vortex back-flow from the region of higher density in front of a displaced particle to the
lower density regions in its wake.4 In the first MD simulation with a "realistic" pair potential for argon,5 Rahman showed that the VAF cannot be accounted for by the Langevin equation which is valid for Brownian motion. He also obtained the frequency distribution of the quasi-oscillations from the Fourier transform (FT) of the VAF (Figure 2b). The value of w,, agrees quite well with the estimate of the cycle time calculated on the assumption that the interval 7,,when the VAF (3) Alder, B. J.; Gass, D. M.; Wainwright, T. E. J . Chem. Phys. 1970, 53, 3813. (4) Alder, B. J.; Wainwright, T. E. Phys. Reu. A 1970, 1 , 18. (5) Rahman, A. Phys. Rev. A 1964, 160, 217.
The Journal of Physical Chemistry, Vol. 91, No. I, 1987
Murthy and Singer
passes through zero, corresponds to 1 / 4 of the complete quasioscillation ( T = ~1.3 ps). ~ ~ ~ ~ An entirely different way of obtaining the cycle time will be discussed in section VIA. For LJ argon the decay time of the VAF near the triple point is -2 ps. The effect of the cohesive part of the potential is small but not negligible (as it is for the radial distribution function;6 for a wide range of state points the VAF of the repulsive r-I2 potential is similar but lies slightly above that of the full LJ potential.' This is expected since cohesive' forces decelerate or prevent the motion of adjacent particles in opposite directions and thus reduce the diffusion constant. The VAF for the COM of diatomics simulated by 2-center LJ potentials has the same qualitative features as that exhibited by monatomic fluids, except for a second shallower minimum near the triple point.* B. Further Analysis of the VAF. I. Rahman's Description of Diffusion. To gain insight into the mechanism of diffusion Rahman analyzed the VAF as f01lows:~let the displacement vector of a particle in the time interval T (-0.5 ps) be Ar(7) = r(T) r(0) and let the {-axis be parallel, and the {- and ?-axes be normal to this direction. v(O).v(t) = S,(t) + R,(t)
3. Transfer of Momentum to Neighbors. When the VAF curve passes through zero after t = T ~ the , initial momentum has been transferred to the surrounding particles. If diffusion is prevented by freezing the nearest neighbors in a typical liquid configuration, the minimum of the VAF deepens so as to reduce the time integral (eq 2.1) to zero. In addition, the VAF begins to deviate from that of particles in ordinary (unconstrained) liquid well before t = T ~ This shows that the relaxation of the cage around a particle occurs 011 the same time scale as the quasi-oscillation of the central particle." The transfer of momentum from a particle to its neighbors has been studied by means of the velocity cross-correlation function (VCF)
24
where S,(t) = ur(0)uqt); R,(t) = uf(O)uf(t)
+ u"0)u"t)
+
S,(t) is further resolved into S,+(t) S;(t) according to whether d ( 0 ) is parallel (+) or antiparallel (-) to Ar(s). The VAF R,(t) on the components S,+(t) on the one hand, and SJt) other, are very different. S,+(t) passes through zero after T~ = 0.7 ps to a very shallow minimum and then slowly rises to zero. The curve for S;(t) + R,(t) passes through zero after 0.2 ps and decays after reaching a deep minimum (-0.33). The time integral over this curve is virtually zero whereas that over S,+(t) is strongly positive. Rahman refers to the two components as "slipping" and "rattling"; diffusion is mainly due to the former. This characterization is confirmed by the Fourier transforms whose maxima are at 2 and 7 THz, respectively. If the interval T is chosen equal to the average time in which the elongation of nearest-neighbor cage changes from 0 to its maximum value, Ar(7) is found to be roughly parallel to the elongation. Diffusion may therefore be described as slipping in the direction favored by the shape of the cage, accompanied by rapid rattling which contributes little to it. 2. The Effect of the Local Density. The resolution of the VAF of a soft-sphere fliud [V(r)= ( u / r ) ] I saccording to the number n of neighbors within the radius 1.480 has yielded interesting results.1° The mean value ti is 11.8. The VAF for n = 9 or 10 is strongly oscillatory, exhibiting more than one minimum and maximum. For n = 13 or 14, there is only one minimum after which the curve rises monotonically through a broad shoulder to zero. The velocity cross-correlation function between the central particle and its neighbors is positive (see below). If the sphere of neighbors is divided into a front and rear half in relation to the direction of ul, one finds that for n < ri the correlation of v1 with neighbors in the rear half is stronger than with those on the front half, while the opposite is true for n > ri. The suggested explanation is that, depending on whether n is smaller or larger than ri, the particles in the nearest shell will develop a radial velocity directed inward or outward. If outward, they recede from the central particle and the chance of the reversal of its momentum is reduced; if inward, the converse is true. Moreover, in a region of lower density adjacent particles may oscillate more easily.
+
(6) Andersen, H. A.; Chandler, D.; Weeks, J. D. Adu. Chem. Phys. 1976, 34, 105.
(7) Kushick, J.; Berne, B. J. J . Chem. Phys. 1973, 59, 3732. (8) Singer, K.; Singer, J. V. L.; Taylor, A. J. Mol. Phys. 1979,37, 1239. (9) Rahman, A. J . Chem. Phys. 1966, 45, 2585. (IO) Endo, Y . ; Endo, H. J . Chem. Phys. 1984, 80. 2087.
where the second particle lies within the first shell ( 0 . 9 ~5 r12 < 1.5g).I2 The VCF reaches a maximum at T~ when the central particle has completely lost its initial momentum. During this period -70% of the momentum is transferred to nearest, and 30% to next nearest neighbors (1.50 Ir2 < 2 . 4 ~ ) .Another important result is the fact that -90% of the transfer is along the r12vectors; it is longitudinal rather than transverse. Collective current fluctuations in LJ liquids have been analyzed in terms of the longitudinal and transverse components of the FT in space or space and time.I3*l4 The single-particle motion is obviously coupled with collective modes; thus, a long-wave fluctuation will move adjacent particles in the same direction whereas for a wavelength equal to twice the nearest-neighbor distance, neighbors move in opposite directions. Slow density fluctuations may account for the long negative or positive tail of the VAF. it is, however, difficult to link the description of the single-particle motion to the collective modes, and a discussion of the latter is beyond the scope of this article. The direct comparison, however, of the mean squared displacement (MSD) of a particle with the collective MSD or the COM of its neighbors lying within the first-shell of radius 1 . 5 (at ~ t = 0), throws light on the mechanism of diffusion. Balucani et a1.lSfind that the relative slopes of the MSDs of the single particle rl, of its position relative to the COM and of the COM of the cage L,, are 5:3.2:1. of the cage rl - b, This shows that there is a positive correlation between the displacements of a particle and of its cage and that the MSD of the particle and the collective MSD of the COM of its cage are of comparable magnitude. It is also found that the rotation of the cage contributes to the MSD of the particle. The fact that particles move in groups or clusters is convincingly borne out by recording the labels of nearest neighbors over a period of time. For L J argon near the triple point, more than 60% of neighbors remain so during the decay time of the VAF (-2 ps) and clusters of 2-4 persist over 4 times that period. At lower density and higher temperature the turnover is much more rapid; but even at p* = 0.626 and T* = 1.135, small clusters persist over 2 ps.16 C. The VAF ofLinear Molecules. For nearly spherical linear molecules such as N2 the VAF is similar to that of spherical particles. When the anisotropy (measured by the ratio of bond length to atomic diameter 1* = l / ~is) greater than 0.5 (F,, Cl,, CO,, CS,), the minimum is shallower and disappears not far below the density of the triple point.* As expected, the mobility parallel to the bond axis at t = 0 is greater than at right angles to it. This was first shown by Barojas et aI.l7 by resolving the velocity at t = 0 into components which are parallel and normal to the unit bond vector u = 1 / 1 at t = 0: v"(t) = [v(r)-u(O)]u(O) and vl(t) = v(t) - d ( t )
(2.3)
(11) Dean, D. P.; Kushcik, J. J . Chem. Phys. 1982, 76, 619. (12) Balucani, U.; Vallauri, R.; Murthy, C. S. Phys. Left. 1984, 84, 133; J . Chem. Phys. 1982, 77, 3233. (13) Ailawadi, N. K.; Rahman, A.; Zwanzig, R. P.hys. Reo. A 1971, 4 , 4. (14) Levesque, D.; Verlet, L.; Kiirkijarvi, J. Phys. Reu. A 1973, 7, 1690. (15) Balucani, U.;Gori, M.; Vallauri, R. Chem. Phys. Left. 1985, I 1 9, 152. (16) Singer, K., Murthy, C. S., unpublished results. (17) Barojas, J.; Levesque, D.; Quentrec, B. Phys. Reu. A 1973, 7, 1092.
.
The Journal of Physical Chemistry, Vol. 91, No. 1, 1987 25
Feature Article The corresponding VAFs, Le., Cwll(t) and C,TL(t)are shown in Figure 2c. For Cwll, T~ is larger and the negative lobe is shallower than for CwL. As inspection of the curves suggests, DIl is larger than D , when they are calculated from the integral over the VAFs. Yet the alternative calculation of DIland D , from the slopes of the MSDs (Figure 2d) shows them (after rather long initial interval) to be equal-as they should be since the diffusion constant in an isotropic liquid does not depend on the molecular orientation at any particular timesi8 The long life of the correlation between orientation and translation for an extremely anisotropic interaction potential, manifested in the integrals for Dlland D,, has been attributed to local structural anisotropy, Le., an alignment of molecules so that those that keep their original orientation find it easier to slip through the approximately parallel neighbors.lg The randomization with time of the distribution of the angle between v(O)v(t) is discussed in a recent paper.”
1.2
1 .o
+
(18) Tildesley, D. J.; Madden, P. A. Mol. Phys. 1982, 48, 129. (19) Adams, D. J. Chem. Phys. Lett. 1982,85, 131. (20) Lynden-Bell, R. M.; Hutchinson, D. J. C.; Doyle, M. J. Mol. Phys. 1986, 58, 307.
(21) Alder, B. J.; Wainwright, T. E. J . Chem. Phys. 1959, 31, 459. (22) Glasstone, S.; Laidler, K.; Eyring, H. Theory of Rate Processes; McGraw-Hill: New York, 1941. (23) Barker, J . A.; Henderson, D. Rev. Mod. Phys. 1976, 48, 587; cell and free volume theories are reviewed in section VIII, p 661. (24) Franchetti, S. Nuovo Cimenro SOC.Ital. Fis., B 1980, 58B, 135.
h I
0.8
0.6
111. Analysis of Translational Trajectories
A . The Distribution of Free Paths in Hard-Sphere Systems. After setting out the principle of the M D method,*’ Alder and collaborators used M D simulation in a series of fundamental studies of hard-sphere, hard-disk, and square-well systems. Of particular significance in the present context is a paper which deals with the distribution of the lengths of free paths.’ We quote from the Abstract: “Mean free paths for solids and dense fluids are found to be within 2% of those predicted by dividing the kinetic theory free paths by the relative probability of two molecules in contact. These mean free paths are much smaller than a mean intermolecular spacing or even the diameter of a free volume. For hard-sphere and hard-disk fluids the free path distributions are found to be monotone decreasing, nearly exponential, and when scaled by the mean free path, nearly density independent for all fluid densities.” The free-path (FP) distribution [shown in Figure 3a] is incompatible with theories which postulate activated jumps of the order of magnitude of the interparticle spacingz2 because the distribution is unimodal and because free paths of sufficient length do cot occur. Cell theories23 are ruled out because for such models the FP distribution would not decay monotonically. It is shown that “large jumps” (of the order of magnitude of the nearest-neighbor spacing) do not occur even in a solid from which 1% of particles are removed. It has been suggested that these arguments may not apply to liquids in which molecules interact with realistic pair potentials.24 The results reported below show that the distribution of amplitudes of quasi-oscillations in the LJ liquid are qualitatively similar to FP distribution in the HS fluid. In the following sections we analyze the trajectories of 1- and 2-center LJ liquids simulating argon, FZ,and C12. Most of the 1-LJ data are taken from ref 2. Those for the diatomic liquids are new. Table I lists the potential parameters and state points and some computational data. Further details are given in ref 2 and 25. B. Quasi-Oscillations Determined by Changes of Sign of the Velocity. Let t , , t 2 , ..., t , be times at which the x component of the velocity u, of particle i equals zero. The time interval between the mth and the (m n)th incidence of u, = 0 (3.1) is T,,(,) = tm+,,- t,. The corresponding coordinate interval is (3.2) t n ( m ) = lx(tm+n) - x ( t m ) l
8
Q 1
0.4
0.2 rlA
\fLX-
I
0
2
1
3
4
rlA
c
b
0.4
0.2
0
4
0.6
Figure 3. (a) Reduced free-path distribution in HS systems. Triangles refer to the fluid (V/Vo= 1.5), the continuous curve to the solid (V/Vo = 1.48) where 5% of the particles have been removed, and the crosses to the solid V/Vo = 1.5; from ref 1. (b) Unnormalized probability densitiesflt,) (dots) a n d f l b ) (crosses) and P([,’>[) (dots inside circles) and P ( A r 3 A r ) (crosses inside circles) for F2.
Unless stated otherwise, the unit of length is the LJ parameter 0.
Since the velocity changes sign on passing through zero, successive “steps” are in opposite directions. The displacment after t w o steps starting at t , is t 2 ( m ) = Itl(m) - tl(m+l)l (3.3) Liquids in equilibrium are stationary systems; the mean values of r i m )and En(,) are therefore independent of the index m which specifies the time origin:
(F,’”)
=
(tn),
( r i m ) )=
(7,)
(3.4)
The same criteria are adopted t o characterize the rotation of the linear molecules; u, is replaced by the component of the angular
26
Murthy and Singer
The Journal of Physical Chemistry, Vol. 91, No. 1, 1987
TABLE I 1 Characteristics of the Cartesian Components
(a) Quasi-Oscillations‘ system
criterion
(Tl)/PS
=0 v, = 0 v, = 0 u, = 0
0.45 0.34 0.32 0.34
0,
F2 CMI) Cl,(II)
(€2)/.
(€I)/.
0.14 0.1 1 0.1 1 0.14
(Sl)/.
0.14 0.10 0.11 0.15
($1
0.29 0.21 0.22 0.27
)/deg 13 11 14
/deg
($2)
(01 ) P e g
19 16 21
37 32 39
(b) Quasi-Librationsb
system F2
criterion w,
=0
w, = 0 wx = 0
CMI) Cl,(II)
(TI)/PS
(vl)/deg
(72)/deg
( 4)/deg
0.34 0.32 0.34
19 15 19
18 15 20
38 31 39
‘The relative fluctuations s ((A2) - ( A ) 2 ) 1 / 2 / ( A )are s ( T ~ )N 0.65-0.80; N
s(&)
N
S(t2)
E
(€I1)/.
(c;,r)/.
0.07 0.08 0.10 1.0-1.2;
SI)
(Sl1)/0
0.11 0.12 0.15 0.6-0.9.
bs(~l)
0.21 0.22 0.27 0.7;
S(al) E
S(a2)
1.0; ~ ( 0 , )N 0.8-0.9.
velocity w, or by zi, where u is the unit vector along the molecular axis. The angular intervals (arcs) are denoted by 7. The mean values ( T ~ ) (E1), , and ( E 2 ) and the mean lengths of step trajectories (3.5) are given in Table IIa. In addition we list the analogous quantities for the diatomics, the angular intervals ql, and path lengths O1 defined in section VA. The following points emerge: (1) The mean coordinate intervals for one step, Le., the amplitudes of the quasioscillations (E1)are small compared with the nearest-neighbor distance: 0.14-0.24 for LJ argon compared with 1.1, and 0.11-0.14 compared with 1.4 for the L J diatomics. (2) (6,)N ( E 2 ) ; the mean displacements after one and two steps are roughly equal in all cases. For an oscillation about a fixed point, (Ez) would be zero, and for an itinerant oscillation about a slowly moving center, (E2) would be small compared with (El). It follows that the center of oscillation or of the cage moves with a velocity which is comparable to that of the oscillating particle. The near-equality ( E2J N ( Ezlrl ) is found to persist up to at least n = 9. (3) The fluctuations of all quantities are large; mean values and relative fluctuations are of similar magnitude. (4) The amplitude ( E , ) , which is a rough measure of the diameter of the cage, increases roughly as the square of the molar volume. The corresponding ( r l ) increases because the thermal velocity Iu,I r ( f l ) / ( r l increases ) more slowly (T1/*). (5) The accelerations change sign more frequently than the velocities; the ratio is 1.5-2 for L J argon and r 2 for the diatomics. For a harmonic oscillation the ratio is 1. (6) The ratio of path length (SI)to amplitude (E1) is r 2 . If the mean displacements along the x, y, and z axes during ( r l ) were equal, the rectilinear displacement in 3D would be 31/2(E1). The additional path length is attributed to curvature. Distribution of Coordinate and Time Intervals. The probability density f(&) and the distribution function P(FI’>E1) =
fflEl’)
dEl’/xmflEl’) dE1’
alternate with large jumps (“large” = similar to the nearestneighbor spacing). The probabiiity densityf(rl) and distribution function P(T1’&1) for F2 is similar to that of Ar(I);Z f(rl) rises from zero to a maximum at 2(r1)/3; the incidence of 7,2 3( r l ) is small. The qualitative difference between Art)andflfl) stems from the fact that very small 6, are associated with very small momenta, Le., comparatively long time intervals between turning points; very small r,-values are therefore very improbable. C. Right-Angle Turns. An alternative criterion for the subdivision of trajectories is V(t1).V(tl+711) = 0
(3.7)
which determines the intervals r l l in which a particle loses the momentum vector which it had at t , , that is, intervals during which the direction of the velocity turns through 90’. The corresponding intervals (now in 3D space) have the length Arl = lr(t1+rlL) - r(t,))
(3.8)
The distance covered after two successive right-angle turns is Ar2 = I r ( t l + ~ 2 i )- r(tl)l
(3.9)
The relevant statistics are given in Table IIIa. The main features are the following. The time intervals ( rL) are slightly smaller than the ( r l ) for quasi-oscillations. Since one is dealing with displacements in three rather than one dimension ( Ar2)/ ( A r , ) is larger than 1; it is ~ 1 . in 5 all cases. The mean displacments (Arl) are larger than the mean amplitudes ( E , ) because the initial velocity is not zero but sampled from the Maxwell distribution, and because the motion is 3-dimensional. The ratio of path length to displacement ( S l ) / ( A r l )is in all cases ~ 1 . 2 . The last column of Table IIIa lists the first zero times ro of the VAFs. At high density they are smaller than ( i1 ) but tend to if a negative lobe is absent. The different behavior is not surprising: ro is the time at which positive and negative values of v(t).v(t+r) cancel each other, whereas ( r I L is ) the average time at which eq 3.7 applies. The formal definitions are
-
(3.6)
for F2 are shown in Figure 3b. The decay is monotonic and approximately exponential. For F2, P ( [ , ’ 2 ( , ) = exp(-8.3E1) whereas for Ar(I), although similar, the decay at small is faster than can be fitted by an exponential curve (Figure 2 of ref 2). The following numbers show that displacements of order of magnitude of the nearest-neighbor distance (1.1 = 8( El) for Ar(I), and 1.4 = 12(5,) for F2) are very improbable indeed: €1
5d m’ 2 5,) p(tl’2
ArU) F2
= 2.5(t1) 0.05
0.1
€1
= 4(&) 0.01 0.026
€1
= 6(h) 0.001 0.004
The curvefl&) is rather similar to the distribution of free paths in the HS system (Figure 3a). Here as there the results provide strong evidence against any models in which small oscillations
The difference arises from the position of the Dirac &function outside the summation over i in (3.10) and inside in (3.11). The time and space intervals resulting from the substitution of vll and v* in eq 3.7 are respectively larger and smaller than those for v shown in the table; the differences are smaller than one might
The Journal of Physical Chemistry, Vol. 91, No. 1, 1987
Feature Article
21
TABLE III: Intervals Corresponding to Right-Angle Turns of the Velocity Vector
(a) Translational Motion' variable in (3.7)
(71 )/PS
W I ) / U
ArU)
system
V
F2
V V
0.40 0.30 0.30 0.34
0.21 0.15 0.17 0.22
CI2(I)
c12m
V
(Ar2) /. 0.30 0.23 0.25 0.33
(SIL)/U
ro(VAF)/ps
0.26 0.18 0.20 0.27
0.33 0.22 0.22 0.52
(b) Rotational Motionb
system
variable in (3.7)
(71 )/PS
W
0.28 0.28 0.30
F2 C12(1) CI,(II)
W
w
( A 4 )/deg 28 24 30
'Relative fluctuations defined as in Table I1 are 0.75-0.85 for s(At9,) u 0.7-0.9; ~ ( 0 , ) 0.8-1.0; s defined as in Table 11.
(71).
s(Arl) z s(Ar2)
(AW/deg
(&)/deg
70WAF)/PS
40 33 44
31 27 35
0.20 0.22 0.23
= 0.55-0.85, and 0.8-0.9 for (SI). bs(71)= 0.6-0.8; s(Agl)
The right-hand side of the table shows gradients estimated by the difference quotients
TABLE I V Convergence to the Diffusion Constant
(a) Translational motion n
Ar(1)
D,O
6.6
C12(I)
Di b,, D,
C12(II)
en
F2
Bn D,
2 -
4
a0
1.9 2.5 0.95 1.35 1.6 2.2 2.6 3.5
2.9 2.2 0.95 1.1 1.7 2.0 2.9 3.3
2.0 2.0 0.75 0.8 1.5 1.6 3.2 3.2
(3.12) 1
A2/ Bl/ Azn B2,, A2, B2,, A2,
B2,,
2 ~
1.9 1.9 1.0 0.95 1.8 1.55 3.4 3.0
1.9 1.9 0.8 0.9 1.6 1.55 3.0 3.3
1
2
0.75 0.65 1.1 1.0 1.95 1.99
0.7
(b) Rotational Motion n
D,'
F2
b/ D, b, D, b,
CIZ(1)
C12(II)
2
4
m
0.6 0.9 1.0 1.35 1.6 2.1
0.65 0.75 1.0 1.15 1.8 2.0
0.7 0.7 1.0 1.1 1.9 1.9
A2,g
B2; Ain
B2, A2,
Bin
0.6 1.05 1.1 1.95 2.05
OD,, = ([2/7,)/2. bD, = (Ar;/7,l)/6. CDefined in (3.12). dDefined in (3.13). = (77;/7.)12/8. fD, = (A0n/7,L)p/24. g.* Defined in (5.7) and multiplied by (1/2)2; 7, A0 in radians and D in 1o - ~m2 s-I.
P ( ( r I L ) '1 rI1) = exp(-(rl - 0.1)/0.18)) D. The Approach to the Diffusion Constant. The relationships D = lim { ( x ( t )- ~ ( 0 ) ) ~ / 2=t lim ) ( ( l r ( t )- r ( 0 ) I 2 ) / 6 t ) I d = lim - - ( ( x ( t ) 4dt (-a,
The gradients of the MSD in the sequence of double steps generated by either criteria (3.1) or (3.7) are constant after n = 2; four steps are therefore sufficient to provide a fair estimate of the diffusion constant. This is surprising because 4 ( r 1 ) or 4 ( r I L ) are smaller than the decay time of the VAF. The agreement with DVAF is normally within the combined margin of error. E . The Diffusion Constant in Terms of Variances and Covariances of Amplitudes of Oscillation. The itinerant Cartesian coordinate oscillations constitute a stochastic process which is a "persistent" random walk in which positive and negative steps of variable lengths and duration alternate. The MSD after 2n steps is
=
(52:)
((,$I(')
- SI(')
+ ... - , $ 1 ( 2 n ) ) 2 )
(3.14)
Because of the stationary property ((t1("))2)
=
(,$,2)
and
(,$l(i),$l(i+j))
= (,$l(l)~lo'+l))
(3.14) can be written in the form
expect from the VAFs reported in ref 8. Distribution of Intervals. The probability densities f ( A r , ) for the lengths a n d f ( r l L ) for the duration of right-angle steps are similar for the 1- and 2-center LJ systems. Figure 3b showsf(Arl) for F2. In contrast to the monotonically decreasingf(Sl), f(Arl) rises from 0 to a peak at ~ 0 . (7Ar,) and has virtually vanished at =5 ( A r l )(=0.8). The difference between the distribution of quasi-amplitudes and right-angle steps arises from the fact that for the former, the initial velocity is 0, whereas for the latter the probability of a very small initial velocity is negligible. The distribution of the time intervals corresponding to eq 3.1 and 3.7 on the other hand is similar; f ( r l L )(not shown) has a maximum at 2/3( r l l ) followed by a smooth decay. The distribution function decreases after a short Gaussian portion exponentially:
f-m
(3.13)
t--
I d - ~ ( 0 ) )=~lim ) --(rl(t) t-m 12dt
- r(O)I2)
prompt the question as to the number of quasi-oscillation or right-angle turn steps in which the asymptotic limit is attained. The results are shown in Table IVa). On the left-hand side the ) n = 1, 2 and the asquotients ( , $ 2 2 / 7 2 n ) and ( A r 2 2 / 6 r 2 n for ymptotic value are given.
by introducing the variance varl =
-
(,$1)2
(3.15)
and the correlation coefficients Pj = ((€i(l)€i"+')) and going to the limit n
-+
D = varl(l - 2p,
m,
- (h2))/varl
(3.16)
one obtains 2p2 - ... + ...)/ 2 ( r 1 )
(3.17)
where ( ~ 2 has ~ ) been replaced by 2n( rl ). In practice, our data do not permit reliable estimations of the very small pj for j > 2. Successive approximations to the diffusion constant, Le., in which the pj are neglected after j = 0, 1, 2, respectively
Do= v a r l / 2 ( r l )
4 = varl(1 - 2~1)/2(71) D2 = varl(l - 2p1 + 2p2)/2(rl)
(3.18)
are given in Table V, which also lists values for var,, p , , p2, and D,. Except for the 1-center LJ systems at lower density and higher temperature, pI is positive and therefore reduces the diffusion constant; it is a measure of the quasi-elasticity of the cage which
Murthy and Singer
The Journal of Physical Chemistry, Vol. 91, No. 1, 1987
28
TABLE V: Step Variances, Correlation Coefficients, and Diffusion
1.0
(a) Translational Motion system
varl/10-2a2
pI
P2
Ar(I)2b
2.03 4.19 8.09 1.17 1.31 2.38
0.083 -0.0006 -0.012 0.111 0.073 0.028
-0.0033 0.00036 0.018 -0.036 -0.020 -0.011
Ar(II)b
Ar(III)b F2 CI2U)
CI,(II) system F2 C12(I) C12(II)
var,/rad2 0.117 0.080 0.148
DOa 2.6 4.6 7.4 1.4 2.3 3.9
Dl
D2
DVAF
2.2 4.6 7.6 1.07 1.95 3.7
2.2 4.6 7.9 1.0 1.85 3.6
1.95 4.7 7.9 0.9 1.55 3.35
(b) Rotational MotionC p1 p2 Dod D, 0.134 0.029 0.85 0.65 0.134 0.027 1.4 1.0 2.0 0.083 0,017 2.4
D2
Dmd
0.7 1.1 2.1
0.7 1.0 2.0
m2 SKI. From ref 2, Table 3, where p I and p2 are erroneously multiplied by IO2. rotational diffusion constant D values have been multiplied by 1*/4 so as to make them comparable in magnitude with the translational diffusion constants. dFrom data in ref 10.
,\
i
t
.9
‘7: .5
-/
-.l
-.3
1 I
1
0.1
accounts for the back-scattering. The change of sign nearer the critical point may be associated with the hydrodynamic enhancement through vortex b a ~ k - f l o w . ~ For the I-center L J systems, the values of D, are within 10% of D,; in the case of the 2-center LJ systems, D , and even D, overestimate the diffusion constant by 10-20%. If all correlation coefficients after p 1 were neglected, the expression for the difference quotient (eq 3.13) would be reduced to D,. In this approximation D is determined by three parameters: varl,pl, and (T,);these are obtained from statistics over the interval 2~~which is appreciably smaller than the decay time of the VAF. The alternating random walk or quasi-oscillation with correlation confined to adjacent steps appears to provide a reasonably adequate description of diffusion in 1-center and (to a lesser extent) 2-center LJ liquids. Some information about the shape of 2-step sections of the sequence of right-angle turns can be obtained from the values of (ArI2)and (Ar2,). Since (Ar,’) = 2(Ar12)+ 2 ( A r , ( ’ ) ~ A r , ( ~ ) ) and (cos y ) =
(Ar2,)
-
2(ArI2)
(Ar,)2 we obtain the following values for \k = cos-’ ((cos7 ) ) ;\k for Ar(1) = 99’; Ar(I1). 93’; Ar(III), 96’; F2, 92’; Cl2(I), 99’; C12(11), 100’. A semicircular orbit consisting of two coplanar successive right-angle turns (quarter-circles) would lead to \k = 90’. For two coplanar quarter-circles curved in opposite directions, \k = 0’. Although the steps are not coplanar or form part of a circle, the angle \k 90’ is a strong indication that there is no change of sign of the curvature, i.e., that successive steps form part of an irregular spiralling rather than “snaking” motion. Further details on the curvature of trajectories are given in ref 2. IV. Rotational Autocorrelation Functions
The most significant are the orientational AFs Ct(2) = P,(cos a ( ? ) ) where Pe(t) is the Legendre polynomial of order C, and cos a(r) = u(O).u(t); and the (normalized) angular momentum or velocity A F (AVAF) C,(t) = CJ(t)= ( w ( O ) . u ( t ) ) / (w(O),)
(4.2)
where w is the angular velocity vector. The corresponding angular acceleration or torque AFs are obtained by substitution of cb for w in (4.2). The relaxation of the direction of a linear molecule is described by
c,(t)= (u(o).u(t))
= (COS a(?))
(4.3)
I
I
I
I
0.5
0.9
, 1.3
I ‘5
TIME I ps
Figure 4. AF CI= PI ((u(0)-u(t))for gaseous and liquid CO modelled by Stockmayer or Stockmayer + quadrupole potenitals; from ref 26. These functions were first computed in MD simulations of liquid CO based on a Stockmayer potential (LJ + dipole) with and without an added quadrupole moment.26 Nonpolar diatomic liquids (e.g., N2, F2, C12,Br2) have subsequently been simulated by 2-center LJ “atom-atom’’ potentials sometimes modified by a quadrupole m ~ m e n t . ~ ~ , ~ ~ , ~ ~ - ~ ~ The rotational AFs are more sensitive to the interaction potential than the VAF. For almost spherically symmetrical potentials or at low density, the rotation is merely damped and C,(t) decays through one or more peaks corresponding to complete rotations. If the anisotropy and density are large, the hindered rotation degenerates to libration, and Cl(t)decays monotonically and, after an initial period, exponentially (see Figure 4). The decay time of C1in the latter case is several times larger than that of the VAF. The exponential behavior can be accounted for in terms of a sequence of small uncorrelated steps. This is the basis of the Debye theory according to which Ce(t) = exp(-D,C(P + l)?) where D, is the rotation diffusion constant. The AVAF shares the major features with the VAF: reversals of the COM velocity correspond to reversals of the angular velocity; if predominant, they produce a negative lobe after T ~ .T~ can be loosely interpreted as of the period of libration. If the molecular anisotropy is weak, or the density low, the decay is monotonic.* The long tail of the VAF is absent. To account for the mechanism of rotational diffusion a number of simple models have been p r o p o ~ e d . ~The ~ . ~most ~ successful of these is perhaps the “J-diffusion” model,31according to which free rotations through random angles are separated by collisions in which the angular momentum is thermally randomized. Neither this nor any of the other known simple models can be made to fit the MD simulation results for dense diatomic liquids.28 More elaborate models suffer from the presence of parameters which cannot readily be related to physical proper tie^.^^,^^ (25) Singer, K.; Taylor, A. J.; Singer, J. V. L. Mol. Phys. 1977, 33, 1757. (26) Berne, B J.; Harp, G. D. Adu. Chem. Phys. 1970, 17, 63. (27) Cheung, P. S. Y.; Powles, J. G. Mol. Phys. 1975, 30, 921. (28) Steele, W. A.; Streett, W. B. Mol. Phys. 1980, 39, 279. (29) Murthy, C. S.; Singer, K.; Vallauri, R. Mol. Phys. 1983, 49, 803. (30) Invanov, F. N. Sou. Phys. JETP (Engl. Transl.) 1964, 18, 1041. (31) Gordon, R. G. J. Chem. Phys. 1966, 44, 1830. (32) Coffey, W. T. Mol. Phys. 1979, 37, 473. (33) Coffey, W. T.; Evans, M. W.; Evans, G.J. Mol. Phys. 1979,38,477. ( 3 4 ) Berne, B. J.; Montgomery, J. A. Mol. Phys. 1976, 32, 363. (35) Steele, W. A.; Streett, W. B. Mol. Phys. 1980, 39, 299. (36) Evans, M. W.; Ferrario, M.; Grigolini, P. Mol. Phys. 1980, 39, 1369. (37) Evans, M. W.; Grigolini, P. Mol. Phys. 1980, 39, 1391. (38) Ryckaert, J. P.; Bellemans, A,; Ciccotti, G. Mol. Phys. 1981, 44, 979.
Feature Article
The Journal of Physical Chemistry, Vol. 91, No. I , 1987 29
A review of the extensive literature on translation-rotation coupling (T-R) is beyond our scope. Some aspects are discussed in section IIC and in ref 35 and 38. The presence of a correlation between Pt (cos a(?))with functions of the COM displacement r ( t ) - r(0) has been established and analyzed for some model systems.
the COM diffusion with the rotational diffusion of the atoms. The difference quotients derived from both types of MSDs after two and four steps give reasonable estimates of the rotational diffusion constant. D. The Rotational Diffusion Constant in Terms of Variances and Covariances of the Librational Amplitudes. The arguments which lead to (3.14)-(3.16) are also applicable to quasi-librational steps leading to successive approximations to D,. The results are shown in Table Vb. It can be seen that the approximation D,and the corresponding 3-parameter stochastic model work very well.
34,36337
V. Analysis of Roational Trajectories A. Quasi-librations about a Fixed Axis. Just as the criterion u, = 0 serves to subdivide the COM trajectories, so the criterion u,
=0
(5.1)
divides the rotation about the x axis into alternately opposite librational amplitudes or steps, referred to as 7. Alternatively, one may use the criterion zi, = (w x u), = 0 (5.2) Since IwI = lul, these criteria are equivalent as far as the contributions of the mean squared step amplitudes and hence the diffusion are concerned. Equations 3.3, 3.4, and 3.5 now become v2(m) (Itn(m))
=
= Iv1(m) - vl(m+l)l (vn),
(5.3)
(~n"') = (7,)
(5.4)
and the path length on the unit sphere 0, =
Im+l;l -1vx2 + vy2 + q,211/2dt,
(7, =
The rotational analogue is
vl)
-
=0 (5.6) determines time intervals ( T , ~ and ) angular intervals ( 0 ) during which the trajectory on the unit sphere turns through 90'. The mean values are given in Table IIIb. The time intervals are a little shorter and the angular intervals longer than the values of the Cartesian component librations in Table IIb. The path lengths are only 10-12% greater than the distance (on the unit sphere) between end points; the angle traversed in two successive steps is =looo. C. Corwergence to the Diffusion Constant. As in section IIIC we examine the convergence of the ratios W(O)*O(T,')
( v 2 , 2 / ~ 2/ 2n )(libration) and (62,2/~2,)/6 (right-angle turn) to the asymptotic value D, and compare the estimate of the rotational diffusion constant so obtained with the difference quotients (v2n+z2
1 (Oi2V) = -(Ioiq2) kT
(5.5)
The mean values for v l , q2, T , , and O1 are shown in Table IIb. The last three columns contain the mean COM coordinate displacements which occur during the intervals generated by the criterion (5.1). The results are in several respects similar to those given in Table IIa for the translational trajectories. The mean angular amplitudes of 15-19' correspond to atomic displacements of 0.08-0.1 (=0.23-0.27 A). The path lengths are about twice as large. The relative fluctuations are 0.7-1.0. Accelerations change sign 1.7 times more frequently than velocities and the corresponding amplitudes have the ratio (vl(acc))/(v,(vel)) z 1/2. The remarks concerning anharmonicity [section IIIA] therefore apply to librations as well as to oscillations. The similarity between librations and oscillations in our systems extends to the distributions of both time and coordinate intervals. The distributions f(vI) and P ( v , ' 1 q 1 ) are similar to Atl) and P({,'IF,) shown in Figure 3b: the decay is a little faster than exponential, and value of vl > 4 (7,)are rare. B. Right-Angle Turns. The criterion
Azn =
VI. Energy-Related Characteristics A. Curvature of the Potential Field. Assumptions about the local potential field form the basis of the cell or free-volume theories of liquids. An important characteristic which also determines the isotope separation factor and the first-order quantum connection of the free energy in the Wigner-Kirkwood expansion is the mean Laplacean of the 1-particle potential energy, or equivalently, the mean-squared force 1 (Oi2v) = ,,(lV,V1*, (6.1)
- 7122)
4(71)
where Oi is the rotational gradient operator. Equations 6.1 and 6.2 determine the mean translational and rotational Einstein frequencies wEt
= (V?V/m)'/z,
These entries in Table IVb have been multiplied by (1/2)2 and the angles measured in radians so as to facilitate comparison of
WE'
= (OL2V/I)'/2
where I is the moment of inertia. The histogram of curvature f(V:V) for Ar(1) consists entirely of positive entries; a small number of negative values occur at lower density and higher temperature (Ar(II1)). The curves are shown in Figure 8 of ref 2. Figure Sa shows the distribution of f(uEf) andflu$) for F2. The curves are similar to those of Ar, but the histogram for uE' contains a small fraction of imaginary entries corresponding to unstable librations. Although the quasi-oscillations and -librations are very anharmonic, the frequencies obtained directly from the cycle periods 2 ( ~ ] are ) reasonably close to the values derived from the Laplaceans (eq 6.1) Ar(1)
Ar(I1)
7.0 7.7
5.9 6.9
TI), THz wE[from (6,1)], THz
Ar(II1) 5.0 6.1
F2(tr) F2(rot) 9.2 9.8
9.2 9.5
The predominance of positive Laplaceans implies that in the dense liquid almost all molecules at any time sit in bowl-shaped potential wells. This appears to contradict the finding that the accelerations change sign 1.5-2 times more frequently than the velocities, which suggests that potential barriers are being crossed. If this were so, the particles should be frequently near the top of potential energy hills where V t V or 0;V are negative. This paradox is produced by the false model of an approximately static local field. In reality the force field acting on a particle changes because of its own motion and that of its neighbors. The rate of change of the force on particle i is given by d -Fi = -x(ki - kj).ViViV(rij)= - ?j)-V,F,j (6.3) dt 1 J where ri.EViFijand -Erj.ViFij j
J
(5.7)
(6.2)
are respectively the contributions of the motion of the ith and the other particles to the rate of change. We have calculated the averages DF; = ( lri-VixFijl)and DF? = ( ICrj.ViFijI) J
J
,
,
Murthy and Singer
30 The Journal of Physical Chemistry, Vol. 91, No. 1, 1987 WL
4;;"
I THz
7;6
19.4
13;2
more rapidly than AKE, but its long thin tail persists well beyond the decay of AKE,. The inset of Figure 5b shows that the CFs AKE - APE and AKE, - AKE, are positive, though the latter rather is weakly so. This indicates that the exchange of energy with the neighbors masks the interconversion of P E KE, and KE, KE,, since the corresponding CFs would otherwise be negative.
-
-
VII. Conclusions c
-$ "
2.3
9.5
5.9
13.1
16.7
r
u E lTHz
TIME1 Pl
i 5. (a) Unnormalized probability density&$) (dashes) and flu$) (continuous curve) for F2 calculated from the Laplaceans of the potential. Ordinate scale: arbitrary units. (b) AFs of APE (full curve), AKE, (dashes),AKE,(dots) for F2.Inset: CFs of APE - AKE (full curve) and AKE, = AKE,- AKE, (dashes). F
For Ar(I), -(II), and -(III) the ratios DFp/DFt are 0.85, 0.92, and 0.96; this shows that the motion of the "cage" contributes almost as much as the self-motion to the rate of change. The frequent accelerations and retardations are mostly caused by changes in position and depth of the potential well, rather than by inversions of its curvature. B. Fluctuations of the Kinetic, Potential, and Total Energy. A plot of the changes in KE, PE, and T E (kinetic, potential, and total energy) over 4 ps of Ar Iz shows them to be rather irregular; the KE fluctuates more rapidly and the fluctuations in T E are quite large. In the case of diatomics the translational and rotational kinetic energies (KE,, KE,) are distinguished. The normalized A F of the fluctuation AA and the C F between AA and AB ( A A = A - ( A ) , AB = B - ( B ) ) are defined by
(AA(O)Wt))
and
(A"w) (( LL-W)~) (AW)2))1/2
Some of these functions for F2are plotted in Figure 5b. The decay of the AFs is monotnic; that of APE is much slower than those of AKE, and AKE,, largely on account of a sharp change in slope between 0.15 and 0.2 ps; up to -0.4 ps, AKE, decays
1. Translational and rotational trajectories of molecules in 1and 2-center Lennard-Jones liquids show no features which point to a simple physical model of the motion. 2. The mean amplitudes of the quasi-oscillation and -librations generated by projection of the translational and rotational trajectories onto a fixed axis are by almost an order of magnitude smaller than the atomic diameters or interparticle spacings and decrease with the molecular anisotropy. This also applies to the steps generated by right-angle turns of the translational and angular velocity vectors. The smallness of these steps and direct examination of the nearest-neighbor lists show that clusters of molecule persist over periods which are multiples of the relaxation time of the velocity autocorrelation function. 3. The distributions of the above-mentioned amplitudes decrease monotonically, rather like the distribution of free paths in a hard-sphere fluid.3 This confirms that in the LJ liquids, as in HS fluids, diffusion does not occur through an alternation of small oscillations and "large" jumps. 4. The fact that translational and angular accelerations change sign 1.5-2 times more frequently than the corresponding velocities shows that the quasi-oscillations and librations are anharmonic. 5. The quasi-oscillations and -librations are equivalent to an alternation of positive and negative steps of variable length and duration, constituting a "persistent" random walk in which the correlation between the lengths of steps rapidly decreases with their separation (in the sequence). If only correlations between adjacent steps are taken into account, one obtains a stochastic model which is fairly adequate for most of the systems considered. The diffusion constant is determined by three parameters: the mean duration and variance of the step lengths, and the correlation coefficient between adjacent steps. This information is obtained by averaging over time intervals which are shorter than the decay time of the VAF. 6. For the systems considered, the translational and, to a lesser extent, the rotational Laplaceans of the potential energy are predominantly positive; almost all potential wells acting on the particles are bowl-shaped. 7. The cages confining molecules in the liquid are highly mobile; the rate of change of the force acting on a particle is produced to roughly equal extents by the effect of the motion of the particle in a fixed cage and the effect of the motion of the cage on the fixed particle. 8. Energy fluctuations show no surprising features; the interconversion of potential and kinetic energy is small compared with the fluctuations of the total energy. The cross-correlations between kinetic and total energy and that between rotational and translational kinetic energy are therefore positive, and the latter is very weak.
Acknowledgment. This work has been supported by a grant from the Scientific and Engineering Research Council (U.K.). C.S.M. thanks Professor Peter J. Rossky for his encouragement and Dr. Jiirgen Schnitker for the use of his stereoscopic plotting program.
