Molecular dynamics simulations of a supercooled monatomic liquid

Jeffrey R. Fox, and Hans C. Andersen. J. Phys. Chem. , 1984, 88 (18), pp 4019–4027. DOI: 10.1021/j150662a032. Publication Date: August 1984. ACS Leg...
0 downloads 0 Views 1MB Size
Download PDF
J. Phys. Chem. 1984,88, 4019-4027

4019

Molecular Dynamics Simulations of a Supercooled Monatomic Liquid and Glass Jeffrey R. Fox? and Hans C. Andemen* Department of Chemistry, Stanford University, Stanford, California 94305 (Received: February 13, 1984)

Using a modification of the constant temperature-constant pressure molecular dynamics method, we have performed simulations of the cooling of a monatomic liquid at constant pressure to study the glass transition and relaxation processes in the low-temperature supercooled material. On cooling, the liquid has a great tendency to crystallize, but if cooled quickly enough it undergoes a transition that can be regarded as a glass transition that is broader and that occurs at higher temperature than the corresponding transition in the laboratory. The transition temperature depends on the pressure and on the rate at which the material is cooled. We present extensive results for the density, enthalpy, and self-diffusion coefficient as a function of temperature for various pressures. The supercooled material has a slow structural relaxation process that presumably falls out of equilibrium at or near the observed transition. This process leads to a small but definite history dependence of the properties of the supercooled material. The pressure dependence of the structural relaxation rate is very different from that of self-diffusion. The activation enthalpy for structural relaxation is temperature dependent.

1. Introduction Molecular dynamics and Monte Carlo computer simulations have been used to study supercooled liquids in the hope of learning about the glass transition and the structural and dynamic properties of glasses. Several model fluids, such as the hard-shere fluid,'" the Lennard-Jones and the Gaussian core fluid,1° and models of sodium silicate glasses," fused salts,12 boron tri~ x i d e , ' ~and J ~ BeF215~16 have been in~estigated.'~The behavior of simulated supercooled liquids bears some similarity to that of laboratory glass-forming liquids. On cooling, [or example, the coefficient of thermal expansion and the specific heat of the Lennard-Jones fluid fall from liquidlike values to solidlike values over a range of temperature in the supercooled liquid in a way qualitatively similar to the behavior of real supercooled fluids near their glass transition. The radial distribution function of the supercooled amorphous Lennard-Jones material is similar to experimental radial distribution functions for metallic glasses; in particular, the second peak is split into two subpeaks.8 The self-diffusion coefficient of the Gaussian core fluidlo apparently vanishes as it is supercooled to a low enough temperature. Finally, the hard-sphere fluid, when it is compressed to densities higher than that of the fluid in equilibrium with the solid, exhibits an apparent anomaly in its equation of state that may be related to a glass Much of the literature on this subject is discussed in recent reviews.18-20 A question of concern in much of this work is how can one decide if the material being simulated is a glass. A related question is how can the location of the glass transition be determined from simulations. An even more fundamental' question is whether material simulated on a computer can ever undergo any process analogous to a glass transition, considering that the time scale of computer simulations is so short compared with experimental time scales. Angell and ~ o - w o r k e r shave ~ ~ ,concluded ~~ that the short time scales of computer simulations lead to a "transition" that is shifted to such high temperatures and broadened in such a way that it bears little resemblance to the laboratory glass transition and hardly deserves to be described as a transition. Despite the limitation to relatively short times, computer simulations are still likely to be useful in developing our understanding of glasses because their use is independent of any preconceived (and hence probably incorrect) notions about the structure of glassy materials and is based only on reasonable and justifiable assumptions about the fundamental interatomic and intermolecular interactions in the material. The usual molecular dynamics and Monte Carlo methods are quite satisfactory for studying the properties of equilibrium systems. When applied to nonequilibrium systems, such as glasses, they have a number of undesirable features. The usual molecular Present address: Chemical Engineering Science Division, National Engineering Laboratory, National Bureau of Standards, Boulder, CO 80303.

dynamics method simulates the motion of a fixed number of atoms in a fixed volume. When a real material is cooled at constant pressure, however, the volume typically decreases, and it is generally acknowledged that the decrease in volume can play a significant role in causing the decreased mobility associated with the

(1) L. V. Woodcock, J . Chem. SOC.,Faraday Trans. 2,72, 1667 (1976). (2) L. V. Woodcock, J . Chem. Soc., Faraday Trans. 2, 74, 11 (1978). (3) L. V. Woodcock, Ann. N . Y. Acad. Sci., 371, 274 (1981). (4) A Rahman, M. J. Mandell, and J. P. McTague, J . Chem. Phys., 64, 1564 (1976). (5) J. P. McTague, M. J. Mandell, and A. Rahman, J . Chem. Phys., 68, 1876 (1978). (6) H. R. Wendt and F. F. Abraham, Phys. Reu. Lett., 41, 1244 (1978). (7) J. H. R. Clarke, J . Chem. SOC.,Faraday Trans. 2, 75, 1371 (1979). (8) F. F. Abraham, J . Chem. Phys., 72, 359 (1980). (9) J. R. Fox and H. C. Andersen, Ann. N . Y. Acad. Sci., 371, 123 (1981). (10) F. H. Stillinger and T. A. Weber, J . Chem. Phys., 70,4879 (1979). (11) T. F. Soules, J . Chem. Phys., 71, 4570 (1979). (12) L. V. Woodcock, C. A. Angell, and P. Cheeseman, J . Chem. Phys., 65, 1656 (1976). (13) T. F. Soules, J . Chern. Phys., 73, 4032 (1980). (14) M. Amini, S. K. Mitra, and R. W. Hockney, J. Phys. C, 14, 3689 (1981). (IS) S. Brawer, J. Chem. Phys., 72, 4264 (1980). (16) S. Brawer and M. J. Bishop, J . Chem. Phys., 75, 3522 (1981). (17) In addition to the Monte Carlo and molecular dynamics methods, several other computer algorithms have been proposed and used for generating structures of amorphous glassy materials. See, for example, D. J. A d a m and A. J. Matheson, J . Chem. Phys., 56, 1989 (1972); J. A. Barker, M. R. Hoare, and J. L. Finney, Nature (London), 257,120 (1975); G. A. N. Connell, Solid State Commun., 16, 109 (1975); L. Von Heimendahl, J . Phys. F, 5, L141 (1975); J. L. Finney, Nature (London),266,309 (1977); D. S. Boudreaux and J. M. McGregor, J . Appl. Phys., 48, 152 (1977); D. S . Boudreaux and J. M. McGregor, ibid., 48, 5057 (1977); D. S. Boudreaux, Phys. Reu E , 18,4039 (1978); K. Maeda and S. Takeuchi, Phys. Status Solidi A 49, 685 (1978); K. Maeda and S. Takeuchi, J . Phys. F, 8, L283 (1978); T. Egami, K. Maeda, and V. Vitek, Phil.Mag. A , 41, 883 (1980); S. Kobayashi, K. Maeda, and S. Takeuchi, J . Phys. SOC.Jpn., 48, 1147 (1980); S. Kobayashi, K. Maeda, and S . Takeuchi, Jpn. J . Appl. Phys., 19, 1033 (1980); T. Fujiwara and Y . Ishii, J . Phys. F, 10, 1901 (1980); R. Harris and L. J. Lewis, Phys. Reu. E , 25, 4997 (1982); F. Lancon, L. Billard, J. Laugier, and A. Chamberod, J . Phys. F, 12,259 (1982). Some of the methods create dense arrays of particles with high coordination number by adding one particle at a time to a cluster according to a well-defined prescription. Others start with such a structure or with a random array of particles, assume some form for the interatomic potentials, and perform a steepest descents minimization of the potential. In some cases, these latter methods give very different structures from those that are obtained by molecular dynamics calculations with the same potentials, indicating that the initially generated structures are dynamically metastable on a very short time scale at temperatures that allow a small amount of atomic motion and/or that the prescription for creating the initial structure generates structures atypical of those achieved by systems undergoing dynamical motion (M. Grabow and H. C. Andersen, manuscript in preparation). (18) D. Frenkel and J . P. McTague, Annu. Reu. Phys. Chem., 31, 491 (1 . 980). (Id) C. A. Angell, J. H . R. Clarke, and L. V. Woodcock, Adu. Chem. Phys., 48, 397 (1981). (20) C. A. Angell, Ann. N . Y. Acad. Sci., 371, 136 (1981).

