-
-
J. Phys. Chem. 1989, 93,6975-6979
with L finite, in the limit g O+ the liquid phase could break up into multiple domains by the creation of more than one interface, in order to take advantage of the entropy gain associated with the increased wandering of the domains. (A similar argument explains the lack of long-ranged order in the one-dimensional Ising model.) Alternatively, consider a canonical system with a single liquid domain containing two interfaces in an infinitely long column with a symmetric potential 40(z) = 40(-z) that tends to locate the center of the domain near z = 0. Clearly, this system also develops diverging zero wave vector fluctuations as the field strength tends to zero. We now discuss some of the additional complications that one encounters in an analysis of a system with finite L and L,. The classical picture is most naturally realized in a finite system with short-ranged wall potentials at z = fL,/2 chosen to induce phase ~ e p a r a t i o n .It~ suffices to impose a “hard wall” condition d0(z) = m for lzl > L , / 2 and to incorporate into 40(z)for z near the lower wall at z = -L,/2 an additional short-ranged attractive term that favors wetting by the denser (liquid) phase. Otherwise, as before, we take 40(z) = mgz with g arbitrarily small. For such a finite system with L, and L of O(&,),a nearly classical picture holds provided we use a canonical ensemble description: where the number of particles, N, is fixed and chosen so that the Gibbs dividing surface is at z = 0, consistent with the potential 40(z). With this setup, we expect to find an “intrinsic” interface whose properties are essentially unchanged as g O+,
-
6975
Le., as the potential becomes arbitrarily weak in the interfacial region. In principle, with L, finite there are contributions in (4.5) from ~ p(z) vanishes the top and bottom boundaries of our ~ y s t e m .Since for Izl > LJ2, the following boundary terms are generated by (4.5): C0(-Lz/2,zz)plw- C0(L,/2,z2)pw. Here Co(zl,zz)denotes the integral over rz of C(z1,z2,rI2) and plw and pw are the densities at the bottom and top walls. Using the definition (4.3), we expect that a given density change near the wall should correspond to a potential essentially localized near the wall, since the wall damps out long-ranged correlations. Thus the boundary corrections should be negligible for z2 near 0 in the interfacial region. The same “shift” argument used to derive (4.5) can be repeated to give a more general derivation of (3.8), valid for finite systems and in the canonical as well as the grand ensemble. In the grand ensemble, with L, infinite and L of O(&,),(3.8) holds true, not because of long-ranged “horizontal” correlations, but because of long-ranged “vertical” correlations caused by fluctuations in N as g 0.’ In the canonical ensemble, with L, and L of O(fb), boundary corrections from the top and bottom walls, yielding terms like H0(*L,/2, z2), allow (3.8) to continue to hold as g - 0.’ Note that because of conservation requirements in the canonical ensemble and the small system size, the density displacement induced by a change in potential near the wall could have substantial effects on the density profile in the interfacial region. Thus, boundary corrections are more important in (3.8) than in (4.5). -+
Measures of Effective Ergodic Convergence in Liquids Raymond D. Mountain Thermophysics Division, National Institute of Standards and Technology,+ Gaithersburg, Maryland 20899
and D. Tbirumalai* Department of Chemistry and Biochemistry and Institute f o r Physical Science and Technology, University of Maryland, College Park, Maryland 20742 (Received: March 1 , 1989; In Final Form: May 4, 1989)
The recently introduced measure for ergodic convergence is used to illustrate the time scales needed for effective ergodicity to be obtained in various liquids. The cases considered are binary mixtures of soft spheres, two-component Lennard-Jones systems, and liquid water. It is shown that various measures obey a dynamical scaling law which is characterized by a single parameter, namely, a novel diffusion constant. The time scales for ergodic behavior are found to be dependent on the particular observable being considered. For example, in water, the diffusion constants for the translational and rotational kinetic energies and for the laboratory frame dipole moments are very different. The implications of these results for the calculation of the dielectric constant of polar liquids by computer simulations are discussed.
I. Introduction The ergodic hypothesis is one of the central concepts in the development of the statistical mechanical theory of thermodynamics.’ This hypothesis states that time averages and ensemble averages are identical for a system in thermodynamic equilibrium. The ergodic hypothesis provides a connection between the results obtained from analyzing the trajectory of a many-particle system, generated in computer simulations of the dynamics of molecular motions, to the results calculated by using the more abstract notion of an ensemble. An ensemble characterizes the various equilibrium states of the system. In computer simulation studies, estimates of these observable properties are obtained as time averages of functions, called phase-space functions, of the coordinates and momenta of the particles making up the system. Letf(t) be a phase-space function whose time average is a physically observable quantity. The ergodic hypothesis asserts that Formerly the National Bureau of Standards.
0022-3654/89/2093-6975$01 SO10
r--? lim L l r od t f ( t ) =
v)
(1.la)
where the ( ) indicate an ensemble average. The appropriate ensemble for constant-energy, constant-volume molecular dynamics simulations is the microcanonical ensemble. For a system with a given Hamiltonian H(p,q) and total energy Eo,averages are expressed as
cf) =
p r w m q ) - Eo)f(r)/Sdr
W ( p , q )- Eo)
(l.lb) and the integration is over the phase space r of the system. It should be emphasized that there does not appear to be any uniform definition of ergodicity in the literature.2 As was ori(1) A good discussion of ergodic theory in statistical mechanics is provided by: Wightman, A. S. In Statistical Mechanics at the rum of the Decade; Cohen, E. G . D., Ed.; Mercel Decker: New York, 1971; pp 1-32.
0 1989 American Chemical Society
6976 The Journal of Physical Chemistry, Vol. 93, No. 19, 1989 ginally used by Boltzmann, ergodicity implied that the trajectory of a many-particle system passes through every allowed point in phase space, regardless of the initial condition chosen. This condition for ergodicity is too strict and is seldom if ever used in practice. It is well-known that the time required for a N-particle system to cover all points in the allowed region of phase space scales as exp(N), and thus ergodicity in the sense used by Boltzmann is not of interest. Because the observation time T , , ~ is always much less than exp(N), the question of real interest is whether on time scales comparable to 7,p the trajectory effectively samples the allowed region of phase space. For effective ergodicity to be obtained we imagine that the allowed phase space is divided into several smaller regions. This coarse graining of phase space means that one is looking at the properties of the trajectory over a finite part of the phase space. If the trajectory uniformly samples the coarse-grained phase space, then the system is considered to be effectively ergodic. Notice that our use of effective ergodicity is similar to the concept of mixing. The analysis of the mixing property (the spreading from a small finite region of phase space into the entire allowed region under time evolution) for a generic Hamiltonian system is an unsolved problem. Recently, we introduced a dynamic measure to probe the minimum time scales required for achieving effective ergodicity in supercooled liquids and in g l a s ~ e s . In ~ this paper we use this measure to investigate the minimum time T needed for effective ergodicity to be obtained in various liquid systems. It is well-known that ergodic behavior is realized for systems with the property that any two trajectories in phase space, which are initially close, diverge exponentially as time advances. This mixing property is sufficient to ensure ergodic behavior. It is not known whether this is also a necessary condition for ergodic behavior. It is possible that the structure of phase space can be such that one has a complex, multivalley structure. The various valleys can be separated by bottlenecks. In such cases (as often happens in supercooled liquids and glasses) one has islands of stochasticity; Le., two nearby trajectories will diverge with time only if they belong to the same valley. In stating this, we are assuming that T is smaller than the time scales needed to overcome the bottleneck in phase space. In liquids, on the other hand, the energy Eo is usually higher than the largest barriers encountered in the multidimensional potential energy picture. For this case the stochasticity can be inferred if certain global stability criteria are met. These criteria, expressed in terms of Lyapunov exponents which indicate the rate of growth of the distance between trajectories, are computationally more difficult to determine than are the solutions of the equations of motion.s As a result, it has generally not been possible to make quantitative statements about times scales for ergodic convergence in condensed matter systems. The method introduced by Thirumalai et al. in connection with the structural glass problem3 is useful in assessing the time scales needed for effective ergodicity to be established. They defined a measure of the mixing property of two trajectories in phase space which is relatively easy to compute during a molecular dynamics simulation. This measure was specifically designed for Hamiltonian systems with a large number of degrees of freedom. Mountain and Thirumalai6 have subsequently shown that in the absence of multivalley structure it is sufficient to examine the properties of a single trajectory in order to determine whether it has ergodic behavior. This has the advantage of removing the need to examine two separate trajectories in the phase space of the system of interest. The guiding principle of this approach is that of the statistical equivalence of the atoms or molecules which make up a fluid in thermal equilibrium. (2) See for example: Ma, S.K. In Sratisticol Mechanics;World Scientific Publishing: Singapore, 1985; Chapter 26. ( 3 ) Thirumalai, D.; Mountain, R. D.; Kirkpatrick, T. R . Phys. Reu. A 1989 -.- - , 39 - . , -7567 - -- . (4) See for example: Regular and Srochusric Motion; Lichtenberg, A. J., Leiberman, M. A., Eds.; Springer-Verlag: New York, 1983. ( 5 ) Posch, H. A,; Hoover, W. G . Phys. Rev. A 1988, 38, 473. (6) Mountain, R. D.; Thirumalai, D. J. Chem. Phys., in press.
Mountain and Thirumalai
In this paper, the ergodic convergence properties of three fluid systems are examined by using this new approach. These systems are a binary mixture of soft spheres, a binary mixture of Lennard-Jones atoms, and liquid water. The ideas needed to determine ergodic convergence are discussed in section 11, and the explicit quantities needed are introduced. The systems examined are discussed in section 111 along with the details of the simulations. Our results are discussed in section IV and are summarized in section V. 11. Ergodic Convergence Concepts The essential idea needed to develop a quantitative measure ~ refers of ergodic convergence is called statistical ~ y m m e t r y .This to the condition that, for a fluid in thermal equilibrium, all particles of a given type have identical average characteristics. This equivalence is implicit in the statement that an equilibrium, onephase fluid is homogeneous. For example, the average kinetic energy of a particle is (3/2)kBT,the product of Boltzmann’s constant kB and the temperature T , independent of where the particle happens to be located in the fluid. Other properties, such as the average energy per particle, may depend on the type of particle but not on the location of the particle. When these conditions are satisfied for time averages over an interval, T,we say that ergodic convergence has been achieved. The question is now how to make these remarks quantitative. Suppose that we are interested in a property,f, which can be obtained as a suitable average of a phase function F. Let us introduce two time-averaged quantities (2.la) which is the time average of property fassociated with particle j , and
fit) =
1N
G,yj(t)
(2.lb)
which is the average of property f taken over all particles. A measure of the ergodic convergence for property f is a metric for
f
-
If the system is ergodic, all of the&(?)%approach the common value f i t ) for long times and Okf) 0. The vanishing of the metric for long times is the indication that the system is ergodic. Thirumalai and Mountain3 used an energy metric, O,, where e is the energy associated with an individual atom in the fluid. They found that the time dependence of O,(t) could be simply represented as Oe(t)/Qe(o)
1 :
1/DJ
(2.3)
for times greater than a short initial time interval. The energyspace diffusion coefficient De is a function of the density, temperature, etc., of the fluid. It should also be possible to use other properties, such as the kinetic energy or various components of the stress tensor, which can be separated into contributions from individual atoms or molecules, as the basis for a measure of effective ergodic convergence. It is not known under which circumstances these measures are equivalent indications of ergodic behavior. In the following, we examine the kinetic energy and the total energy of simple mixtures as propertiesf. We also examine the translational and rotational kinetic energies of water molecules as test propertiesf. Some preliminary results on the use of the laboratory frame components of the dipole moment of individual water molecules will also be presented. 111. Simulation Details The computer simulation technique of molecular dynamics has been used to generate the necessary information to evaluate various measures of effective ergodic convergence. The methods are standard and are discussed in more detail elsewhere.’**
The Journal of Physical Chemistry, Vol. 93, No. 19, 1989 6977
Measures of Effective Ergodic Convergence in Liquids Three separate systems have been examined. The first is the soft-sphere binary studied by Thirumalai and Mountain in their work on the structural glass problem.' Here we consider an equimolar mixture of spheres interacting via the inverse power potential d)Adrij)
= c(uAB/rij)12
I
(3.1)
where the subscripts i, j and A , B = 1 , 2 refer to a particle of species A at position ri and to a particle of species B at position rj and rij = Iri - ,'I. The diameters were taken to be u I I = 1 , uI2 = 1.05, and uZ2= 1 . 1 . The mass of type one particles was taken to be ml = 1 , and the mass of type two particles was taken to be 2 m l . The system consisted of 500 particles in a periodic cubic cell such that the density was 0.85 p a r t i ~ l e s / ( o , , )and ~ the temperature was 1.16c/kB. The second system is an equimolar Lennard-Jones fluid mixture interacting via the pair potentials +Adrij) = 4 ~ [ ( ~ A B / r i-j )('~~A e / r i j ) ~ I
6oo 500
t Figure 1. The inverse total energy metric (solid line) and the inverse kinetic energy metric (dashed line) for the Lennard-Jones mixture. The universal form is evident. The time is in units of T .
(3.2)
where we have used the same values for the diameters, the density, and masses as were used for the soft spheres. The temperature was 2.02c/kB. The third system in water interacting via the SPC (simple point charge) model. This interaction is described fully in ref 9 and 10. The system consisted of 108 molecules at a density of 999 kg/m3 and a temperature of 3 1 2 K. In each case, the simulations started from well-equilibrated configurations. Constant-volume, constant-energy conditions were imposed. The equations of motion were integrated by using the Beeman algorithm." In the case of water, an iterated version of this algorithm was used.* The time steps were such that the energy was conserved to within 1 part in 5000 during the run with no adjustment. When considering a mixture, it is necessary to slightly modify the definition of Qkt) since the long time values ofjJt) can depend on the species. In practice we used
. ^ a, v
c
0.0
-3
where the subscript a refers to the species of particle j and (3.4)
is the average over the Nu particles of type a. In this way we ensure that O L t ) goes to zero at long time for an ergodic system. Three quantities have been used as the quantity$ These are the kinetic energy of a particle, kj = (1/2)mju?, the energy of a particle
and paj, the a t h component of the dipole moment of particle j . For water, we also examine the rotational kinetic energy
1
-1
1
ENERGY Figure 2. The distribution of time-averaged energy of the particles, e,, is displayed on the left-hand side of the figure and the time-averaged kinetic energies, kj, are displayed on the right-hand side of the figure. The averaging time if 107, and the energies are in units of B introduced in eq 3.2. The system is a mixture of Lennard-Jones particles. The circles joined by solid lines represent results for the type 1 particles, and the triangles joined by dashed lines represent results for the type 2 particles.
ergodic. In our previous work3 we had argued that the two-trajectory energy metric, d ( t ) which is defined in ref 3, is governed by a single parameter referred to as the energy-space diffusion constant. The plots of n(O)/Q(t)as a function of time reveal the same scaling form as does the energy metric, i.e. Q(t)/Q(O) = l/Det
where Zju and w , ~are the ath principal moment of inertia and the a t h component of angular velocity of molecule j in the molecule-fixed reference frame. IV. Results A . Dynamic Behauior of Q ( t ) . In all cases, Q ( t ) decays smoothly to zero, indicating that the systems become effectively (7) Thirumalai, D.; Mountain, R. D. J . Phys. C 1987, 20, L399. Mountain, R. D.; Thirumalai, D. Phys. Reu. A 1987, 36, 3300. (8) Mountain, R. D. J . Chem. Phys. 1989, 90, 1866. (9) Berendsen, H. J. C.; Potsma, J. P. M.; van Gunsteren, W. F.; Hermans, J . In Infermoleculor Forces; Pullman, B., Ed.; Reidel: Dordrecht, Holland, 1981. (10) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J . Chem. Phys. 1983, 79, 926. (1 1 ) Beeman, D. J . Compuf.Phys. 1976, 20, 130.
3
(4.1)
for t greater than a short transient time. The slope taken from plots of Q(O)/Q(t)as a function of time can be used to obtain the energy-space diffusion coefficient De (or for other quantities, the appropriate-space diffusion coefficient). The scaling behavior exhibited in eq 4.1 suggests that the variables constituting the metric Q ( t )behave as normally distributed random variables. The behavior in eq 4.1 indicates that the variance decreases as l / t . It is also instructive to display the distribution of the e,'s at various times as it illustrates how the statistical symmmetry of the particles is established. From an examination of eq 4.1 and the distribution of time-averaged energies, ej(t), an estimate of the time scales for ergodic behavior for the quantities can be inferred. The results for the soft-sphere mixture are examined in some detail elsewhere? Here we concentrate on the Lennard-Jones mixture and present preliminary results for water. B. Lennard-Jones Mixture. The energy and kinetic energy metrics for the Lennard-Jones mixture are displayed in Figure
6978
The Journal of Physical Chemistry, Vol. 93, No. 19, 1989
Mountain and Thirumalai
0.2.
100 r
,
q
80
~
+
z
_
1
I
,
d
60 10 T
J
3
10
0
5
20
ENERGY Figure 3. Distributions for the Lennard-Jones mixture at SOT. The symbols are the same as in Figure 2.
0 3.
10
50
60
'70
t , ps Figure 5. Reciprocals of the metrics for the translational kinetic energy of S P C water (squares) and for the rotational kinetic energy (circles) as a function of time in picoseconds. The lines are only a guide to the eye. The universal scaling behavior represents an average trend of the curves.
100
t
0.21 c
30
80
l
~i ._ 0
h
4
Y
v
n..
60
c:
v
\
0
h
E c 401 20 1
0
ENERGY
0
10
0
Figure 4. Distributions for the Lennard-Jones mixture at 1007. The symbols are the same as in Figure 2.
1 . From this plot we determine that De = 4.5 and Dk = 6.3 in units of 7-l. If one adopts the criterion that ergodic convergence is obtained when n(t)/n(O) 0.01, then the time needed to satisfy the ergodic hypothesis is T 207-307. To see that this is a reasonable criterion, we examine the distributions of the k,(t) and the e,(?) for t = 107, t = 507,and t = 100s in Figures 2-4. At 10s the distribution of the e i s is not well resolved while at 507 it is. Extending the averaging time to 1007 has little effect on the distributions other than to sharpen them, a result in accord with the l / t decay of O ( t ) . It is important to realize that the distribution for the total energies is bimodal for this mixture while the distribution for the kinetic energies has a single peak. The two-peak structure is a reflection of the different environments experienced by the two types of particles. Thus, we see that quantities,f, which are functions of the coordinates of the particles can provide more information about the state of the system than can the kinetic energy. A single-peak distribution for the k,(t), which narrows with increasing time, is expected for all systems with conservative forces, whether the system is in equilibrium or is in a metastable, but stationary, state such as a glass. Quite similar results were obtained for the soft-sphere fluid state. For that system we obtain a value of De = 3.97-' and a value of Dk = 3.57-I. Ergodic convergence was achieved for 'T = 2 0 ~ 3 0 7 . The distributions also showed behavior much like those of the Lennard-Jones mixture. C. Liquid Water. In this preliminary investigation we have posed the question of how long an interval is needed to establish the equality of the translational kinetic energy averages, k,(t), and the corresponding rotational kinetic energy averages, K j ( f ) , a condition which would be needed if the mean-square fluctuations in these quantities are to be reliably estimated.
-
20
30
40
50
60
t , ps Figure 6. The reciprocal of the translation-rotation metric for SPC water as a function of time.
'j -~
" 0
10 0'
30 10 50 60 'io t . ps
Figure 7. Reciprocals of the metrics for the components of the dipole moments of S P C water molecules as a function of time. These components are determined in the laboratory reference frame.
The metrics for the translational and rotational kinetic energies are displayed in Figure 5 as functions of time. From the corresponding scaling behavior, the energy-space diffusion coefficients are determined to be Dk = 1.0 ps-' and D, = 1.4ps-l. We also have computed a translation-rotation metric 1N n k , ( t ) = G z [ k j ( l )- K j ( t ) I 2 (4.2) is shown in Figure 6. The energy-space diffusion coefficient for this metric has a value 1.4 ps-l.
I
Measures of Effective Ergodic Convergence in Liquids
The Journal of Physical Chemistry, Vol. 93, No. 19, 1989 6979
The metric using the laboratory-frame components of the dipole moment of individual molecules reflects this. This is shown in Figure 7. The diffusion coefficient in the dipole moment space has a value of approximately 0.2 ps-'. The relatively smaller value of the coefficient measuring the rate of ergodic convergence of the dipole moment components is an indication that the sampling of rotational configurations proceeds slowly in water. This is due to the presence of hydrogen bonding. These relatively long lived bonds make structural rearrangements a slow process.12 This value has implications in the calculation of the dielectric constant using molecular dynamics. If one wants a reliable estimate for the dielectric constant, the sampling of the net dipole moment of the system has to be adequate.13 Using the diffusion constant of 0.2 ps-l for the exploration of the various dipole moment fluctuations, one can infer that the simulations have to be carried out for times on the order of 200-500 ps.14
considered) only if restricted averages are considered (restricted averages are obtained when the phase space is restricted to configurations belonging to a single valley).15 These conclusions have been found in several numerical studies of nonlinear coupled oscillators such as the Fermi-Pasta-Ulam systems.16 For ergodic systems the long time behavior of Q(t) given by eq 4.1 may be obtained by noting that Q ( t ) effectively measures the fluctuations offat time t.17,1*The measure Q(t) can be written as
(5.lb) where fifi(t) = fi(0- A t )
V. Conclusions In this paper we have introduced measures that can provide an estimate for time scales needed for ergodic behavior to be obtained in liquids. The applications to various systems have revealed several interesting features that should be useful in numerical simulations of the dynamics of liquids. The major implications of our study are summarized here. On the basis of numerical studies of the ergodic convergence of several physical quantities, we have suggested that the relaxation of the metric can be adequately represented by the universal scaling function given by eq 4.1. This behavior is observed after a quite short initial transient time. The minimum time scale, 7, for effective ergodicity to be obtained is given in terms of a single parameter, namely, the appropriate diffusion constant. Thus, after a relatively short simulation, De can be computed easily and hence 7 can be estimated trivially from the scaling behavior of the metric. The universal scaling behavior l / t is related to a generalized central limit theorem indicating the reduction in variance for independent random processes. It is also apparent that different physical quantities become effectively ergodic on different time scales. This is clear from the differences in the effective diffusion coefficients calculated for liquid water. Conversely, one could have situations in which certain quantities may appear ergodic on the observation time scale while others may be nonergodic on the same time scale. This is definitely the case in glassy systems or in general in systems which have multivalley free energy structures. In these situations the system appearieffectiveliergodic (depending on the quantity (12) Rapaport, D. C. Mol. Phys. 1983, 50, 1151. (13) Neumann, M. J. Chem. Phys. 1983,82, 1567. (14) Anderson, J.;Ullo, J. J.; Yip, S.J. Chem. Phys. 1987,87, 1726. They found that 100 ps was the minimum time required for good convergence.
(5.1~)
and 1N
C(s,-sz) = -Cfifi(s1) Nj-1 fifi(s2)
(5.ld)
The right-hand side of eq 5.ld is roughly equivalent to the result obtained by performing an ensemble average. Using the weak law of large numbers, one can identify C(sl-s2)as the equilibrium time correlation function. With this it follows that eq 5.la for long times becomes Q(t) = l / D / (5.2a) where
[
D, = lim 1 / 2 fC(s) ds] I---
0
(5.2b)
Notice that the equivalence between ensemble and time averages has been used to obtain eq 5.2. However, the measure Q(t) is a useful diagnostic even when ergodicity is broken. In fact, the behavior of n(t) different from that given by eq 5.2a can be used as a signature for the failure of erg0di~ity.l~ Acknowledgment. It is a pleasure to dedicate this paper to R. Zwanzig, whose advice and encouragement has been invaluable to us. This work was supported in part by National Science Foundation Grant CHE-86-57356. Registry No. Water, 7732-18-5. (1 5) Palmer, R. J. Adu. Phys. 1982, 31, 669. (16) Levi, R.; Pettini, M.; Ruffo, S.;Vulpiani, A. J . Slur. Phys. 1987, 48, 539, and references therein. (17) We thank R. Zwanzig for pointing this out to us. (18) See also: Zwanzig, R.; Ailawadi, N . K. Phys. Rev. 1969, 282, 280. (19) Thirumalai, D.; Mountain, R. D. J . Stu?.Phys., in press.
