4916
J. Phys. Chem. 1987, 91, 4916-4922
The fact that hard spherocylinders can have both a nematic and a smectic phase may explain the problems that VieillardBaron” experienced in his search for the isotropic-nematic transition in a system of 616 hard spherocylinders with L I D = 2. Vieillard-Baron found in his Monte Carlo simulations that extremely long preparatory runs were required to generate an equilibrated (isotropic) fluid configuration, starting from the crystalline solid. In retrospect, it seems likely that VieillardBaron’s simulations were plagued by critical slowing-down effects due to the vicinity of both the isotropic-nematic and the nematic-smectic phase transitions. Whether the system studied by Vieillard-Baron does in fact have thermodynamically stable nematic and smectic phases is still an open question. Simulations on spherocylinders with LID ratios other than 5 are currently under way. To conclude, let us return to Onsager’s model system consisting of long, thin spherocylinders. As mentioned earlier, this model plays a special role in the statistical mechanical theory of liq(24) Dowell,
F. Phys. Rec. A 1983, 28, 3526.
uid-crystal formation because it is the only exactly solvable model that predicts a transition to the nematic phase. However, in the strongly aligned nematic phase the “second virial” theory is no longer valid (see footnote 6). Actually, a breakdown of the Onsager approximation at high densities is to be expected, because at sufficiently high packing fraction thin spherocylinders should form a crystal, rather than a liquid crystal. However, the present simulations on spherocylinders with a finite L I D ratio (both aligned and nonaligned) suggest that, before they freeze, Onsager’s spherocylinders may well form a smectic liquid crystal. From the work of Rahman, I learned that the “complex” behavior of many-body systems can often be explained in terms of very simple intermolecular interactions. It is gratifying to see that this approach is fruitful, even for systems as complex as liquid crystals.
Acknowledgment. I gratefully acknowledge many stimulating discussions with Bela Mulder, Henk Lekkerkerker, and Alain Stroobants. Mike Allen contributed most to the development of a vectorizable molecular dynamics code for nonspherical hard-core molecules which was used in this work.
Frequency-Dependent Specific Heat in a Simulation of the Glass Transition Gary S. Crest* Corporate Research Science Laboratory, Exxon Research and Engineering Company, Annandale, New Jersey 08801
and Sidney R. Nagel* The James Franck Institute and Department of Physics, The University of Chicago, Chicago, Illinois 60637 (Received: December 4, 1986; In Final Form: February 19, 1987)
Recent experiments have shown that there is a frequency dependence to the specific heat of liquid glycerol as it is cooled through the glass transition. We show that in a simulation of a simple Lennard-Jones fluid the same phenomenon can be observed. There is a large difference between the behavior of the constant pressure and constant volume specific heats. The former shows a pronounced drop as the temperature is lowered while the latter shows very little signature of the transition at all
Introduction The glass transition occurs when a liquid is cooled sufficiently far below its equilibrium freezing temperature so as to become an amorphous solid. To reach this glass state the liquid must be supercooled rapidly in order to prevent the nucleation of the crystalline phase. The nature of the glass transition is still open to debate. For example, there is ongoing controversy over whether all classical liquids can form a glass if cooled rapidly enough or whether only liquids which have certain properties can be glass formers. It is also unsettled whether the glass transition is a true phase transition or whether it is simply a kinetic freezing of the liquid. Experimental attempts at studying this phenomenon have been hampered by one of the essential characteristics of the glass transition itself. The liquid becomes so viscous as the transition temperature, Tgris approached that it becomes progressively more difficult to keep the sample in equilibrium. Thus the time scale determined by the experiment depends either on the cooling rate or on the frequency of the probe being used and static susceptibility cannot be measured near TB. On the other hand, theoretical investigations of the glass transition have not been able to agree on the microscopic mechanism which is responsible for the phenomenon. Should the same mechanism be responsible for the glass formation in polymers as well as in metals or does one need 0022-365418712091-4916$01.50/0
different theories for different types of systems? It may seem paradoxical that a computer simulation of the glass transition could provide much help at sorting out these various issues. After all, glass transitions are inextricably bound up with long relaxation times whereas computer simulations can only be run for relatively short times. The slowest cooling rate used in a computer simulation is much faster than can be achieved in even the most rapid experimental quench. The simulation is thus forced to study glass formation very far away in temperature from what would occur in an ideal, infinitely slow, cooling rate experiment. Another problem with simulations is that they usually are on systems with very simple interparticle interactions. It is much more difficult to simulate a liquid of complicated molecules interacting via long-range forces. Despite these drawbacks computer simulations have provided a number of insights about the glass state and about the transition itself. Rahman and co-workers’ studied the properties of a frozen Lennard-Jones fluid and showed that it behaved much like a solid. They also investigated the nucleation rates for various potentials and showed that in some cases the sample could remain amorphous during very long simulation^.^^^ It was also found4 that the ( I ) Rahman, A.; Mandell, M. J.; McTague, J. P. J . Chem. Phys. 1976, 64, 1564.
0 1987 American Chemical Society
Simulation of the Glass Transition structure of the Lennard-Jones model glass was very similar to that found in diffraction studies of real systems. More recently Fox and Andersenss6 showed that the Lennard-Jones liquid displays many of the same phenomena that are found in glass transitions observed experimentally in the laboratory. Among other things a drop in the slope of the enthalpy vs. t e m p e r a t ~ r e ~and -~ a dependence on cooling rate of the low temperature density was found.s,6 This latter effect was the first time in a simulation that the cooling rate was found to effect the properties of the low temperature glass. These simulations show that the Lennard-Jones fluid does exhibit all of the signs of a glass transition. We believe that this implies that all classical liquids, if cooled sufficiently quickly, will form a glass. After all, one would expect that argon (which is well simulated by a Lennard-Jones interaction) would be one of the hardest systems to keep from crystallizing. Since the glass transition is seen in a sample with only short-range central forces, it becomes clear that molecular complexity is not necessary for glass formation and therefore it becomes valuable to study glass formation on the simplest possible system. In this paper we will show that molecular dynamics simulations show other signs of the glass transition. A recent experiment7,* showed that there is a frequency-dependent specific heat at the glass transition of glycerol. We show here that the same frequency dependence can be found in a simulation. Specific heat measurements play a fundamental role in the study of glass transitions because it is the basis of the Kauzmann para do^.^ Kauzmann showed that if a drop in the specific heat did not occur at a temperature not far below where the system fell out of equilibrium then the entropy of the liquid would become much less than that of the corresponding crystal. Thus the study of the specific heat is perhaps the most valuable probe of the glass transition: the specific heat is constrained by the Kauzmann paradox and it is the only probe that should couple equally to all the modes. In the experiment the frequency-dependent specific heat was obtained by measuring the ac temperature oscillations in response to an oscillation of the heat current at the same frequency. It may seem contradictory to have a frequency-dependent specific heat since it is usually supposed to be a static thermodynamic quantity. However, Oxtobylo has shown that the same results are obtained by using the naive notion of a frequency-dependent specific heat as are found from his more careful analysis in terms of a static value of cp and the explicit inclusion of slowly relaxing modes. As we show below the frequency-dependent specific heat can be obtained from the temperature (or kinetic energy) autocorrelation function in a simulation. As in the expeiiment the specific heat is a complex quantity with both a real and an imaginery part: c = c ' + ic". The Model and Method
In our simulations we worked with a 500-particle system with periodic boundary conditions. Since our simulation runs were very long we needed to decrease the chance of nucleation. Consequently, we made our sample out of a binary mixture of particles with two different radii but with the same mass. The smaller particles (which made up 20% of the sample) had a radius uss= O . 8 q where us and uIIare the radii of the small and large particles, respectively. The depth of the potential was the same for all pairs of particles. In order to ensure energy conservation to better than a few parts in lo6 over the course of a long run two conditions (2) Hsu,C. S.; Rahman, A. J . Chem. Phys. 1979, 70, 5234. (3) Hsu, C. S.; Rahman, A. J. Chem. Phys. 1979, 71, 4974. (4) Wendt. H. R.; Abraham, F. F. Phys. Rev. Lett. 1978, 41, 1244. (5) Fox, J. R.; Andersen, H. C. Ann. N.Y. Acad. Sci. 1981, 371, 123. (6) Fox, J. R.; Andersen, H. C. J . Phys. Chem. 1984, 88, 4019. For a review of many simulations of the glass transition see Angell, C. A,; Clarke, J. H. R.; Woodcock, L. V . Adu. Chem. Phys. 1981, 48, 397. Nose, S.; Yonezawa, F. J . Chem. Phys. 1986,84, 1803, have seen similar behavior in their simulations of the glass transition. (7) Birge, N. 0.;Nagel, S. R. Phys. Reu. Lert. 1985, 54, 2674. (8) Birge, N. 0. Phys. Reu. B 1986, 34, 1631. (9) Kauzmann, W. Chem. Rev. 1948, 43, 219. (10) Oxtoby, D. J . Chem. Phys. 1986, 85, 1549.
The Journal of Physical Chemistry, Voi. 91, No. 19, 1987 4917 TABLE I: Temperature T/t, Enthalpy H/t, and Average Density pu3 for the Initial Cooling Run' 0.81 0.72 0.60 0.50 0.40
-2.482 -3.019 -3.578 -4.101 -4.590
0.8370 0.8833 0.9427 0.9895 1.0326
0.30 0.20 0.10 0.02
-5.023 -5.379 -5.708 -5.957
1.0663 1.0903 1.1120 1.1269
Each system was run at constant pressure for 607 before the temperature was lowered by 0 . 1 ~The pressure was Pu3/c = 0.5 and the coupling to the heat bath was I'lm = 0.5.
TABLE 11: Temperature T / c , Enthalpy H/t, Average Density pu3, Time Step A t / 7 , and Number of Time Steps M for Each System for the Constant Pressure Microcanonical Runs (r = O ) a All7 M PO3 T I6 HI 1.02 -1.371 0.722 0.0030 64 000 0.80 -2.5204 0.841 0.0025 160 000 -3.6442 0.940 0.0025 0.59 96 000 0.52 0.39 0.32 0.17 0.02
-4.0141 -4.6441 -5.0051 -5.4893 -5.9694
0.982 1.037 1.067 1.098 1.126
0.0025 0.0025 0.0025 0.0030 0.0050
256000 320 000 288 000 256 000 96 000
"Pressure Pu'lt = 0.5
were necessary. First, the interatomic potential between two particles a distance r apart was chosen to be a Lennard-Jones potential modified to have a smooth first derivative. U ( r ) = 4 ~ S ( r ) [ ( u ~ ,-/ r(uj,/r)'j ) ~ ~ - ( u j j / r C ) l+ 2 (u2,/rJ6] with S ( r ) = 1.0 for r
< rl
= 1.0 - ( r - r1)2(3rc- rl - 2 r ) / ( r c- r l ) 3 for r, < r
= 0 for r
< rc
(1)
> rc
and rc = 2.3ui, and rl = 1 . 9 ~ , ~This . smoothing function, S ( r ) , was used by Stillinger and Rahman" in their study of water. The second necessary condition was that the time step, At, had to be 2 the kept small, typically At/T = 0.0025 where T = u , , ( m / c ) ' / is fundamental unit of time in our simulation. This is a factor of two smaller than is used in most other simulations of LennardJones liquids. The simulation was first run at high temperatures to equilibrate the liquid. We cooled the liquid slowly starting at T / c = 0.8 and ran for 607 before dropping the temperature by 0 . 1 ~ .These runs were done under conditions of constant pressure following the method of Parrinello and RahmanI2 except that the box was kept rectangular in order to be able to perform constant volume simulations on the liquid at a latter time. A pressure of P = 0.5t/q3 was applied to approximately make up for the virial pressure lost due to the truncation of the potential. The mass of the wall was chosen to be the same as a single partide mass. As we shall see later the mass of the wall plays an important role in the dynamics and a light mass was used so that the wall motion would occur at a high frequency and not interfere with the low-frequency phenomena in which we are interested. In the cooling runs the sample was weakly c o ~ p l e d ' ~to, 'a~ heat bath by using the equation of motion:
-sui -
+ ii.,(t) (2) yhere Ui is the potential for the ith particle, r is the friction, and mFi =
W i ( t ) ,a Gaussian white noise term, is the random force of the (11) Stillinger, F. H.; Rahman, A. J . Chem. Phys. 1972, 57, 1281. (12) Parrinello, M . ; Rahman, A. Phys. Rev. Lett. 1980, 45, 1196. See also Andersen, H. C. J . Chem. Phys. 1980, 7 2 , 2384. (13) Schneider, T.; Stoll, E. Phys. Reu. Lett. 1978, 17, 1302. (14) Grest, G. S.; Kremer, K. Phys. Reu. A 1986, 33, 3631.
4918
The Journal of Physical Chemistry, Vol. 91, No. 19, 1987
Grest and Nagel
la2
1.1
* I-. . \
0
-6 0.0
I
I
0.2
I
0.4
I
0.8
0.6
T/E 6.0 I
I
I
I
T/E I
I
-51-0
D
1.01 0.0
I
I
I
I
I
0.2
0.4
0.6
0.8
1
!
TIE Figure 1. (a) Enthalpy H / c and (b) density pa3 vs. T / c for the initial cooling run. (c) Constant pressure, cp, and constant volume, c,, specific heats vs. temperature from the microcanonical runs. (d) Diffusion constant D vs. temperature for the constant volume microcanonical runs.
heat bath acting on each particle. The parameters for the initial cooIing run are listed in Table I, where the data were averaged over the last 487 for each temperature. After we had generated the states at each temperature by the cooling run, we made long, r = 0 and Wi = 0, microcanonical runs at both constant volume and constant pressure. For the constant pressure runs the lengths of the three sides of the box were now allowed to vary independently12but the angles were held fixed at 90°. The wall mass was left at the same value as was used in the cooling run. The constant volume runs were started from the final state obtained a t each temperature in the cooling run. Thus the volume was different from one temperature to the next and the results for the constant pressure and constant volume specific heats should refer to approximately the same state. The details for each of these simulations are summarized in Table I1 for the constant pressure runs and in Table I11 for the constant volume runs. Note that the runs in the vicinity of the glass transition temperature are quite long so that we could determine the low-frequency behavior. In order to calculate an autocorrelation function to a time t M , the total length of the run was at least 82, and often longer. The correlation function was averaged over the maximum number of initial states taken 6At apart. In Figure 1 we show various characteristics of the liquid as a function of temperature. In parts a and b of Figure 1 the enthalpy and density are shown for the initial cooling run. These results are very similar to those obtained by Wendt and Abraham4 and by Fox and A n d e r ~ e n . ~In. ~Figure I C we show both specific heats, cp and c,, for the microcanonical runs. It is quite apparent that the constant volume specific heat seems much less sensitive to the glass transition than does cp which drops by roughly 40% as the temperature is lowered through Tg. In Figure Id, the diffusion constant is shown for the constant volume runs. The diffusion
TABLE 111: Temperature T / e , Energy E / c , Average Pressure Pa'/€, Time Step A t / 7 , and Number of Time Steps M for Each System for the Constant Volume Microcanonical Runs" 1.06 0.8 1 0.59 0.52 0.40 0.3 1 0.17 0.02
-1.3714 -2.5 199 -3.6404 -4.0130 -4.6438 -5.0053 -5.4895 -5.9693
0.72 0.56 0.52 0.57 0.55 0.5 1 0.47 0.60
0.0030 0.0025 0.0025 0.0025 0.0025 0.0025 0.0030 0.0050
64 000 224 000 64 000 288 000 384 000 288 000 160 000 96 000
'While the volume was constant for each run, the density did not remain fixed as the temperature is lowered, instead the density of the system was chosen from Table I. seems to stop around T / t = 0.4, the same temperature at which there is a break in the curves for enthalpy and density. Correlation Function In these simulations we were interested in trying to see if there is a frequency-dependent specific heat near the glass transition of a Lennard-Jones liquid. We can calculate this dynamic specific heat since we know that the static value (Le., the o = 0 value) of the specific heat is simply related to the fluctuations of energy and enthalpy. We let
4dt) = (keT)-2((E(t) - E ) (E(O)- E ) ) / N
(3)
where E(?)is the total energy of the system at time t , E is the average energy and N is the number of particles. The brackets indicate an average over many starting times. By construction $ , ( t = O ) = c,(w=O) where c,, the constant volume specific heat,
The Journal of Physical Chemistry, Vol. 91, No. 19, I987 4919
Simulation of the Glass Transition
4.0 h
is measured in units of k g . We can then define
I
I
I
I
I
Constant Pressure
c f f v ( w )= w x m & ( t ) cos wt dt c f V ( w )= &(t=O)
- w l m & ( t ) sin wt dt
(4)
It is easy to show that this definition of c, obeys the KramersKronig relations and has the correct value in the static, w = 0, limit. Our simulations were under microcanonical conditions where the energy is kept constant and the temperature (which is just the average kinetic energy per particle) is allowed to fluctuate in time. Under these conditions it is also possible to calculate the specific which is given by
c(w=O)
= ( 2 / 3 - K(t=O)]-’
(5)
3.0
=-2 .-
2.0
b
8
a
I
1.o
where K(r) = N ( ( T ( t )- Q ( T ( 0 ) -
Q)/P
(6)
and where T is the average temperature. The frequency-dependent specific heat is then given by 2/3 kf(w) cf(w) =
cf’(w) =
+
[2/3
+ kf(w)]* + [Kff(w>]*
-K ”( w ) [ 2 / 3 + K f ( w ) l 2 + [k”(w)]*
0.0
. ...
t-’I/-1 W
0.0
0.1
I
I
I
0.2
0.3
0.4
0.5
tiT 4.0 R
I
I
I
I
I
Constant Volume
(7)
where
k’(w) = J m k ( tcos ) ut dt = -K(t=O)
+ w j w K ( t )sin w t dt 0
kf’(w) = J m K ( t ) sin wt dt = - w J m K ( t ) cos w t dt The constant pressure (constant volume) specific heat is obtained from the K ( t ) correlation function measured in a constant pressure (constant volume) run. K ( t ) , the kinetic energy autocorrelation function, for both constant pressure and constant volume, is shown for various temperatures by the solid lines in Figure 2. The dashed lines show the velocity autocorrelation function VAC = (i5(2t)4(0)) / ( ~ ~ ()0 )
T = 0.40 T = 0.31
(8)
At very low temperature, T / t = 0.02, there is almost no difference between the two correlation functions whereas at high temperatures they are very dissimilar. We can understand easily this low temperature behavior where the system is harmonic. In this case we can prove that k”(w) is just proportional to the density of normal modes as is the Fourier transform of the VAC. Thus at these temperatures, aside from a factor of two in the argument, the velocity and kinetic energy autocorrelation functions must be proportional to each other since they are both the Fourier transform of the same function. (The factor of two in the argument comes from the fact that K ( t ) is the product of four velocities whereas the VAC is the product of only two.) A similar argument was given for the temperature of a system after a quench.I8 At low temperatures, in the harmonic solid, there are three ways of calculating the density of normal modes, two from the kinetic energy and velocity autocorrelation functions, respectively, and a third by directly diagonalizing the dynamical matrix. We have checked that these three determinations all agree. At higher temperatures the two correlation functions behave quite differently. The VAC retains its basic shape to high temperatures: it quickly drops from 1 at t / =~ 0, goes negative, and then decays to zero by approximately t / r = 0.3. The kinetic energy curve changes much more rapidly with temperature. At (15) Lebowitz, J. L.; Percus, J. K.; Verlet, L. Phys. Reu. 1967, 153, 250. (16) Cheung, P. S. Y . Mol. Phys. 1977, 33, 519. (17) Ray, J. R.; Graben, H . W. Mol. Phys. 1981, 43, 1293. (18) Grest, G. S.; Nagel, S . R.; Rahman, A.; Witten, T. A. J . Chem. Phys. 1981, 74, 3532.
T = 0.02
0.0
0.1
0.2
0.3
0.4
0.5
tiT Figure 2. Kinetic energy autocorrelation function K ( t ) / K ( O )(solid lines) and the velocity autocorrelation function (C(2r).C(0))/(c(0)2) (dashed lines). Results for K ( t ) were obtained by averaging over the maximum number of initial states taken 6At apart. The total length of the run for each case is listed in Tables I1 and 111.
-
intermediate temperatures, near Tg, K(t) remains positive and does not decay to zero for very long times. For example at T / t 0.4 we do not observe the decay even out to t / T = 50. This indicates the presence of some very low frequency relaxation phenomena that do not couple to the velocity autocorrelation function but do couple to K ( t ) (which is the velocity squared autocorrelation function). At high temperatures K ( t ) again decays to zero fairly rapidly. We also note that the overshoot of K ( t ) through zero is washed out much more rapidly in K ( t ) than in the VAC. The oscillations that are seen in K ( t ) at constant pressure, are due, we believe, to the motion of the walls. We have checked this and have seen that the period of these oscillations varies as the wall mass is changed. The other features, such as the long decay times at intermediate temperatures, are not effected by changing the wall mass. In Figure 3 we show on a logarithmic scale the cosine Fourier transform of K ( t ) for both a constant pressure and a constant
4920 The Journal of Physical Chemistry, Vol. 91, No. 19, 1987 0.10
(a)
I
1
~
Grest and Nagel
(a)
Constant Pressure
Constant Pressure
4.0
T=0.39
0.087
I
I
I
I
0.05
0.04
0.02
0.00 -3.0
1 .o
-1.0
3.0
In W T n in
u.
0.02
t I
0.00 -3.0
I -1.0
1
1 .o
I
3.0
21
In WT Figure 3. A 32 768-point cosine Fourier transform of K ( t / 2 ) vs. In UT for (a) T / c = 0.39 at constant pressure and (b) T/c = 0.40 at constant
COT
volume volume run near T / e = 0.4. At low frequencies these curves get quite large implying that there are contributions at frequencies too low for us to measure. This behavior was not evident at lower temperatures where the peak of the Fourier transform was at WT 10 and K ( t ) went to zero for In WT < 0.1. The sharp drop in the constant pressure curve just below In WT = 3 is due to the wall motion. We have kept the wall mass light so that this feature will be at high frequency and not interfere with the data collected below In WT = 0 which is where we will look for the slow relaxation modes responsible for the frequency-dependent specific heat. We have done a simulation with a heavier wall mass and found that this drop in the transform moved to lower frequency. At high frequency we also see that the data are quite noisy. We believe that this is due to finite size effects in that the separation between adjacent normal modes is larger than the frequency width of a mode determined from the anharmonicities. We have smoothed the data in Figure 3 six times in order to reduce the amount of noise.
-
(b)
Constant Volume
I
I
I
Specific Heat
In Figure 4 we show the specific heat for a variety of temperatures for both constant pressure and constant volume on a linear frequency scale. The data has again been smoothed three times as it was for Figure 3. The upper (lower) line is the real (imaginary) part of c(w). The deep hole seen in c"&w) is again due to the wall motion. There is a clear difference between c,(w) and cP(w). As the temperature is increased there is a much greater low-frequency contribution to c,(w) than to c,(w). Indeed c,(w) seems to be basicly pushed to higher frequencies as the temperature is raised. At low temperatures C''&W) and c",,(w) both increase gradually from zero, go through a maximum near U T = 30, and then decrease to zero at high frequency. The real parts, e',(@) and c:(w), are the Kramers-Kronig transforms of this behavior. At lowfrequency they have the classical Dulong and Petit value, cp = c, = 3, and at high frequency they approach 312, the asymptotic free-particle value. Since we are simulating the entire motion of
I
0.0 0
T = 0.81
50.0
I
0.0
T = 1.06
50.0
100.0
COT Figure 4. Specific heat vs. UT for eight selected temperatures (measured in units of c) for (a) constant pressure and (b) constant volume. The upper lines in each frame are for the real part of c and the lower lines are for the imaginary part. Results were obtained by first carrying out
The Journal of Physical Chemistry, Vol. 91, No. 19, 1987 4921
Simulation of the Glass Transition
5.0
I
I
I
(a) T=0.31
II I 2
1
4.01
e
CI
2 .-’c
3.0
.-0 E
2.0
t
(a) T=0.32
g 2.0 v)
I
$1 1
1.0 -
1.o
:,I
,.” -----i ___.__
I
---------I
o.o
3.0
0
_1
0.0 -3.0
5.0
1 .o
-1.0
5.0
3.0
In WT
r;n
1
I
I
(b) T=0.39
I
(b) T=0.40 4.0 -
4.0 e
2
.-E
3.0
0
6z 2.0 1.o
0.0 -3.0
-1.0
1.o lnwT
In W T
I
I
I
(C) ..
.-0 .-E0 8
n
0
5.0
3.0
T=0.52
I
-
4.01
4.0
1
(c) T=0.52
e
i
3.0 3.0
-
2.0
-
g
v)
1.0-
0.0 -3.0 6.0
2.0
v)
_ _ - ‘-’-, ,,:‘.u* ,-
.I
L O /, -
I
-1.).
_ _ _ - - - -1.0 -, I
0.0
I
1.o In W T I
3.0
-3.0
5.0
5.0
3.0
In WT I
5.0
- --7
(d) T=0.80
5.0
1.o
-1.0
4.0
I
1
I
(d) T=0.81
“ I
2
i
i
3.0
n
1.ot
,,i
I
1
0.0
0.0
-3.0
In W T
1
-3.0
I
-1.0
----.-,-p
1.o
3.0
.
5.0
In WT
Figure 5. Real (solid line) and imaginary (dashed line) part of the specific heat c,(w) vs. In UT for four selected temperatures near the glass transition temperature. To go to lower frequencies, we carried out a 32768-point transform and used eq 7 to obtain c,(w).
Figure 6. Real (solid line) and imaginary (dashed line) part of the specific heat c,(w) vs. In W T for four selected temperatures near T,. To go to lower frequencies, we carried out a 32 768-point transform and used eq 7 to obtain c , ( w ) .
this system, our highest frequencies must be higher than the highest frequency normal mode. Thus for the T / E= 0.02 data the entire structure is due to the density of these normal modes. However, at higher temperatures we do see that there are effects
of the glass transition becoming visible in the data. At T / c = 0.39 for example, c’&u) has a vertical drop at the lowest frequency. At the same time c”Jw) rises almost vertically as well. This
4922
The Journal of Physical Chemistry, Vol. 91, No. 19, I987
behavior is seen again at T / c = 0.52 and to a smaller extent at constant volume at T / c = 0.40 and T/c = 0.52. At higher temperatures the values of c’&w) and c’Jw) decrease more gradually from the lowest frequency. We can see this behavior more clearly if we plot the specific heat on a logarithmic scale of frequency (as is done in real experiments). In Figures 5 and 6 we show the data for four temperatures around the glass transition for constant pressure and constant volume, respectively. At T / c = 0.32 the low-frequency behavior of cp is flat with only a hint of an increase at the lowest values. For T/c = 0.39 there is a visible slope to both the real and imaginary parts. Although the peak is not yet in view these data indicate that there is one somewhat below our lowest frequency. The low-frequency data are the tail of this peak. We note that the w = 0 value of the specific heat is 5.2 and our lowest frequency point only shows a value of 4.5. This indicates that there is considerable weight at lower frequencies which we are not able to see because of the limited length of our simulation. At T/c= 0.52 we see the peak in c”&w) centered approximately at In W T = -0.5. The data for c”&u) clearly are decreasing to zero as the frequency is lowered (rather than increasing as it did at lower temperatures). Finally, at T/c= 0.80 the peak has moved to even higher frequency and is beginning to merge with the peak due to the density of states of phonons. At the lowest values of frequency the data for both c’,(w) and c”&w) are almost flat, supporting our interpretation that there is little weight left at low frequency. The data for the constant volume simulations tell much the same story except that the magnitude of the effect is very much less. There is much more evidence of the glass transition in the constant pressure than in the constant volume data. Clearly this implies that most of what is happening occurs in (cp - c,) = TVa2/KT(where a and KT are the thermal expansion and the isothermal compressibility, respectively). Thus the compressibility and thermal expansion would also be important quantities to measure dynamically as we have done for the specific heat. Conclusions
We have shown that there is a large contribution to the specific heat from low-frequency modes as a liquid undergoes a glass transition. It was surprising to us that there was such a great discrepancy between the constant pressure and constant volume data. For constant volume the glass transition ,was practically invisible whereas it was quite apparent when the pressure was held constant. We also note that the wall mass does effect the dynamics in a constant pressure simulation. With this data we are now in a position to evaluate the role of simulations in understanding the glass transition. It has always been difficult to know when the simulation was no longer in
Grest and Nagcl equilibrium and thus it was somewhat problematic to understand if the simulation was an accurate description of the true physical situation. In our simulation we observe the same frequency-dependent specific heat as is seen in laboratory experiments when the sample is in equilibrium. The computer simulations are then long enough for the sample to be equilibrated when we see the frequency dependence in that data. There also seems to be corroboration to the view19 that the glass transition occurs at different temperatures depending on the cooling rate. What is important is that at each temperature the system be given time enough to undergo many oscillations of the fluctuating specific heat (or compressibility etc.). At a faster cooling rate the temperature is higher at which the time scale determined by the cooling rate does not allow for these oscillations to take place. The computer can therefore observe the same phenomena as in a real experiment except that it must be viewed on a much faster time scale. The fact that a Lennard-Jones system can show the same phenomena in its glass transition as does glycerol indicates that the frequency-dependent specific heat observed in that material is not dependent on the special nature of the glycerol molecule. On the contrary we have shown that the same effect can be seen in perhaps the simplest of all glass formers. There is a great simplification if one can study the glass transition in simple liquids rather than in complex ones. We hope that in the future it will be possible to gain even more information from these studies by investigating correlation functions that depend on spatial as well as on time variables. For example, we would hope that the function
K(q,t) =
ijl2(0) ij,*(t) expGq(?,(O)
-
?,(t))l
(9)
iJ
(where C, is the velocity of particle J at time t ) would give information about whether there is any length scale associated with the enthalpy fluctuations near the glass transition.
Acknowledgment. This work grew out of a very close collaboration with Anees Rahman that has lasted many years. We are grateful to him not only for teaching us his techniques and clarifying many of our thoughts, but also for inspiring us with his infectious enthusiasm and his delight in seeing the physics appear out of the computer simulations. We owe him a debt as well for many of the ideas that appear in this work. We also thank Norman Birge, Yoon Jeong, Paul Dixon, Gene Mazenko, David Oxtoby, and Tom Witten for many important insights and stimulating conversations. This work was supported in part by NSF Grant DMR 85-04172. Registry No. Glycerol, 56-81-5. (19) Angell, C. A . Ann. N . Y . Acad. Sci. 1981, 371, 136.