0022-365418412088-4019$01.5010 0 1984 American Chemical Society

4020

The Journal of Physical Chemistry, Vol. 88, No. 18, 1984

onset of glass formation. There are constant pressure versions of the Monte Carlo method that get around this problem, but the atomic ”motion” in the Monte Carlo methods is not true dynamical motion, and one might wonder whether the structural changes taking place in a Monte Carlo material can mimic those in a real material, especially if the latter involve collective motion of many atoms. For dynamics simulations in which the temperature of the system must be changed, a typical procedure is to reduce the speeds of all the atoms simultaneously. It is then found that the system heats up as potential energy is converted into kinetic energy, and more energy must be removed later, leading to an oscillatory and discontinuous behavior of the temperature as a function of time. Since the properties of real glasses are dependent upon their history, it is unsatisfying to subject the simulated material to such a grossly unphysical history. Recently, we have developed a set of molecular dynamics methods that eliminates many of these problems,21and we have presented some preliminary studies of a supercooled monatomic fluid made with these r n e t h ~ d s .The ~ new methods allow us to simulate the interaction of an atomic fluid with an external heat bath and an external pressure bath, while preserving the periodic boundary conditions that must be used to eliminate surface effects. In particular, the external temperature and pressure can be varied in any desired way, and the response of the atomic fluid can be monitored. This allows us to design simulations whose protocols closely follow those of real experiments, since in laboratory experiments on glass-forming materials the experimenter usually controls the external pressure and temperature, not the volume and internal energy of the material being studied. In this paper, we present the results of simulations designed to study the behavior of a monatomic fluid as it is slowly cooled at constant pressure, with special emphasis on how the properties of the low-temperature material depend on the cooling rate and on how relaxation rates vary with temperature and density in the supercooled regime. 11. Interatomic Potential and Simulation Method These calculations used a truncated and shifted Lennard-Jones potential u(r) = 4 t [ ( u / r ) I 2- ( ~ / r )+~c]

=O

r

< rc

r>r,

where c is a constant chosen to ensure the continuity of the po~ chosen for r,. This potential tential at r = r,. A value of 1 . 5 was has the virtues of being of finite range, which increases the speed of the computation, and of having exactly the same forces as the Lennard-Jones potential for r C r,. In presenting the results, we use a set of units in which the unit of energy is t, the unit of length is u, and the unit of mass is m, the mass of an atom. The unit of time, 7,is defined as 7

=[m~~/e]‘/~

The unit of pressure is e/u3, and the unit of temperature is elk, where k is Boltzmann’s constant. The usual choices of t/ k and u for modeling argon are 120 K and 3.4 A. To model argon with the truncated and shifted potential above, the same value of u is appropriate, but a different value of should be chosen since the depth of the truncated and shifted potential is 0 . 6 8 0 ~not t . One choice would be to make the actual well depth correspond to 120 K. This gives t/k = 176.5 K. This is the value we will use. (Other more detailed arguments would give somewhat different values, but all would be larger than 120 K.) The resulting value of 7,the unit of time, is 1.770 ps. We have performed these simulations using the method of ref 21. See Appendix B for more details of the implementation of the method. The calculations were performed with a Floating Point Systems FPS-120B array processor with a PDP 11/34A minicomputer acting as a host. The inner loops of the dynamics program (force (21) H.C. Andersen, J . Chem. Phys., 72, 2384 (1980).

Fox and Andersen calculation and integration of the equations of motion) were written in APAL, the assembly language for the array processor. The number of atoms in all the simulations reported here was 500 or 512, and the equations of motion were integrated at a rate of about 100 time steps per minute. The time step used was 0.017, which corresponds to a physical time step of about 0.018 ps. With this hardware and software it is economically practical to perform simulations that are much more extensive than those previously performed on such systems. The calculations reported here are some of the longest dynamics runs ever reported. This system enabled us to investigate the effect of cooling rates that are much smaller than those used in previous studies of supercooled materials (but still they are much faster than those employed in laboratory quenches of liquid glass-forming materials). 111. Types of Simulations Performed

Several types of simulations were performed. The first is continuous isobaric cooling of the material starting from the dense liquid regime. This was implemented by equilibrating a dense liquid at a temperature and pressure in the range where the equilibrium material is a liquid and then changing the heat bath temperature to zero, keeping the external pressure fixed. The stochastic collision rate, which governs the rate at which the system exchanges energy with the surrounding bath, was adjusted to achieve the desired cooling rate. The slowest such continuous cooling runs were those that cooled the material to very low temperatures in a time of the order of 10007, about 2 ns. For continuous cooling runs, each calculation was broken into time intervals of duration t,, and, in each interval the temperature, enthalpy and volume were averaged. The second type of simulation was stepwise cooling. This was implemented by starting with an equilibrated liquid and dropping the temperature of the bath instantaneously by a finite amount, ATb, usually 0.05, and using a physically reasonable stochastic collision rate. The external temperature was held fixed for a time interval of t E , typically 607, then decreased instantaneously by a finite amount, held there for a time of t,, and then dropped again, etc. For each temperature, the configuration at the end of the tE equilibration period was used as the starting point for another dynamics run of duration tD, usually 607, during which the temperature, enthalpy, and volume were averaged. Moreover, the mean-squared displacement of particles from their starting positions were monitored during this t D data collection period and used to calculate an apparent self-diffusion coefficient by fitting a straight line to the mean squared displacement as function of time for the last half of the data collection period. These apparent self-diffusion coefficients are true self-diffusion coefficients only if the system came into equilibrium or metastable equilibrium at the temperature of the calculation during the equilibration period, tE, and if the velocity correlation function decays to zero in a time of tD or less. (The latter condition is equivalent to the statement that the mean squared displacement will continue to grow for all times a t a rate equal to the slope observed during the last half of the data collection interval.) For each of the stepwise cooling runs reported below, the system can presumably equilibrate at high enough temperatures and cannot equilibrate at low enough temperatures. The third type of simulation was “catastropic cooling”. Starting from an equilibrated liquid, the bath temperature was suddenly changed by a large amount and the calculation was continued for an “equilibration” period of tE and a data collection period of tD.

