J . Phys. Chem. 1993,97, 10478-10485
10478
Dipole Time Correlation Functions of the Stockmayer Fluid in the Microcanonical and Canonical Ensembles Vilia Ann Payne, Maria Forsyth,? Jiri Kolafa,* and Mark A. Ratner. Department of Chemistry and Materials Research Center, Northwestern University, Evanston, Illinois 60208
Simon W. de Leeuw Laboratory for Physical Chemistry, Nieuwe Achtergracht 127, 1018 WS Amsterdam, The Netherlands, and Department of Computer Science, Coates Hall, Louisiana State University, Baton Rouge, Louisiana 70803-4020 Received: May 14, 1993; In Final Form: July 1 , 1993e
Computer simulations of Stockmayer fluids were performed to generate dipole time correlation functions (TCF) at three temperatures and three dipole moments in both the microcanonical and canonical ensembles. The effect of Nos6 constant-temperature dynamics on time-dependent quantities is discussed, and empirical results are given to show that the choice of thermal inertia parameter influences the speed with which a system moves through its phase space. The time correlation functions from the simulations were analyzed in terms of current theories for dipolar systems. A functional form is proposed to cover both the longtime and short-time behavior of the time correlation functions of dipoles. The relationship between this functional form and the dielectric function of the Stockmayer system is also discussed.
N
I. Introduction This paper reports a series of simulations using Stockmayer particles, which have permanent point dipoles imbedded in spheres of dimensions specified by Lennard-Jones parameters. Stockmayer fluids are simple systems, but neither molecular dynamics simulationsnor solution theory has progressed to the point where their behavior is completely understood. This study focuses on the time correlation functions of Stockmayerfluids. Transport coefficientsand dielectricconstants, although not time-dependent themselves, are conveniently calculated using time correlation functions. These functions also appear in the analysis of a wide variety of transport processes.' For a one-component system of N dipoles each with dipole moment ~ [ ( tthe ) relevant functions are the single particle and collective dipole time correlation functions
where the collective dipole is defined as M(t) = CElr,(t).i W andpzaredefined as ( (M(0))2)and ( (C(~(O))~),respectively. The cross-correlation function BX(t)between two different particles is related to the @.,(t) and @,(t) functions Author to whom correspondence should be addressed. address: Department of Materials Engineering, Monash University, Wellington Rd., Clayton, Victoria 3 168, Australia. t Present address: E. Hala Laboratory of Thermodynamics, Institute of Chemical Process Fundamentals, 165 02 Praha 6-Suchdol, Czech Republic. *Abstract published in Advance ACS Abstracts, September 1, 1993. 7 Present
N
(4)
This study uses the dynamics formulated by Nos6 for the purpose of carrying out canonical simulations,2 which superseded the conservation of temperature by the frequent rescaling of ~elocities.~Nos6 offered proof that if the extended system produces configurations in the microcanonical ensemble, then the real system will give correct results for the canonicalensemble, at least as far as the calculation of static averages is concerned.* However, the proof that Nos6 dynamics will produce canonical ensemble averages of time-dependent quantities has never been given. Hoover has shown that with Nos6 dynamics in single oscillator systems the canonical ensemble may never be reached because the number of degrees of freedom is not enough to guarantee ergodi~ity.~ Hoover has offered an alternate formulation of constant-temperature dynamics employing a thermodynamic friction coefficient.M A generalized version of constanttemperature dynamics using an extended phase space which appears to provide ergodic behavior for even very simple systems has been developed.' Another approach which employs a chain of variables and may be easier to use than the Kusnezov method has also been published.* Recently, a alternative to the original Nos6 dynamics has been proposed by Litnie~ski.~ The simulation results of both Nos6 and NosbHoover dynamics will be the same for ergodic systems such as those under study here. The choice of whether to employ Nos6 or NosbHoover dynamics is a matter of personal preference. For Nos6 dynamics,
0022-3654/93/2097-10478$04.00/0 0 1993 American Chemical Society
Time Correlation Functions of the Stockmayer Fluid
The Journal of Physical Chemistry, Vol. 97, No. 40, 1993 10479
TABLE I: Reduced Units Used in the Simulations' parameter unit of Mass: particle Mass unit of distance: Lennard-Jonesu unit of energy: Lennard-Jonese time time step radial frequency, wavenumber, and frequency
symbol
m nu
kbT = C
U
t
At w, cm-I,
Hz
particle no. density and concentration moment of inertia temperature dipole moment
cc
thermal inertia: energy-time2
Q
r e d u d units m = 1.0 a* = 1.0 e* = 1.0
real units 40.0amu 3.405A 120 K 1.66 X 10-21 J/particle 0.997 kJ mol-' 2.15596 X 1&l2 s 4.31192 X s 2.46 cm-I 0.0738 X 1OI2Hz 256 49.3 A3/particle 33.7mol L-1 19.3 X l e g cm2 162 K 240 K 360 K 0.4044D 0.8087D 1.2131 D 5.58 X 104 K (ps)2 7.72X l(r3 J s2 464 kJ (ps)2 mol-'
= 1.0 A P = 0.002 w* = 1.0 t*
N = 256 p*
= 0.80
I+ = 0.025 kbp
1.35
k b p = 2.00 k b p = 3.00 p* = 0.5 cc* = 1.0 p* = 1.5 Q* = 100
Their values in real units are also shown for comparison. the best value to use for the Q parameter is unknown, whereas, in NosbHoover dynamics, the value of the relaxation time is unknown. With the NosbHoover formulation, the implementation of the dynamics is simpler, but the elegance of a holonomic system and Hamiltonian structure is lost.'
II. Simulation Methodology Classical molecular dynamics experiments were carried out on 256 or 500 Stockmayer particles using either typical constantenergy dynamics or Nos6 constant-temperature dynamics. The Hamiltonian of a Nos6 system may be expressed as
+
/--*-------$
p.
4
4.200
> 1
-0.400 0.0
-u
0.1
0.2
0.3
0.4
0.5
0.6
Time I Reduced Units
where Vral and Krul are the potential and kinetic energy of the real system, Tis temperature, k b is the Boltzmann constant, and f is the number of degrees of freedom of the real system. Q is the thermal inertia parameter, and s andp, are the Nos6 parameter and the corresponding momentum. Simulations were performed on a CRAY Y-MP f 832 and a STARDENT ST 2500. Long-range forces were treated with the Evald summation technique10Jl using the periodic image convention3and "tinfoil" boundary conditions with the external region having a dielectric ced of infinity. Production simulations cover a real time span of 430 ps. The absolute value of bi(t) was constrained to be constant during each simulation. Time correlation functions were calculated to 4.3 ps. Errors in static quantities were estimated using the blocking method.12 Units reduced with respect to 6-1 2 potential Lennard-Jones parameters of the fluid were used in the program and in derived properties. All simulations were confined to the one-phase region of the Stockmayerfluid phasediagramas documented in the1iterat~re.l~ Table I gives a summary of the parameters used to represent the Stockmayer fluid. The cross-correlation function BXwas calculated from 9,and 9cproducedduring each simulation. Curve fits of time correlation functions were obtained from the Macintosh KaleidaGraph program. The general curve fitting used unweighted least squares to tit data to a given function form. The particle independenceof the results was investigated using 500 particles in a test simulation. Static and time-dependent results were the same as simulations with 256 particles. Such systems are reported to behave with macroscopic dielectric properties for particle numbers of 256 and up.14 Tests using 108 particles carried out during the initial stages of this research
Figure 1. Single(*.I) and collective (ac)time correlationfunctions (TCF) for the microcanonical and canonical ensemblhs for the system with p* = 0.5 and Ts = 3.00. The size of the associated standard deviation error bars is shown in Figure 2.
showed a poor signal-to-noiseratio in the time-dependent results. Results apply to simulationsperformed in the canonicalensemble, except where otherwise stated.
III. Results and Discussion Nine simulationsin the canonicalensemble using three different dipole moments and three different temperatures were carried out on a Stockmayerfluid system of 256 particles. Thesimulations were repeated in the microcanonical ensemble. Figure 1 compares the results from the two ensembles for one set of conditions. Under all conditions the results from the two ensembles could not be distinguished within experimental error. Figure 2 shows the standard deviationerror bars, calculated from the results of seven production runs, for the same conditions in the canonicalensemble. The time correlation functions from the canonical ensemble experiments were analyzed to characterize dipole decay in Stockmayer fluids. The dipole moments used in the simulations were limited to low to moderate values, since Kusalik has extensively investigatedsimulations of dipolar particles with high dipole moments.15 The analysis of the behavior of particles with low dipole moments requires comparatively more attention to short-time behavior. Another series of experiments was carried out to investigate the effect of Nos6 dynamics for the system with T* = 3.00 and cc* = 0.50 both away from and at equilibrium. In these experiments the Q parameter was varied between 1 and 5000 reduced units. Figure 3 shows the effect of varying Q on the
10480 The Journal of Physical Chemistry, Vol. 97, No. 40, 1993 1 .Ooo
0.800 0.600 T.C.F.
0.400 0.200
-0 200
O.Oo0
1 c
1;;
-0400
' ! ; ' ; t i - ; !
000 005
010
! ! ! ! ! !
! ! ! " ! ! " ! ! '
0 1 5 020 0.25 030 Time I Rdueed Units
040
035
Figure 2. Standard deviation error bars for the system in the canonical ensemble with p* = 0.5 and Ts = 3.00. The error bars were calculated from the results of seven different simulations of the system, each of which was 430 ps long.
.........Q'
=
1000
7.0
I +--t
6.0
5.0 Energy 4.0 kJ mole-'
3.0
2.0 1.o
0.0 ! 0.0
!
'
~
2.0
! !
:
:
!
4.0 6.0 Time I ps
I
so
,
,
,I 10 0
Figure 3. Effect of three differentvalues of the thermal inertia parameter Q for the system with g* = 0.5 and Ts = 3.00. The total energy of the
real system is plotted versus time for systems initially in a lattice configuration.
energy in the real part of the system for three of the simulations beginning from a nonrandom configuration in which the particles are initially placed on a cubic lattice. The behavior of the systems was observed moving toward equilibrium: the instantaneousstatic quantities undergo damped harmonic oscillations until finally showing random noise modulation with no significant drifts in the averages (slight downward drifts in the extended system energy are normal). A. Nose Dynamics. The effect of Nos€ dynamics is much more obvious during equilibration than during production. As can be seen from Figure 3, the random fluctuations in the total system energy as the system relaxes toward equilibrium have larger amplitude and smaller frequency with a larger value for Q. However, a thermodynamicformula relates the averagesquare energy fluctuation ((6E)Z)in a canonical system at equilibrium to the constant-volume heat capacity C, of that system, which is a characteristic of the real part of the system and should not depend on the choice of Q:16
At long time during the simulations with large Q (on the order of 103 reduced units), the fluctuations become larger and eventually approach the amplitude of those in the simulations with smaller Q. The large Q dramatically slows down the rate at which the system samples the available phase space. Nos6 dynamics do not approach the microcanonical ensemble for very large Q, contrary to what a casual inspection might suggest. Instead, the system evolves more slowly through the canonical phase space. This interpretation is supported by an examination
Payne et al. of the system Hamiltonian: the simple and proper way to convert the algorithm from the canonical to the microcanonicalensemble is by setting the s parameter to unity and its derivatives to zero at each time step. On the other hand, a small Q produces sharp spikes in the static properties during equilibration. The exchange of energy between the real and extended system becomes so rapid with very small Q that the real particles no longer exhibit Newtonian behavior. The selection of an appropriate value for Q, which formally has units of energy-time2,Z is a matter of trial and error. The main criterion for the choice of a value for Q is that a simulation in the canonical ensemble should progress both smoothly and rapidly through the available phase space. The value must be large enough so that the real system displays Newtonian dynamics over short times, but not so large that temperature relaxation does not occur over the time span of the simulation. In the simple Stockmayer fluid, the precise choice of Q is not critical. In most of these simulations discussed in this paper, Q was set to a value of 100 reduced units. During someequilibrationswith Nos€ dynamics,whether from the lattice or a different starting point, sudden exponential increases in the s parameter and the extended system energy have been observed. A suggested method for reachingequilibrium has been found: the s parameter and its derivatives are reset at the first steps of one to three short runs. The phenomenon is probably due to energy trapped in the 1 degree of freedom of the s parameter. NosC has given a similarexplanationfor a simulation which did not reach equilibrium.* The influence of Nos€ dynamics on the time-dependent quantities in an equilibriumsystem is a more subtle issue. Pollock and Alder performed simulationsin the microcanonical ensemble, generating smooth curves for aCand 9,at a variety of dipole moments." The simulations reported here were performed in the canonical ensemble and then repeated in the microcanonical ensemble. The time correlation functions obtained agree with Pollock and Alder's results and with each other within experimental error, as Figures 1 and 2 show. Furthermore, the simulations in the canonical ensemble with different values of Q reveal no pattern relating the dipole behavior to the size of Q. The N o d dynamicsmight be predicted to influencea constanttemperature simulation in the following case: the temperature, which is controlled by Nos€ dynamics,has its own time correlation function, and if its fluctuations are on the same time scale as a real time-dependent quantity of interest, the artificial nature of thedynamics may taint the results. Whether or not the numerical results of the two ensembles agree, it cannot be known for certain how well the NosC dynamics mimics real constant-temperature dynamics without an analytical proof. In these simulations, observation shows that the fluctuations of the Nos€ parameter are on a much faster time scale than the decay of aCand a,, and it will be tentatively assumed that theNosCresu1t.s arecomparable to real constant-temperature dynamics. The conclusion is that the two dipole functions are not sensitive to the simulation ensemble, which in itself is an indication that this system is sufficiently large to be considered macroscopic. B. Time Correlation Functions. The evolution of both and aCis due solely to the rotational reorientation of the point dipoles in the system. 9,necessarily is affected by cross-correlations between dipoles. Rothschild makes several points which will be important in this analysis.l* First, the dipoles will behave at very short times as free spinning rotors without the influence of intermolecular forces. A Taylor series expansion shows that at short time 9,may be approximated as ( QZ ) t2
9,=1--
2
(7)
where t is time and (@) is the mean-squared angular classical velocity. The relation between ( Q2) and the moment of inertia
Time Correlation Functions of the Stockmayer Fluid
The Journal of Physical Chemistry, Vol. 97, No. 40, I993
TABLE II: Theoretical Radical Frequencies (a*) Corresponding to the Initial Decay of Free Rotors at the Temperatures (kb2')and Moments of Inertia (I*) Used in the Simulations' kbT* I* n* 1.35 0.025 10.39 2.00 0.025 12.65 3.00 0.025 15.49 e The asterisks indicate the values are reported in reduced units.
TABLE III. Fits of the Time Correlation Functions. The parameters' o*,l/u*T*, and 1 / ~ *Calculated from the Fits of the asTime Correlation function to the Given Equations (a) For the Equationb @.,(t)= cos(ot) exp(-(a2
dipole moment(p*) temp(T+) 0.5
Z of the particles is19 (Q2)
1 .o
= 2k,T/I
Taking Tas the average temperature in a simulation leads to the theoretical values for Q shown in Table 11. A negativevalue for a timecorrelation function at intermediate times indicates that the dipoles have on average rotated more than 90' from their initial positions. Furthermore, at long times 0, will decay to zero.18 The long time decay is expected to be exponential. This behavior is often called Debye relaxation. In most cases, including the ones discussed here, exponential decay at long times is due to the Markovian nature of the process in which the development of the system after any time tA depends only on the state of the system at time t ~ . l 8 The relation between single and collective dipole relaxation has been treated using analytic approximations. An analysis for polar liquids was recently reported by Chandra and Bagchi using the generalized Smoluchowski equation;20-21they obtained (9)
where T.,, ~ , land , rC2arethe single dipole relaxation time, collective transverse polarization dipole relaxation time, and collective longitudinal dipole relaxation time, respectively. The authors predict T., and ~ , to l have approximately the same magnitude, while 7 ~ is2 smaller, corresponding to faster relaxation. Expressions in terms of the pair distribution function are given for the constants A, A I ,and Az. These conclusions agree qualitatively with other literature results.Z2 Of course, these are not the first papers to analyze the separate single-particle and collective contributions to correlation functions. That discussion has been going on for d e c a d e ~ ~ 3 . 2 ~ A great amount of work remains to be done to bridge the gap between analytical theories and the time correlation functions obtained from simulations or spectroscopy experiments.18 In particular, any treatment based on Langevin-type descriptions with &function memory will not deal correctly with the (inertial) short-time dynamics; such treatments generally lead to exponential short-time decay. Time correlation functions are often discussed in terms of exponentials such as the following E ( t ) : n
Aa a= 1
Ea(t) = exp(-d,t)
(1 1b)
However, the following function form J ( t ) provides a more appropriate form for 9,:
10481
1.5
1.35 2.00 3.00 1.35 2.00 3.00 1.35 2.00 3.00
+ (t/7)2)'/2)
u*
l/u*r*
1/7*
8.76 10.67 12.87 7.57 9.97 12.79 2.0 7.49 11.13
3.47 3.33 4.20 6.49 4.70 5.03 9.23 10.0 8.29
9.01 12.7 14.7 8.76 13.2 15.6 11.9 11.3 14.9
W*
eo
10.39 1.76 12.49 1.48 15.08 1.29 10.68 5.55 12.70 3.71 15.56 2.58 10.67 20.1 13.03 9.70 15.74 5.60
(b) For the EquationC
a,(t)= cos(wt) exp(-(u4 + (r/7)4)1/4) dipole moment ( p * )
temp(TS)
u*
l/U*T*
I/r*
W*
P~
0.5
1 .o 1.5
1.35 2.00 3.00 1.35 2.00 3.00 1.35 2.00 3.00
6.06 7.55 6.6 4.61 6.37 8.13 1.7 4.27 6.47
1.75 2.27 2.38 2.24 2.75 2.99 2.22 2.73 3.69
12.1 16.8 20.4 12.7 16.4 18.9 10.5 14.4 20.9
7.60 1.23 9.76 1.14 9.60 1.02 7.05 1.83 9.26 1.61 11.08 1.42 5.13 3.43 7.58 2.31 10.91 1.83
e The asterisks indicate that the values are reported in reduced units. W* = [ ( o * ) ~ + ( 1/:*(7*)2]1/2is calculatedfor comparison with ll* from Table11. The static dielectricconstant g calculated foreach simulation is also shown. The correlation factor g, calculated from @,(O) is also reported.
J(t) is in some respects a better functional form for a time correlation function than in exponential such as E ( t ) . First, J ( t ) is an even function, which is required for all real time autocorrelation functions.25 The cosine term is necessary to describe negative regions of the time correlation function. The u term provides for the correct linear free rotor decay at short times required by eq 7; at times exceeding 47, J ( t ) essentially becomes an exponentially decaying cosine. Thus, J ( t ) is physically reasonable, as it provides for the two limiting cases of free rotor decay of rotating dipoles and dipoles with strong interparticle forces. The argument can be made that the short-time and longtime behavior are separable, and only the latter is important. However, this process requires a somewhat arbitrary division of the data and can lead to the neglect of interesting short-time effects. The fits obtained for 0#and 9,are reported in Table I11 and shown in Figure 4. Data were used in the time region where the noise was small. The values reported are for unweighted data. Fits calculated with weighted data showed little difference in the results. The single-particle time correlation function 9,should be fitted to eq 12 with n = 1. In the three cases of the smallest dipole moment, it should be noted that the fits consistently rise more rapidly to the second maximum than the real function,
Payne et al.
10482 The Journal of Physical Chemistry, Vol. 97, No. 40, 1993
@~
[;=0.51
0.0010
1\ )
0.4
6 $
0.0012
+\
0.0008
i
T'=2'00
T'=3.00
0.2
-2.00
O.OW6
\
O.ooo4 @x
0.0002
O.oo00 4,0002
0.ooo4 4.oo06
0.2
0.0
0.4 0.8 lime I Reduced Units
1.o
0.8
0.0
0.1
0.2 0.3 0.4 lime I Reduced Units
0.5
0.6
T'-2.00 T'=3.00
0.2
-
# \ &4oo . o
f
0.002
0.0 1
0.2 -4.4
0.000
, : : I ;
.
- ....
-.
--
-----i
; I ; : : ; : : ' ; : : :
Figure 4. Time correlationfunctions(a, left side) ip, and (b, right side) ipx shown for all v&cs of the dipole moment and temperature. @, is normalized by definition, while the magnitude of *x indicates the relative magnitude of intermolecularinteractions. The functions obtained from simulation are shown as full lines, while the calculated fits reported in Table 111 for each function are shown as dashed lines.
which may be attributed to inertial effects having relatively more importance with the comparatively small dipole moment, since interparticle interactions scale with ~ 2 . From equating the short-time behavior of eqs 7 and 12b, one obtains
2k,T - W 2 + 71 --
Z
(13)
UT
Equation 13 places a constraint on the parameters w, u, and 7 for as.The theoretical inertial frequency S2 = (2kbT/1)1/2 at each temperature in the case of free rotors is shown in Table 11. The value 51' = (02 (l/u~2))1/2calculated from the results of the fits is shown in Table 111. The close agreement between S2 and $2' indicates that the quality of the asfits is good. Equation 13 does not hold for ax,since the cross-correlationeffect will not lead to free rotor decay even at short times; however, values for $2' for 51, fits are also reported for comparison. It should be remarked that eq 13 may be used to constrain one of the three parameters in eq 12. Fits for as were also calculated by constraining w and determining u and 7 ; the quality of the fits was not substantially changed. The interparticle interactions scale as p2p/kbT, which will be reflected in the dipole moment and temperature dependence of the time correlation functions. The parameters w , UT)-^, and r l
+
from the asfits in Table I11 indicate the speed of rotation of the dipoles, the importance of short-time decay, and the importance of long-time exponential decay, respectively. Clearly, w should increasewith increasingtemperature and decreasewith increasing dipole moment. The parameter w has an approximately linear dependenceon F / 2 as , suggestedby eq 13. The parameter UT)-^ increases with increasing dipole moment but shows no strong correlation with temperature. In the free rotor regime (no collisions) that UT)-^ represents, decay is governed by the rotational kinetic energy of the particles. The parameter T-1 shows the tendency to increase with increasing temperature but shows no strong correlation with dipole moment. The longtime behavior of asis exponential decay caused by collisions; thus, it is reasonable that the size of r1be temperature-dependent. It is reasonable to expect that axshould fit to a equation similar to the one used for a,, where T is smaller than that for as, correspondingto collective longitudinal dipole relaxation as in eq 10. This correlation function must be even.25 Because of the finite time required for propagation of interactions through space, the second derivativemust also be zero accordingto the constraint
(%%) =o at at when i and j denote different particles. Therefore, the lowest
nonvanishing time derivative of 9,(t=O) is in fact the fourth derivative. This might imply that a generalization of eq 12 in the form
J(t) =
a-1
n = 1,2,3, ...
g, = M z / N p 2 (16) for a system with external dielectric constant eat approaching infinity. The correlation factor,g is dependent on eat, and a general formula for periodic boundary conditionshas been given by de Leeuw, Perram, and Smith.27 Under the condition that cut is equal to EO of the system being simulated, g, is replaced by the Kirkwood correlation factor gk. The relationship between 9, and g, is seen by comparison of eq 16 with eq 4:
+1
(17)
The amount by which g, exceeds unity, which represents collective decay, has the largest magnitude for the largest dipole moment and decreases as temperature rises. Equation 17provides an alternate method for the calculation of g,, which should be independent of the number of particles in the simulation. The values of g, calculated from %(O) are reported in Table 111. Higher values of g, for somewhat less mobile systems (as would be expected) have been reported in the 1iterat~re.l~ The results obtained in this paper definitely support the theory of Chandra and Bagchi. The long-time behaviors of 9,and 9, are single exponential, which corresponds to single exponential decay of and double exponential decay of aCpredicted by eqs 9 and 10. The parameters 7 of asand T of axcorrespond to ~~1 and ~ ~ respectively, 2 , in eqs 9 and 10. These two relaxation times are different in each case except for p+ = 1.5, Ts = 1.35, where they are within experimental error. The reported values for 7' of axare larger than those for r1 of Os,correspondingto a shorter relaxation time for the collective longitudinal relaxation, as prdicted by Chandra and Bagchi. The values reported here for the two relaxation times are of the same order of magnitude, so the large difference in size predicted by Chandra and Bagchi for ~~1 and ~~2 is not observed. C. Dielectric Function and Dielectric Relaxation T i e . The behavior of the frequency-dependentdielectric function may be obtained from 9c. Caillol, Levesque, and Weis28.29 have derived the following equations for nonpolarizable ionic solutions where M, is the solvent polarization and J is the ionic current E(W)
= 1 + =((M:) 3v
+ iw(ww),+ (WJ),)
--Ef
--
-
p 0 . 5 T*-1.35 -+- p*=l.o r-1.35 p*-lsr-1.35 ~ ' - 0 . 5 T'-2.00 -8- p'=l.O T'-2.00 --C p ' 4 . 5 T'92.00 P*-0.5 T*-3.00 p*=i.o r-3.00 p*-i.5 r-3.00
(15a)
might be a better description for ax. Both eq 12 and eq 15 were used to fit the axtime correlation functions. The fits produced from eq 15 provided better fits than those from eq 12 and are reported in Table I11 and Figure 4. The agreement between functions and their fits is not quite as good as those for the 9,functions. In particular, BXfor p+ = 0.5, Ts = 3.00 has a poor fit. However, in general, the axfunctions may be fit to eq 15 with n = 1. The parameters w , (UT)-', and 7 1 have the same approximate dependencies to temperature and dipole moment as discussed above for as,with the exception of the relationship between UT)-^ and dipole moment being less strong. 9,is directly related to the correlation factor g,, defined asz6
g, = ( N - 1)@,(0)
-
The Journal of Physical Chemistry, Vol. 97, No. 40, 1993 10483
Time Correlation Functions of the Stockmayer Fluid
(l8a)
4.2 4.2
0.2
0.0
0.6
0.4 c'
I
0.8
1.0
1.2
Le
Figure 5. Cole-Colediagramsplotting imaginary versus real parts of the dielectric function for all values of the dipole moment and temperature. A few points on each function are marked with symbols to distinguish the lines. Some deviation from smooth curvei is due to the inclusion of time correlation function points with relatively large error bars during the calculation of the dielectric function.
(AB),
I Jomdr e'"'(A(t)
B(0))
( 18b)
For nonpolarizable dipolar fluids the equations taken a simplar from14826
co=l+ ~ ( w= ) 1
+ (eo-
4PdM2> 3NkbT
1)(1
+ iw(@.,>,)
( 19b)
The static dielectric constant EO calculated during each simulation is shown in Table IIIa. The dielectric function of a system may be thought of as describing the interaction of a system with a harmonic electric field E&". The most convenient formulation calculates the interaction of the dipoles with the field that would be present in the absence of the particles, known as the vacuum field.22 The imaginary part of the dielectric function is proportional to the absorption of energy from the field. In general, the real part has a miximum at zero frequency, while the imaginary part has a maximum at some frequency 1/7d, where 7 d is the dielectric relaxation time.' Td is equal to ~~l in Chandra and Bagchi's theury.20J' The Cole-Cole plots for the dielectric functons from these simulations are shown in Figure 5. E' and e" refer to the real and imaginary parts of the dielectric function, respectively. Time correlation functions were smoothed by the Hamming window3O before the dielectric function was calculated. The dielectric function provides a link between theoretical systems and experimentally-measurable properties. A dielectric relaxation time Td of the system may be calculated from the dielectric function obtained from either a simulation or a laboratory experiment by using a suitablemodel such as the Debye, Cole-Cole, ColeDavidson, or Havriliak-Negami functions.31-34 The Debye function proved inadequate for modeling the dielectric functions discussed here. The other functions mentioned also failed to model the functions with suitable accuracy. Therefore, a dielectric relaxation time for each simulation was estimated from the peak frequency of the imaginary part of the dielectric function. These values are shown in Table IV. Also shown are the values for 7 d predicted from T , according to the Glarum-Powles
10484 The Journal of Physical Chemistry, Vol. 97, No. 40, 1993
Payne et al.
TABLE IF? Comparison of the Single Dipole Relaxation Time and the Dielectric Relaxation Time for Each System. dipole moment ( p * ) temp ( T ) co T,* Td* Glarum-Powles 7d' 0.5
1.35 2.00 3.00 1.35 2.00 3 .00 1.35 2.00 3.00
1.o
1.5
1.76 1.48 1.29 5.55 3.71 2.58 20.1 9.70 5.60
0.111 0.079 0.068 0.1142 0.076 0.0641 0.084 0.0882 0.0670
0.111 0.078 0.074 0.125 0.098 0.072 0.263 0.152 0.086
0.130 0.089 0.074 0.157 0.101 0.08 1 0.123 0.126 0.092
a i , * refers to the relaxation time for a, from Table IIIa. For comparison, the Glarum-Powlcs prediction for Td (eq 20) as calculated from T , is also shown. The asterisks indicate the values are reported in reduced units.
1.2
The calculated dielectric relaxation times show a strong dependence on both temperature and dipole moment. The predictions from the Glarum-Powles relation comparebest with the estimated values of 7 d for systems with low g, values. Chandra and Bagchi refer to a different relation as the GlarumPowles expression,21
$
from W =exp(4t)
0.0 -0.4 -0.2
0.0
0.2
0.4 E'/
With the use of equations relating g, and gk,27 the dielectric relaxation time 7 d predicted by eq 21 was found to be 0.288 for the p* = 1.5, T+ = 1.35 system, which compares better with the estimated value of 0.263 than the value from eq 20. This system has the highest reported g, factor. Equation 21 overestimates the calculated relaxation time for all the other systems. Cole-Cole plots of analytically-generateddata for the functions E ( t ) and J(r) from eqs 11 and 12 demonstrate that the latter but not the former shows a feature which is also seen in the simulation Cole-Cole plots. This feature, referred to as a "foot" in the literature,3*is a change in derivative sign at high frequency. The foot is also seen in the Cole-Cole plots of other Stockmayerparticle simulations14J7and can clearlybe attributed to inertial effects.22J8 Another featureobserved in the present work is the bowing out of the Cole-Cole plots at small frequency so that the following is true:
Bowing out is an effect of the rotation of the point dipoles, which can be shown by the analysis of a dynamical cosine term such as the one which represents this rotation in eq 12. For the simple single exponentialform for the collective time correlationfunction
@.,W= exp(--Yt)
(23a)
an analytical solutionfor the dielectric function using eq 19 shows that bowing out at low frequency cannot occur. (lime') - e'(w=O) = (eo - 1) -< O --WL
w-0
y2
+ w2
(23b)
Now consider an equation having a dynamical cosine term
@.,W= W 8 t ) exp(--Yt)
(244
0.6 0.8 1.0
1.2
Eo
F i p e 6. (a, top) Analytically-generatedtime correlation functions for three functional forms and (b, bottom) the associated Colbcole plots. Units are comparable to those in Figures 4 and 5. Note that only the complete form of eq 12 correctly fits the high-frequency foot and the low-frequency bowing-out seen in the Colbcole diagram from the simulations (see Figure 5).
Bowing out will occur at low frequency if @ is large: (lime') - c'(w=O) = A
(lime') - e'(w=O) A
>0
if and only if f12
> (7' + w2) (244
This analysis shows that bowing out will be most pronounced with small dipole moments (when /3 is large), which is supported by Figure 5. The rapid rotation of molecules with small dipole moments significantly affects short-time correlation function behavior and should not be neglected. Figure 6showsanalyticallygenerated time correlation functions and the associated ColeCole plots for exponential decay, decay according to eq 12, and decay accordingto eq 12 without the inclusion of the cosine term. Only the complete form of eq 12 provides for both the low-frequency bowing out and the high-frequency foot observed in the Cole-Cole plot of the simulation dielectric function. Gordon suggested that dipole correlation functions with negative regions are found only in gaseous systems, where reorientation is freer than in a liquid or solution.25 The Stockmayer fluid is nominally a liquid in the simulationsreported here, but apparently it retains this aspect of gaseous behavior, perhaps because of the spherical shapes of the particles which offer no steric resistance to rotation. IV. Conclusions
This paper describes a series of simulations from which time correlation functions of the Stockmayer fluid were obtained for
Time Correlation Functions of the Stockmayer Fluid both the microcanonical and canonical ensembles. No difference has been observed between the two ensembles for the equilibrium dynamical properties of the systems under consideration, but observations during equilibration suggest that the choice of the thermal inertia parameter Q greatly affectsthe speed of movement through phase space. Single particle (aII)and collective dipole (a,) time correlation functions were generated at different dipole moments and temperatures for systems in the canonical ensemble. A crosscorrelation function @,wasalso calculatedand related to collective longitudinaldipole relaxation in the theory of Chandra and Bagchi. The two time correlation functions 9, and a, were fit to a functional form designed to separate the effects of short-time free rotor decay and long-time exponential decay. This analysis supports the suggestionof Chandra and Bagchi that singleparticle decay and collective decay at long times in dipolar fluids may be described in terms of one and two exponentials, respectively. Negative regions in dipole time correlation functions lead to the inclusion of a dynamical cosine term in the functional form proposed for them. The negative regions suggest that although the Stockmayer fluid is in the liquid phase in the simulations reported here, it has some gaslike properties. The cosine term also accounts for the low-frequencybehavior of dielectric functions of the Stockmayer fluid.
Acknowledgment. We are grateful to S. D. Druger, M. Lonergan, J. W. Perram, J. A. Pople, N. Snider, and especially D. F. Shriver for helpful remarks. We are also grateful to the reviewers for their comments. Computer time on the CRAY Y-MP/832 was supplied by the Pittsburgh Supercomputing Center (DMR910014P). Thisworkwassupported bytheNSF/ MRC through the Northwestern MRL (Grant DMR 8821571) and by the ARO (Grant DAAL-03-90-G-0044). V.A.P. and M.F. gratefully acknowledge an NSF Graduate Fellowship and a Fulbright Postdoctoral Fellowship, respectively.
The Journal of Physical Chemistry, Vol. 97, No. 40, 1993 10485 (5) Hoover, W. G.; Ladd, A. J. C.; Moran, B.Phys. Rev. Len. 1982.48, 1818. (6) Posch, H. A.; Hoover, W. G.; Vesely, F. J. Phys. Reu. A 1986, 33, 4253. (7) Kusnezov, D.; Bulgac, A.; Bauer, W. Ann. Phys. 1990, 204, 155. (8) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2635. (9) Litniewski, M. J. Phys. Chem. 1993, 97, 3842. (10) de Leeuw, S.W.; Perram, J. W.; Smith, E. R. Proc. R . Soc. London A 1980. 373. 27. (1 1)’ de Leeuw, S.W.; Perram, J. W.; Smith, E. R. Proc. R . Soc. London A 1980, 373, 57. (12) Flyvbjerg, H.; Petersen, H. G. J . Chem. Phys. 1989, 91, 461. (13) Smit, B.;Williams, C. P.; Hendrih, E.M.; de Leeuw, S.W. Mol. Phys. 1989, 68, 765. (14) Neumann, M.; Steinhauser, 0.;Pawley, G. S.Mol. Phys. 1984,52, 97. (15) Kusalik, P. G. J. Chem. Phys. 1990, 93, 3520. (16) Chandler, D. Introduction to ModernStatistical Mechanics; Oxford University Press: New York, 1987. (17) Pollock, E. L.; Alder, B. J. Phys. Reu. Lett. 1981, 46, 950. (18) Rothschild, W. G. Dynamics of Molecular Liquids; Wiley: New York, 1984. (19) Halliday, D.; Resnick, R. Physics, 3rd ed.; Wiley: New York, 1978; Chapter 12. (20) Bagchi, B.; Chandra, A. Ado. Chem. Phys. 1991,80, 1. (21) Chandra, A.; Bagchi, B. J . Phys. Chem. 1990,94, 3152. (22) Madden, P.; Kivelson, D. Adu. Chem. Phys. 1984,56,467. (23) Craddock, H. C.; Jackson, D. A.; Powles, J. G. Mol. Phys. 1968,14, 373. (24) Keyes, T.; Kivelson, D. J. Chem. Phys. 1972,56, 1057. (25) Gordon, R. G. Adu. Ma@. Reson. 1968, 3, 1. (26) de Leeuw, S.W.; Smit, B.;Williams, C. P. J . Chem. Phys. 1990,93, 2704. (27) de Leeuw, S.W.; Perram, J. W.; Smith, E. R. Annu. Rev. Phys. Chem. 1986,37,245. (28) Caillol, J. M.; Levesque, D.; Weis, J. J. J. Chem. Phys. 1989, 91, 5544. (29) Caillol, J. M.; Levesque, D.; Weis, J. J. J . Chem. Phys. 1989, 91, 5555. (30) Emst, R. R.; Bodenhausen, G.; Wokaun, A. Principles ofNucleor
IL..~1976. (2) Nose, S. Mol. Phys. 1984, 52, 255. (3) Allen, M. P.; Tildesley,D. J. ComputerSimulationof~iquids; Oxford . . University Press: New York, 1987. (4) Hoover, W. G. Molecular Dynamics; Springer-Verlag: New York,
Magnetic Resonance in One and Two Dimensions; Oxford University Press: New York, 1987; Chapter 4. (31) Havriliak, S.;Negami, S.J . Polym. Sci.: Parr C 1966, 14, 99. (32) Davidson, D. W.; Cole, R. H. J . Chem. Phys. 1951, 19, 1484. (33) Cole, K. S.;Cole, R. H. J . Chem. Phys. 1941, 9, 341. (34) Jonscher, A. K. Dielectric Relaxation in Solids; Chelsea Dielectrics Press: London, 1983, Chapter 3. (35) Deutch, J. M. Faraday Symp. Chem. Soc. 1977,11,26. (36) Cole, R. H. J. Chem. Phys. 1965, 42,637. (37) Powles, J. G. J. Chem. Phys. 1953, 21,633. (38) Lobo, R.; Robinson, J. E.; Rodriguez, S.J . Chem. Phys. 1973, 59,
1986.
5992.
References and Notes (1) McQuarrie, D.A. StatisticalMechanics; Harper & Row: Evanston, ~~