IV. Results When the material is cooled to temperatures well below the freezing point it has a great tendency to crystallize, especially for slow cooling rates. A good fraction of our computer runs had to be discarded for this reason. When crystallization occurred, it was obvious from an abrupt rise in temperature due to release of the latent heat of fusion. The presence of crystals could also be detected by analyzing the pair correlation function of the material. All the data reported below, unless otherwise stated, are for amorphous states.

Simulations of a Supercooled Monatomic Liquid and Glass TABLE I: Density p, Enthalpy/Particle H / N , and Apparent Self-Diffusion Coefficient D as a Function of Temperature T for Stepwise Isobaric Cooling at D = 0.1“ T 0 HIN D 0.5271 0.4697 0.4203 0.3768 0.3285 0.2683 0.2152 0.1690 0.1216

0.4640 0.6556 0.7416 0.7941 0.8352 0.8764 0.8997 0.9135 0.9248

-0.4466 -1.1967 -1.6339 -1.9613 -2.2595 -2.598 1 -2.831 1 -2.9990 -3.1551

The Journal of Physical Chemistry, Vol. 88, No. 18, 1984 4021 TABLE 111: Density p, Enthalpy/Particle H / N , and Apparent Self-Diffusion Coefficient D as a Function of Temperature T for Stevwise Isobaric Cooling at p = 1.0“ D T P HIN 0.7183 0.6579 0.6219 0.5598 0.4972 0.4534 0.4048 0.3501 0.3001 0.2515 0.2033 0.1487 0.0988 0.0501

0.022 30 0.01430 0.005 55 0.001 10 0.000 18 0.000 16

“ T h e parameters governing the COOhg were u = 0.2, fE = 60, to = 60, and AT, = -0.05.

TABLE 11: Density p, Enthalpy/Particle H / N , and Apparent Self-Diffusion Coefficient D as a Function of Temperature T for Stepwise Isobaric Cooling at p = 0.5” T D P HIN 0.6004 0.5522 0.5055 0.4484 0.3993 0.3520 0.3053 0.2518 0.1964 0.1531 0.1046

0.6557 0.6982 0.7454 0.7853 0.8230 0.8571 0.8834 0.9006 0.9195 0.9306 0.9430

-0.2120 -0.5570 -0.9014 -1.2046 -1.5165 -1.8 100 -2.0615 -2.2810 -2.4939 -2.6443 -2.8093

0.5991 0.5555 0.5067 0.4554 0.3995 0.3598 0.2942 0.2534 0.1958 0.1459

0.6571 0.6907 0.7382 0.7837 0.8217 0.8501 0.8851 0.9049 0.9205 0.9318

-0.2702 -0.5 126 -0.8486 -1.1866 -1.5150 -1.7528 -2.0914 -2.2984 -2.5056 -2.6710

0.7470 0.6514 0.5560 0.4580 0.3542 0.2565 0.1504 0.0499

0.5332 0.6093 0.6973 0.7808 0.8533 0.9009 0.9285 0.9500

+0.6076 +0.0642 -0.5470 -1.1656 -1.7858 -2.2622 -2.6374 -2.9568

0.018 50 0.01030 0.004 62 0.001 27

0.6799 0.7098 0.7351 0.7688 0.8083 0.8295 0.8564 0.8821 0.9041 0.9222 0.9342 0.9454 0.9552 0.9642

+0.6424 +0.3490 +0.1288 -0.1945 -0.5645 -0.7 840 -1.0558 -1.3366 -1.5869 -1.8035 -1.9750 -2.1550 -2.3 150 -2.4670

0.034 20 0.029 20 0.021 70 0.01490 0.005 79 0.002 89 0.000 54 0.000 14

” T h e parameters of the cooling process are the same as those in Table I.

TABLE IV Density p, Enthalpy/Particle H / N , and Apparent Self-Diffusion Coefficient D as a Function of Temperature Ton Stepwise Isobaric Cooling at p = 3“ T P HIN D 0.8098 0.7532 0.7060 0.6563 0.6083 0.5510 0.4922 0.4565 0.4043 0.3485 0.3057 0.2588 0.1995 0.1502 0.1003 0.0504

0.027 90 0.019 60 0.011 10 0.005 93 0.001 99 0.000 31

0.8140 0.8295 0.8434 0.8598 0.8748 0.8922 0.9120 0.9241 0.9374 0.9560 0.9644 0.9728 0.9829 0.9907 0.9982 1.0057

2.6887 2.4412 2.2253 1.9783 1.7502 1.4792 1.1832 1.0014 0.7808 0.4992 0.3449 0.1862 -0.0105 -0.1679 -0.325 2 -0.4786

0.048 50 0.032 80 0.025 10 0.022 40 0.018 70 0.009 40 0.004 69 0.002 55 0.000 99 0.000 04 0.000 20

“ T h e parameters of the cooling process are the same as those in Table I. 0.021 90 0.004 47 0.000 44

“Three separate runs. The parameters of the cooling process are the same as those in Table I, except that for the third ATb = -0.1.

Tables I-VI1 give enthalpy, density, and apparent diffusion coefficient data for the material as it is cooled isobarically at a variety of pressures. Figure 2 of ref 9 shows the density data. The density and the enthalpy are approximately linear functions of temperature at both low and high temperatures, but the slopes are smaller in magnitude at low temperature and there is a region in between, which we shall call the transition region, in which these properties have significant curvature. Figure 1 shows the density as a function of temperature in the transition region. This behavior has been observed and commented on by several a~thors.~+’J~ The transition regions are definitely in the supercooled liquid regime. The slopes of the density vs. temperature and enthalpy vs. temperature are related to the coefficient of thermal expansion and the constant pressure heat capacity, and the transition region represents a transition from a high temperature, high expansivity, high heat capacity behavior to a low temperature, low expansivity, low heat capacity behavior. Thus the “transition” is qualitatively like a glass transition, except for its width. For the time being, we use the phrases “transition” and “transition region” without intending to imply that they have anything to do with a glass transition.

TABLE V Density p as a Function of Temperature T for Continuous Isobaric Cooling at p = 6”

T

P

T

P

T

P

0.9063 0.8663 0.7537 0.6573 0.5824 0.4976

0.8983 0.9086 0.9313 0.9515 0.9679 0.9819

0.4232 0.3563 0.3072 0.2604 0.2168 0.1846

0.9925 1.0032 1.0115 1.019 1.0242 1.028

0.1585 0.1315 0.1095 0.0897 0.0754

1.0318 1.0352 1.0378 1.0397 1.0417

‘The parameters of the cooling process are Tb = 0.0, u = 0.08, and t* = 5.

TABLE VI: Density p as a Function of Temperature T for Continuous Isobaric Cooling at p = 8“ T

P

T

P

T

P

0.9123 0.8513 0.7199 0.6207 0.5488 0.4594

0.9471 0.9614 0.9779 0.9942 1.0073 1.0185

0.3866 0.3222 0.27 0.2255 0.1909 0.1584

1.0277 1.0363 1.0425 1.0478 1.0522 1.0555

0.132 0.1066 0.0865 0.071 0.058

1.0582 1.0609 1.0627 1.0643 1.0656

“ T h e parameters of the cooling process are the same as those in Table V.

The data may be analyzed more quantitatively by assuming that above and below some transition region both the enthalpy and the density are strictly linear in temperature. At each pressure, a decision was made as to what temperatures were in the transition region. The data above and below these temperatures were fitted

4022

The Journal of Physical Chemistry, Vol. 88, No. 18, 1984

TABLE VII: Density p as a Function of Temperature Ton Stepwise Isobaric Cooling at p = 10" T 0.8858 0.8477 0.7527

T 0.6769 0.5580 0.4622

4

0.9885 0.9950 1.0083

4

1.0212 1.0346 1.0478

T 0.3679 0.2604 0.2256

4

Fox and Andersen TABLE X Arrhenius Parameters for Apparent Diffusion Coefficients for Various Pressures" P A B

1.0575 1.0693 1.0730

0.1 0.5 1.o 3.0

"The parameters of the cooling process are v = 0.2, A T , = -0.05, t E = 20, and tD = 20.

P

AH

EH

1.1659 1.0890 1.0658 1.0927 1.0969 1.1249

AL

EL

TH

TL

Tt

-0.2141 -0.1960 -0.1526 -0.1235 -0.1008 -0.1091

0.9625 0.9743 1.0135 1.0511 1.0715 1.0977

0.20 0.20 0.20 0.26 0.19 0.37

0.40 0.40 0.46 0.58 0.55 0.68

0.32 0.31 0.33 0.46 0.40 0.61

OThe fits were of the form of eq 4.1. T H is the highest temperature used for the low-temperaturefit. T L is the lowest temperature used for the high-temperature fit. T , is the transition temperature calculated with eq 4.2. TABLE I X Straight Line Fits to the High- and Low-Temperature Results for Enthalpy/Particle on Cooling at Various Pressures, p" P

0.5 1.0 3.0

AH 6.3077 5.4627 4.7924

2.3073 2.278 2.4483 3.2683

"See Figures 3 and 4. Data were fit to D = A exp(-B/T).

TABLE VIII: Straight Line Fits to the High- and Low-Temperature Results for Density p for Cooling at Various Pressures, p " 0.5 -0.8478 1.0 -0.5712 3.0 -0.3133 6.0 -0.2144 8.0 -0.1646 10.0 -0.1539

5.662 2.9294 2.5383 2.9399

BH -4.0434 -3.2651 -1.1717

AI 3.3093 3.2120 3.1407

BI -3.1464 -2.6302 -0.6385

TH 0.2 0.2 0.2

TI 0.4 0.4 0.46

Tt 0.30 0.28 0.32

to straight lines by standard least-squares methods. The choice of which data to exclude is somewhat arbitrary, but at all pressures other than 0.1 a reasonable choice could be made. Typically, if too much data were included in the fit, the deviations of the data from the best straight line were of one sign at the middle of the temperature range and of the opposite sign at the two ends. By progressively eliminating more and more data in the transition region, it was possible to obtain fits that had no such systematic deviations or that agreed with the data to within its stated precision. The results are shown in Tables VI11 and IX. In each region the data were fitted to a formula of the form (4.la) at high temperature and (4.lb) at low temperature, where y is either density or enthalpy, and the coefficients depend on the property fitted, the pressure, and whether the fit is for high or low temperature. If the straight lines are extrapolated into the transition region, their point of crossing can be called the transition point.22 The temperature at which they cross is given by (4.2) The results as a function of pressure are shown in Tables VI11 and IX. The values of the A and B coefficients are smoothly varying functions of pressure, but the quotient in eq 4.2 is very sensitive to errors. The resulting Tt is a noisy but increasing function of pressure. A linear fit to the results is

+

0.1 0.5 1.o 3.0

0.1547 0.09906 0.17385 0.44099

0.592 0.397 0.50014 0.6022

1.018 1.0375 1.015 0.964

"See Figure 5. Data were fit to D = A exp[-E/(V- VD)] by a least-squares procedure by varying A , E , and VD. (See also Table XII.)

the heat capacity is slightly larger than 3, which the Dulong and Petit value for a harmonic solid in these units. The results can be used to test if the data fits the Ehrenfest equation for second-order phase transition^^^ and the Davies-Jones relation" for a glass transition in which a single ordering parameter is sufficient to characterize the excess thermodynamic properties of the liquid over those of the glass and in which the transition takes place at a certain critical value of the parameter. These equations are of the same formz5 dT/dp = TVAcy/ACP

"See Table VIII for more information.

Tt = 0 . 0 2 6 5 ~ 0.2771

TABLE XI: Doolittle Parameters for Apparent Self-Diffusion Coefficients at Various Pressures' D A R V,

(4.3)

The transition temperatures calculated from the enthalpy are equal to those obtained from the density to within the noise. The B coefficient for enthalpy is the constant pressure heat capacity. It can be seen from Table IX that at low temperatures (22) In the language appropriate for the phenomenology of glasses, the crossing point would be called the fictive temperature achieved in the limit of zero temperature.

where Aa and ACp are the differences in thermal expansion coefficient and constant pressure heat capacity for the glass and liquid. At p = 3, the results in Tables VI11 and IX and eq 4.3 give 0.0265 for the left side and 0.0340 for the right. This is reasonable agreement but no definite conclusions can be drawn. If this material has glasslike behavior, we should expect that there are slow relaxation processes in the low-temperature material and that properties should be dependent of the history of the sample as well as the temperature and pressure. Figure 1 compares density data for the p = 0.5 isobar for stepwise cooling at two different rates, stepwise heating and catastrophic cooling. The data all fall close to one smooth curve, indicating that the density is close to being a unique function of temperature and pressure and independent of history. Use of wider variety of histories did reveal a small history dependence to the density. Figure 2 of ref 9 gives an example of hysteresis in the density, with different results being obtained for the density on heating and cooling at various rates. A more detailed study was carried out to see how cooling rate affected the properties of the resulting material at low temperatures. Figure 2 shows the density achieved at zero temperature when the material is cooled isobarically at different cooling rates. At the slowest cooling rates studied, the slope of the curve corresponds to a density increase of about 0.3% for each order of magnitude decrease in cooling rate. The slowest cooling run corresponds to going from liquid temperatures to well below the glass transition in about 1 ns. If we take the liquid temperature to be of the order of lo3 K, which is appropriate for a liquid metal, this would correspond to a cooling rate of the order of 10l2K/s. Very fast liquid quenches in the laboratory are of the order of ~ one believes that the linearity can be extrapolated, lo6 K / s . ~ If it follows that a system cooled at a fast laboratory rate might achieve a density only 1.8% higher, or about 0.993. This is still lower than that of the fcc crystal at zero temperature. If this (23) H. B. Callen, "Thermodynamics", Wiley, New York, 1960. (24) R. 0. Davies and G. 0. Jones, Proc. R . SOC.London, Ser. A, 217, 26 (1953). (25) M. Goldstein, J . Chem. Phys., 39, 3369 (1963). (26) H. S . Chen, Rep. Prog. Phys., 43, 353 (1980).

S i m u l a t i o n s of a Supercooled M o n a t o m i c L i q u i d and Glass

The Journal of Physical Chemistry, Vol. 88, No. 18, 1984 4023 0.03

r

i Diffusion Coefficient

@O

B

1 Oa8

i

0.02

Density

t

0 1

a

% a

0.01

*_

8

0.7.

% ^. &

0.0

0.2

0.4

0.6

0.2

0.4

Tempe r a t u r e Figure 1. Density vs. temperature for supercooled amorphous material near the transition temperature calculated by various procedures a t a constant external pressure of 0.5. Open circles: catastrophic cooling. Closed circles: stepwise heating. Others: stepwise cooling. For all the stepwise cooling and heating runs, tE = tD = 60 and u = 0.2. For the open squares and open diamonds stepwise cooling AT, = -0.05 and for the open triangles stepwise cooling, AT, = -0.1. For the closed circle stepwise heating, AT, = 0.025. For the catastrophic cooling runs, t~ =: t~ = 100 and u = 0.2.

I

0.6

Tempe r a t u r e Figure 3. Apparent self-diffusion coefficient as a function of temperature for p = 0.5. These data are from stepwise cooling, stepwise heating, and catastrophic cooling simulations. See the caption for Figure 1 for the description of each simulation. The curve is a least-squares fit to the Arrhenius equation. See Table X.

lo-'

0

0.970

IO+

0 0

0.965

D e n s i t y a t Zero Tempe r a t u r e 0.960

1 .o

L

lo+

1

1

1

io+

10-1

Cool i n g Rate Figure 2. Density achieved a t T = 0 for continuous cooling a t p = 1 as a function of cooling rate. The abscissa is the stochastic collision frequency per particle, u . The number inside each point is the number of runs a t that cooling rate whose final densities were averaged to get the number plotted. system were actually cooled t h a t slowly, however, it would certainly crystallize. U s i n g t h e m e t h o d described in t h e previous section, w e calc u l a t e d a p p a r e n t self-diffusion coefficients as functions of t e m -

2.0

3.0

4.0

5.0

I n v e r s e Tempe r a t u r e

I

Figure 4. Arrhenius plots of the apparent self-diffusion coefficients as functions of temperature for various isobars. Open diamonds: p = 0.1, stepwise heating and cooling, one run of each. Closed squares: p = 0.5, stepwise heating (1 run), stepwise cooling (3 runs), and catastrophic cooling (1 1 runs). Open circles: p = 1.O, stepwise cooling, 2 runs. Closed triangles: p = 3.0, stepwise cooling, 1 run. The straight lines are fits to the Arrhenius equation. See Table X for the parameters of the fits. perature and pressure. T h e results a r e in Tables I-IV a n d Figures

3-6. F i g u r e 3 has diffusion coefficient d a t a a s a function of t e m p e r a t u r e for t h e p = 0.5 isobar, together with a best f i t t o t h e A r r h e n i u s equation. I n all cases, t h e A r r h e n i u s e q u a t i o n gives

Fox and Andersen

The Journal of Physical Chemistry, Vol. 88, No. 18, 1984

4024

TABLE XII: Doolittle Parameters for Apparent Self-Diffusion Coefficients at Various Pressures" 0.1 0.5 1.o 3 .O

0.2184 0.3684 0.427 0.5604

0.7491 0.8056 0.8062 0.7014

0.9979 0.99 0.9809 0.9503

'See Figure 6. Data were fit to D = A exp[-B/(V - V,)] by a least-squares procedure by varying A and B, but keeping V, equal to the zero temperature volume per particle of a perfect fcc crystal at the same pressure. (See also Table XI.)

Pair Distribution Function

6.0

10

0

I+

20

l/(V-V,) Figure 5. Doolittle plots of the apparent self-diffusion coefficient for p = 0.5. These data are from stepwise cooling, stepwise heating, and

catastrophic cooling simulations. See the caption for Figure 1 for the description of each simulation. The straight line is a least-squaresfit to the data obtained by using the Doolittle equation D = A exp[-B/(VVD)], The parameters are A = 0.09906, B = 0.397, and VD = 1.0375. Table XI contains the corresponding Doolittle parameters for other pressures.

2.01

8

1 .o

2.0

3.0

4.0

R a d i a l Distance Figure 7. Low-temperature radial distribution functions for systems prepared by continuous cooling at a constant external pressure of 1. The stochastic collision frequencies v were 3.2 for the points and 0.04 for

the 0 points. 0.02

0.01

0.7

0.e

1 .o

0.9

Dens it y Figure 6. Apparent self-diffusion coefficients as a function of temperature for various pressure. The curves are least-squares fits to the data by the Doolittle equation in which the Doolittle volume parameter V, is constrained to equal the zero temperature volume of a perfect fcc crystal at the corresponding pressure.

a good representation of the results (see Figure 4). The best-fit Arrhenius parameters for various pressures are given in Table X. In Figure 5 the data forp = 0.5 is plotted vs. l / ( V - VD), where VD is a best-fit Doolittle parameter. Table XI has the best-fit

+

Doolittle parameters as functions of pressure. There is much flexibility in the choice of parameters to fit the data. In general, the Doolittle equation provides a good representation of the data. AngellZ0 has suggested that for soft spheres the Doolittle equation describes self-diffusion data with a V, that is equal to the volume of the perfect fcc crystal. Figure 6 has fits to the Doolittle equation when the Doolittle volume is forced to be that of the perfect fcc crystal at the corresponding pressure. Table XI1 has the parameters of the fits. This form of the Doolittle equation also represents the data well. The calculated self-diffusion coefficients become zero to within calculational precision at the transition temperature or slightly below. However, the smallest nonzero result that we can calculate for simulations of this length is of the order of 10-3cr2~-'.When the values for u and 7 given above are used this smallest selfdiffusion coefficient is of the order of 7 X lo-' cm2 s-l. A real material with this self-diffusion coefficient would probably flow like a liquid on a laboratory time scale. If one trusts the Doolittle or Arrhenius fits to the data as valid extrapolations, the material would attain the viscosity of laboratory glasses only at much lower temperatures than the one at which the self-diffusion coefficient is apparently zero in our simulation. Since the volume depends on the history of the material, it is reasonable to expect that some aspects of the microscopic structure of the material are history dependent. This might appear in the pair correlation functions. Figure 7 compares the pair correlation functions for two low-temperature materials prepared at the same pressure at two very different cooling rates. There is no perceptible difference within the statistical error of the simulation. (Note

Simulations of a Supercooled Monatomic Liquid and Glass that the cusps at 1S o are a consequence of the discontinuity in slope of the potential at that distance.) The pair correlation function is too insensitive a measure of the structural changes that take place during the slow relaxation processes in low-temperature atomic liquids. Egami’s method of analysis” is a promising one for detecting and analyzing these changes. Also Voronoi polyhedron a n a l y s i ~ and ~ ~ ~bond ~ * orientational order analysisz9 may prove useful.

V. Comparison with the Phenomenology of the Glass Transition In presenting the results of the previous section, we have so far avoided giving a theoretical interpretation. In this section we discuss the results in terms of the idea that these simulations are observing a glass transition. Angell and c o - w o r k e r ~have ~ ~ convincingly ~~~ pointed out that “transitions” such as the one we observed in the computer simulation are glass transitions that have been broadened and shifted to a higher temperature than would be obtained if the corresponding experiment were performed in the laboratory. The reason for this shift and broadening is that the cooling rates for computer simulations are orders of magnitude larger than those for typical laboratory experiments on glasses. The basis for Angell’s analysis is the assumption that both the glass transition in the laboratory and the phenomenon observed in simulations is the result of a falling out of equilibrium of one or more degrees of freedom that have long relaxation times at low temperatures. (In this analysis the concept of relaxation time is used as a measure of the inverse rate of the relaxation process, with no implication that the relaxation must be exponential.) According to this analysis, the temperature at which a system falls out of equilibrium on cooling is determined by the rate at which the system is cooled and the value of the temperature-dependent relaxation time. The width of the transition is related to the temperature derivative of the relaxation time. For stepwise cooling, Angell argues that the system should be able to equilibrate at a certain temperature if it stays at that temperature for a time of at least 10 times the relaxation time at that temperature. It will not be able to equilibrate if this condition is not satisfied. For continuous cooling, Moynihan et aL30and Ritland31have derived the following relationship between the glass transition temperature and the cooling rate q: d In q/d(l/T,) = - A h / R (5.1) where q is the cooling rate in units of temperature/time, Ah is the activation enthalpy for the relaxation rate, and R is the gas constant. (In obtaining this result, it is assumed that the relaxation rate, Le., the inverse of the relaxation time, has an Arrhenius temperature dependence for the range of temperatures near the glass transition temperatures. This result can also be obtained from Angell’s analysis by taking the limit in which each stepwise drop in temperature decreases in magnitude to zero and the time spent at each temperature decreases to zero in such a way as to make the cooling rate finite. This derivation also requires assumption of an Arrhenius temperature dependence for the relaxation rate.) These two criteria allow us to deduce information about the temperature and pressure dependence of the relaxation time from the cooling rate and pressure dependence of the transition temperature. Let us use the phrase “structural relaxation” to describe the process that falls out of equilibrium and its relaxation time will be called the structural relaxation time. Angell also maintains that for atomic liquids the structural relaxation time is proportional to the viscosity coefficient and inversely proportional to the self-diffusion coefficient. From experimental results for argon and simulation results for the (27) C. S. Hsu and A. Rahman, J . Chem. Phys., 71, 4974 (1979). (28) M. Tanemura, Y . Hiwatari, H. Matsuda, T. Ogawa, N. Ogita, and A. Ueda, Prog. Theor. Phys., 58, 1079 (1977). (29) D. R. Nelson and J. Toner, Phys. Rev. B, 24, 363 (1981). (30) C. T. Moynihan, A. J. Easteal, and M. A. DeBolt, J . Am. Ceram. SOC.,59, 12 (1976). ( 3 1 ) H. N. Ritland, J . Am. Ceram. SOC.,37, 370 (1954).

The Journal of Physical Chemistry, Vol. 88, No. 18, 1984 4025 Lennard-Jones liquid, he concluded that t,D = 3.74 X

cm2

in laboratory units, where t , is the relaxation time for structural relaxation. When the values of E and (r discussed above are used this corresponds to t,D = 0.00324

(5.2)

in our units. Before comparing this analysis with our results, we should consider in more detail the question of whether the calculated apparent diffusion coefficients are true equilibrium diffusion constants. In order for this to be the case at a particular temperature, the equilibration time tE, typically 607, must be long enough for the system to equilibrate and the data collection time, tD, must be long enough for the slope of the mean squared displacement curve to have approached its final value. Whether or not the system has come to equilibrium at a given temperature can be tested by seeing if the density and enthalpy have values that are consistent with extrapolations of the higher temperature equilibrium data. Whether or not the slope of the mean squared displacement has approached its final value is more difficult to determine. For simulations of dense low-temperature liquid arg ~ n , the ~ * slope achieves its final value in a time of about 1-2 ps when the diffusion constant is 2.4 X cm2 s-l. In our units, this corresponds to a time of about 0.5-1.07 and a diffustion coefficient of 0.037$~-~.Our apparent self-diffusion coefficients have this magnitude in the equilibrium liquid regime. (For example, at p = 1.0, the freezing temperature is about 0.55,9 and D is 0.029 at that temperature and pressure.) We fit the mean square displacements to a linear function for times between 30 and 607 in extracting apparent self-diffusion coefficients. This is sixty times longer than necessary to get equilibrium diffusion constants of the order of 0.037. If the time for the mean squared displacement to achieve its ultimate slope is proportional to the structural relaxation time (and inversely proportional to the diffusion constant itself), then runs of this duration should be able to calculate true diffusion constants as small as 0.037/60 = 0.0006. To be conservative, let us triple this estimate. Thus, we believe that the apparent self-diffusion coefficients are equilibrium self-diffusion coefficients reflecting the true long-time behavior of the mean squared displacement for those states in which we have reason to believe that the system has equilibrated and in which the apparent self-diffusion coefficient is greater than about 0.002. It is clear from our simulation results that eq 5.2 relating the structural relaxation time and the self-diffusion coefficient does not hold for the system studied, even for a different choice of the constant in the equation. To see this, consider the data in Tables 11-IV for stepwise cooling at pressures between 0.5 and 3, excluding the third run at p = 0.5. The protocols of these simulations are the same for all pressures and the same values of t E and t D were used. For a pressure of 3.0, the temperature of 0.45 is well above the transition and the calculated density and enthalpy at this temperature are on the straight line fits to the high-temperature data. (Even the density and enthalpy at T = 0.4 are on the high-temperature fits, suggesting that at T = 0.4 the system is able to equilibrate on the time scale of the simulation.) Thus, there is every reason to believe that the system was able to equilibrate at T = 0.45. In this state, the apparent self-diffusion coefficient is 0.00255, which is greater than 0.002. By the criteria discussed in the preceding paragraph, therefore, the apparent self-diffusion constant is a true self-diffusion constant for this state. For p = 0.5, the density results at temperatures just above the transition indicate that the system is not able to equilibrate at T = 0.40. Since the system is not equilibrated at that temperature, it is not safe to assume that the apparent self-diffusion coefficient is a true self-diffusion coefficient. We can estimate the true self-diffusion coefficient at that temperature by fitting the results for temperature of 0.45 ( 3 2 ) A. Rahman, Phys. Rev. A , 136, 405 (1964)

4026

Fox and Andersen

The Journal of Physical Chemistry, Vol. 88, No. 18, 1984

and above to an Arrhenius form and extrapolating the formula to T = 0.40. This gives an estimated true self-diffusion coefficient of 0.012, which should be close to correct since the extrapolation was over only a small temperature interval. If the structural relaxation time were inversely proportional to the true D,we would expect that it is easier for the system to equilibrate when D = 0.012 than when D = 0.00255, but the opposite is the case. Thus the pressure dependence of the structural relaxation time is significantly different from the pressure dependence of the inverse of the true self-diffusion constant. Let us now examine the temperature dependence of diffusion and structural relaxation. We can obtain the activation enthalpy for structural relaxation by using eq 5.1, which applies to continuous cooling, and compare it with the activation enthalpy for self-diffusion. A least-squares linear fit to thep = 1 data in Figure 2 for the smaller cooling rates gives dpo/d In q = -0.00149 where p o is the density achieved at zero temperature. This is equivalent to EL in eq 4.lb when the property y is density. If we assume that the coefficient of thermal expansion of the low-temperature glassy material, Le., the parameter AL in (4.lb), is independent of the cooling rate,33then we can find the cooling rate dependence of the transition temperature by differentiating eq 4.2 with regard to In q. We find dT,/d In q = (dBL/d In q ) / ( A H- AL) Using this equation, the previous one, and the results in Table VI1 for p = 1 we find d In q/d(l/T,) = -23.51 This implies an activation enthalpy of 23.51 for the structural relaxation. For comparison, the activation enthalpy for selfdiffusion at that pressure is 2.10. (This result was obtained by fitting the data for the high-temperature range in which the calculated diffusion coefficients are known to be equilibrium constants.) This indicates that the structural relaxation rate is a much more rapidly varying function of temperature than the self-diffusion constant. The difference in the temperature dependence of these two relaxation processes might be explained by the fact that the activation enthalpy for diffusion is calculated from high-temperature equilibrium states while the calculated activation enthalpy for structural relaxation applies to the lower temperature states in the glass transition region itself. For most glass-forming materials, activation enthalpies for transport processes are temperature dependent and increase in magnitude as the temperature is lowered.I9 Moreover, the curvature in Figure 2 implies that the activation enthalpy for structural relaxation also increases as the temperature is lowered. Thus, it is possible that the temperature dependence of both the structural relaxation rate and the selfdiffusion process are governed by the same temperature-dependent activation enthalpy. This suggests that a relationship of the form of eq 5.2 might hold for a range of temperatures for a given pressure, but the “constant” on the right side of the equation would have to be a pressure-dependent quantity. Thus the relaxation process that falls out of equilibrium at the transition has a pressure dependence different from that of the self-diffusion coefficient but may have the same temperature dependence. The structural relaxation process has an activation enthalpy that is temperature dependent and increases as the temperature is decreased.

barically cooled through the glass transition. The system studied was a model of a monatomic fluid. The material had a great tendency to crystallize, but if cooled rapidly enough it undergoes what appears to be a broadened glass transition. Hysteresis in the transition is evident from a comparison of heating and cooling runs. Slowly cooled materials achieve a higher density than the more rapidly cooled materials. These results indicate that the material has a slow relaxation process that involves volume changes. The apparent self-diffusion coefficient appears to be well described by both the Arrhenius equation and the Doolittle equation. The volume appearing in the denominator of the exponent of the Doolittle equation is approximately equal to the volume of the fcc crystal at the corresponding pressure, but because of the uncertainty in our results for small values of the self-diffusion coefficients we cannot state the precise relationship between that volume and the crystal volume. The relaxation process that falls out of equilibrium at the glass transition has a very different pressure dependence than does the inverse of the self-diffusion coefficient. The structural relaxation process has an activation enthalpy that is temperature dependent and increases as the temperature is decreased. The question arises of whether we are really seeing a glass transition. In our opinion, the essential feature of the glass transition is that the system crosses from one regime in which the system can achieve internal equilibrium on the experimental time scale to another regime in which it cannot do so. The time scale of our computer experiments is much shorter than that of laboratory experiments on glasses, but this does not invalidate the description of our calculation as a simulation of the glass transition. The important question is whether the degrees of freedom that fall out of equilibrium in the simulation are the same as those that fall out of equilibrium in the laboratory system. We do not yet know what that degree of freedom is, but it is clear from the density data that the material has some property that is history dependent and that is apparently falling out of equilibrium in the transition region. Acknowledgment. This work was supported by the National Science Foundation (Grant C H E 81-07165) and the Center for Materials Research at Stanford University. Appendix A. Algorithm for Integrating Equations of Motion Let ri and p i , i = 1, ..., N, denote the positions and momenta

of the N atoms in the sample being simulated, and let V be the volume of the cubic sample. The constant pressure molecular dynamics method2’ establishes a correspondence between these variables and those of a scaled system according to V=Q

ri = Q’13pi pi = iri/Q1I3

(33) We have not tested this assumption quantitatively for the present system. For analogous calculations for a supercooled Lennard-Jones mixture, this assumption is very accurate (M. Grabow and H. C. Andersen, unpublished calculations).

..., N

i = 1, ..., N

The Lagrangian and Hamiltonian and the corresponding equations of motion for the scaled system were presented in ref 21. The method we used to integrate the equations of motion is based on the velocity form of the Verlet algorithm.34 To introduce the method, let us first consider the velocity form. Consider a differential equation of the form X(t)

= AX(?)]

(A.1)

The velocity form of the Verlet algorithm is based on the following equations:

+ h ) = x ( t ) + hX(t) + h2Ax(t)]/2 + h ) = X(t) + (1/2)hWx(t)l + A x ( t + h)11 x(t

VI. Conclusions and Summary

We have performed constant pressure molecular dynamics simulations that, except for the magnitude of the cooling rates, closely mimic laboratory experiments in which a liquid is iso-

i = 1,

X(t

(A.2a) (A.2b)

where h is the time step in the numerical integration. These equations have errors of order h3. When they are used as an iterative procedure, with the result of the first equation being (34) W. C. Swope, H. C. Andersen, P. H. Behrens, and K . R. Wilson, J . Chem. Phys., 76, 637 (1982).

The Journal ofphysical Chemistry. Vol. 88, No, 18, 1984 4027

Simulations of a Supercooled Monatomic Liquid and Glass

substituted into the right side of the second, their error is still of order h3. These equations can be used to calculate x and 2 at time t h from the values at time t. The global error, i.e., the error made in integrating the equations for a finite length of time, is of order hZ,just as in the usual form of the Verlet alg~rithm.’~ In the present case, this method cannot be used because the accelerations depend explicitly on the velocities as well as the coordinates. Instead of (A.l), consider coupled equations of the form = Ax(t),x(t),v(t),$(t)l (A.3a)

+

.w)

J(t) = g[x(t)At)dt)l

(A.3b)

Here x corresponds to the p variables and y corresponds to Q. Note that the acceleration of y does not depend on the velocity of y . Such equations can be solved with the following series of equations: x(t

+ h ) = x ( t ) + hx(t) + ( 1 / 2 ) h z f l x ( t ) , a ( t ) , y ( t ) , $ ( ? ) ] (A.4a)

$*pp(‘ + h ) = $ ( t ) + (1/2)(g[x(t),n(t),y(l)l + g[x(t + h),x(t),Y(t + h)lJ (A.4c) i(t

$(t

+ h ) = x ( t ) + (1/2)hVTx(t),x(t),y(t),$(?)l + Ax(t + h ) , i ( t + h),y(t),$app ( t + h)lI

(A*4d)

+ h) = $(t) + (l/z)h(g[x(t),n(t),r(t)l + g [ x ( t + h M ( t + h)9Y(t + h)lJ (A.4e)

Equations A.4a, A.4b, A.4d, and A.4e are similar to eq A.2a and A.2b and contain errors of order h3,provided that in eq A.4d the quantityqaPp(t h) is interpreted as the correct y(t h). Equation A.4c gives an approximation for p ( t h ) that has errors of order h2. These equations can be used to solve the equations of motion. If x ( t ) , y ( t ) ,k(t), and j ( t ) are known, then the first two equations give x ( t h ) and y ( t + h ) , with errors of order h3. When they are substituted on the right side of (A.4c) the resulting japp(t h ) has errors of order h2,of the same order as if the exactly correct x(t h ) and T ( t + h ) were used. Then (A.4d) gives x(t + h ) with an error of order h3, the hZ error in $app contributing errors of higher order. Finally, (A.4e) gives an improved value of $(t h), with errors of order h’. The error made on each time step is the same as that of the velocity form of the Verlet algorithm, and so the global error will also be the same, namely, of order

+

+

+

+

+

+

+

h2.

+

Note that eq A.4d contains x(t h ) on both sides of the equation and is an implicit equation for n(t h). Thus, this method is not a general one for equations of the form of (A.3). In the present case, however, the equation is a linear one and can easily be solved explicitly.

+

Appendix B. Details of the Constant Temperature-Constant Pressure Method In the constant temperature molecular dynamics method:’ the particles are subjected to stochastic collisions. We used that method in the present calculations. A parameter of the calculation is the stochastic collision frequency, u, which is the average number of stochastic collisions suffered by a particle per unit time. The (35) G. Dahlquist and A. Bjorck, “Numerical Methods”, Prentice-Hall, Englewood Cliffs, 1974, p 353.

product vhN is the average number of particles that suffer a collision per time step, h , of the numerical integration. This number is typically less than one and can be interpreted as the probability that some particle will suffer a stochastic collision in a time interval equal to the time step. At each time step, a random number uniformly distributed between 0 and 1 is compared with uhN. If the random number is less than uhN, then a Stochastic collision is performed by choosing a particle at random and assigning new velocity components to it that were drawn at random from a Boltzmann distribution at the temperature of the heat bath. [Each velocity component is a randomly signed product of (2kTb/m)1/2and the inverse error function of a random number uniformly distributed on the interval from 0 to 1.1 The stochastic collision at time t , if any, is performed after the positions and velocities at time t are obtained and before any of the positions and velocities at time t + h are calculated. In the constant pressure molecular dynamics method,21a value of the parameter M must be chosen. M represents the inertia associated with changes of volume of the system and hence can be regarded as the mass associated with a piston that moves when the system expands or contracts. (The piston is not a physical one, however, since the expansions and contractions are isotropic and preserve periodic boundary conditions.) The value of M is irrelevant for equilibrium simulations, but for nonequilibrium simulations M affects the frequency of oscillations of the volume. Sustained oscillations of the volume represent sound waves generated by the nonequilibrium process that are trapped in the sample because of periodic boundary conditions. In a volume element of an infinite fluid, such sound waves would propagate away on a time scale on the order of the length of the volume element divided by the speed of sound. This behavior can be achieved in our periodic boundary condition simulations by choosing the parameter M so that the oscillations have the desired period (namely of the order of the length of the box divided by the speed of sound) and subjecting the velocity of the volume to stochastic collisions so as to damp the oscillations on the same time scale. A series of simulations at p = 1 and T = 1 were performed to determine oscillation frequency for the volume as a function of M . From the results, we concluded that M = 0.002 m oW4was appropriate for an assumed speed of sound of the order of 4u/r (approximately lo3 m s-’). The parameters governing the stochastic collision of the volume velocity were chosen so that the oscillations are damped on the same time scale as their period. Stochastic collisjons affecting the rate of change of the volume (Le., the variable Q)are implemented in a manner similar to the Metropolis algorithm for Monte Carlo calculations. A trial value Q is calculated by adding to Q a random number uniformly distributed in the interval from -R to +R, where R is a parameter of the stochastic collisions. If the trial value of Q is smaller in magnitude than the original value, then the trial value is accepted as the new value. If the trial value is larger in magnitude than the original value, then the kinetic energy ( M Q 2 / 2 )associated with the trial value is larger than that of the original value, in which case the trial value is accepted with a probability of exp( - A E / k T * ) , where AE is the increase in kinetic energy and T* is the system temperature modified to included the volume as a dynamical variable: kT* = kT[3N/(3N

+ l ) ] + (1/2)Q2M/(3N + 1)

The stochastic collisions of the piston volume are implemented on every time step with R = 1. This leads to a Brownian motion of the piston, with at most small changes of the piston velocity on each time step.