J. Phys. Chem. 1986, 90, 203-205
203
Constrained Dynamics and Its Application to the Isothermal-Isobaric Simulation of the Supercooled-LiquidIGlass Transition Region Farid F. Abraham IBM Research Laboratory, San Jose, California 951 93 (Received: June 17, 1985; I n Final Form: September 16. 1985)
Constrained dynamics is reviewed and generalized so that specified temporal variations of temperature and/or pressure may be implemented in the simulation of classical many-particle systems. This technique is applied to the study of glass formation of a rapidly quenched Lennard-Jones liquid at a constant and low pressure. The glass transition temperature is found to be approximately one-half the melting temperature for low external pressure.
Introduction We have presented an analysis for isothermal-isobaric Monte Carlo simulations of metastable supercooled liquid and amorphous phases of Lennard-Jonesium and noted several distinct signatures for identifying the glass transition boundary; Le., the density, enthalpy, and pair distribution function dependences on temperature and pressure are different for the two phases.’ A recent theory of the glass transition requires computer-experimental determination of the transition temperature at low pressure.2 Therefore, we have extended our earlier study to include constant-pressure quench experiments at the Lennard-Jones reduced critical point pressure of 0.17. Furthermore, the simulations used a different formulation of conventional molecular dynamic^,^ termed constrained dynamics,+l* for studying isothermal-isobaric processes. For this present study, an extension to constrained dynamics has been invented and is described in the next section. Application of this technique to the Lennard-Jones glass study yields a glass-transition temperature of approximately one-half the melting temperature for an external pressure equal to the critical point pressure.
The Constrained Dynamics Method The molecular dynamics simulation technique yields the motion of a given number of atoms governed by their mutual interatomic interactions, this being calculated by numerical integrations of Newton’s equations of m ~ t i o n : ~ , ~ dri/dt = pi/m
(1)
dpi/dt = Fi
(2)
In the traditional molecular dynamics experiment, the total energy E for a fixed number of atoms N in a fixed volume Vis conserved as the dynamics of the system evolves in time, and the time average of any property is an approximate measure of the microcanonical ensemble average of that property for a thermodynamic state of N,V,E. For certain investigations, it may be advantageous to perform the simulation at constant temperature and/or pressure. A very clever procedure for simulating a constant temperature system has been invented by Hoover6 and, independently, by Evans.’ Taking the relation p?/2m = 3/,NkT
(3)
I
(1) F. F. Abraham, J . Chem. Phys., 72, 359 (1980).
(2) Alf Sjolander, Chalmers University, private communication. (3) A. Rahman, Phys. Rev. A, 136, 405 (1964). (4) F. F. Abraham J. Vac. Sci. Technol., E , 2, 534 (1984). The reader
is referred to the bibliography at the conclusion of this paper. (5) W. G. Hoover, A. J. C. Ladd, R. B. Hickman, and B. L. Holian, Phys. Rev. A, 21, 1756 (1980). (6) W. G. Hoover, A. J. C. Ladd, and B. Moran, Phys. Rev. Lett., 48, 18 18 (1982). (7) D. J. Evans, J . Chem. Phys., 78,3297 (1983). (8) D. J. Evans and G. P. Morriss, Chem. Phys., 77,63 (1983). (9) F. F. Abraham, J . Vac. Sci. Techno/.,E , 2, 534 (1984). (10) D. J. Evans and G. P. Morriss, Phys. Lett., 98A,433 (1983). (11) D. J. Evans and G. P. Morriss, Compt. Phys. Rep., 1, 297 (1984).
0022-3654/86/2090-0203$01 .50/0
one can replace the equation of motion (2) by dpi/dt = Fi - (pi
(4) with the constraint that the time derivative of eq 3 equals zero. Therefore
s’ = (E F,Pi)/C Piz i
i
(5)
A natural extension to Hoover’s novel approach to molecular dynamicsS has been made independently by Evans and Morrisss and by Abraham9 in order to simulate a constant-pressure process. Equations 1 and 2 are replaced byS dri/dt = (pi/m) $ri (6)
+
dpi/dt = Fi - @pi where, in three dimensions 1 dV $ = - 3V dt
(7)
In the spirit of Hoover’s approach of imposing the constant-temperature constraint, we also impose the constraint that the time derivature of the pressure given by the virial expression PV =
y3c(p?/m + riFi) i
(9)
equals zero; i.e., dP/dt = 0. Similar to obtaining an expression for the constant-temperature coupling parameter {, we obtain the following expression for i: in terms of the instantaneous positions, momenta, and derivatives of the forces: i: = Q N / Q D
(10)
where +ij is the two-body central potential between the ith and j t h particles separated by the distance rij = Ir, - rjl. The order of the derivative of $ij with respect to rij is denoted by the superscript primes. Evans and Morriss have derived exact expressions for the equilibrium distribution functions generated by eq 6-12.1° A useful summary of constrained molecular dynamics is given in ref 11. For this glass formation study, we have extended the molecular dynamics technique to account for an imposed time variation of temperature and/or pressure. This extension allows us to do quench, crush, and crunch experiments.’ This is achieved by taking the time derivative of eq 3 and 9 and by not setting dT/dt = p and dP/dt = P equal to zero. The equations of motion are given by eq 6 and I where
0 1986 American Chemical Society
204
The Journal of Physical Chemistry. Vol. 90, No. 1, 1986
Abraham
7.00 6.00 5.00 4.00 3.00 2.00
(P", T") = (0.17, 0.01)
1 .oo 0.00
t
-0.5
'0
E 0.2
0.4
0.6
0.8
1.0
Temperature, kT/E Figure 2. Time-averaged equilibrium density and enthalpy as a function of temperature quench.
i
0.00
.
_
2.00
U
4.00
Distance, r/u Figure 1. The pair distribution functions g(r) as a function of radial distance r for a sequence of temperature quenches initiated from the liquid state (P*,T*) = (0.17,l.O).
QN and QD are given by eq 11 and 12. The extension to two dimensions is straightf~rward.~ The temperature is initialized by normalizing the particle velocities so that the mean kinetic energy corresponds to the desired initial temperature To.From an arbitrary initial configuration (e.g., a crystal lattice), the pressure P, is calculated from the virial expression (9) and changed to the desired pressure Po by using the dynamical eq 6 to 8 and 13 and 14 and specifying P ( t ) over time interval 7 such that Po = P,
+
L'
P ( t ) dt
Of course, for a constant pressure simulation at p0, P ( t ) is zero for t 2 7. The numerical integration of the dynamical equations was implemented by using the IBM Dynamic Simulation Language/VS, and the variable-step, variable-order Adams-Moulton integration method was chosen.I2
The Glass Experiments The interatomic potential was chosen to be Lennard-Jones 12-6 with (€,a)denoting the well-depth and size parameters, respectively. In order to simulate the bulk, the standard periodic boundary conditions were imposed with respect to translations parallel to the faces of the computational cube composed to 512 atoms. In each simulation, we started with the last fluid configuration of a well-equilibrated fluid at reduced temperature T* = kT/t = 1.0 and reduced pressure P* = Pu3/t = 0.17 (critical point pressure). Then, using our extension of constrained dynamics, we quenched to the new (T*,P* = 0.17) using a linear temporal decrease of temperature over a time span of 100 6 t * , where 6t* = 6t (e/mu2)'/2 = 0.002. After the Lennard-Jones system equilibrated ( ~ ( 3 - I O )X lo3 at*, depending on the extent of the quench), (5-10) X lo3 6t* were simulated in order to obtain good statistical averages of the quantities of interest. As illustrated in Figure 1, the pair distribution functions g ( r ) show the general features expected. In the liquid and supercooled-liquid state, two smooth peaks exist: a first prominent peak and a second smaller peak corresponding to the first and second coordination shells of an atom in the liquid, respectively. As the temperature quenches probe deeper into the metastable region, the first peak becomes more pronounced in magnitude and narrower in width, the first minimum decreases in magnitude, and (12) IBM Dynamic Simulation Language/VS, Language Reference Manual, Program No. 5798-PXJ, IBM Corporation Department, G12/141, San Jose, CA 95193.
J. Phys. Chem. 1986, 90, 205-206 the second peak gradually flattens in shape with an eventual ”bimodal splitting” at very low temperatures (T* 5 0.4). The development of the split second peak indicates an amorphous atomic packing, this being the principal structural feature in the experimental PDF of amorphous materials that previous theoretical models have attempted to describe. In Figure 2, we present the time-averaged density p and enthalpy H for the quench experiments. As first noted in our earlier study,’ the dependences of density and enthalpy over significant portions of the liquid/supercooled-liquid branch and “what we now call” the amorphous branch of the (P*,T*) state variables
205
are linear, the slopes of the two branches are different, and the intersections by linear extrapolation of the respective branches are essentially identical; Le., for the temperature quench experiments the intersections for p, H are at Tg*= 0.395-0.4. We identify this point as the “glass transition point”. The discontinuity in slope of the density variation at the glass transition temperature (for a given quench rate) is well-known for laboratory g1a~ses.l~ (1 3) D. R. Uhlmann and R. W. Hopper in “Metallic Glasses”, J. J. Gilman and H. J. Leamy, Ed., American Society for Metals, Metals Park, OH, 1978, Chapter 6.
Dilatancy Associated with One-Dimensional Compression of Materials George H. Duffey Department of Physics, South Dakota State University, Brookings, South Dakota 57007 (Received: July 18, 1985)
Just as the packing of the molecules in a fluid is disrupted by shear, this packing is also disrupted by one-dimensional compression. Evidence for this effect comes from the apparent attenuation of the von Neumann spike in detonation waves. Using the dilatancy as an additional independent variable, we have developed some of the relevant thermodynamics.
In thermodynamic discussions, systems are characterized by Compression waves traveling at nearly constant velocities are macroscopic variables. Relationships among these are deduced observed in detonating explosives. A typical wave consists of a from the energy considerations embodied in the zeroth, first, shock front, a zone in which the reaction supporting the front second, and third laws. The classic theory treats equilibrium proceeds, a Chapman-Jouguet surface, and a continually exsystems and reversible processes.’-3 The newer developing theory panding rarefaction zone.8 When the wave through the reaction deals with nonequilibrium systems and irreversible p r o c e s s e ~ . ~ ~ ~zone is steady, mass, momentum, and energy do not accumulate In this irreversible thermodynamics, the familiar state quantities in any intermediate region. Then the Newtonian force and conservation laws yield simultaneous algebraic equations. When vary with position; but in each small region they are usually the pertinent material properties, the wave velocity, and conditions considered to be related as in an equilibrium system. However, should such relationships not be affected by the in front of the wave are specified, these equations yield the conmacroscopic processes that have occurred and are occurring in ditions at any stage in the steady-state region. Introducing the the element of volume? Thus, Hanley and Evans6%’have conChapman-Jouguet condition for the end of the reaction zone then sidered the effects of a constant shear rate on a fluid system. They determines the detonation velocity. found that the shearing disrupts packing of the molecules, lowering In the crude approximation that the number of moles of gas the density. The effect is called shear dilatancy. per unit mass, the covolume, and the heat capacity ratio are Now, a rapid compression, as we see in a shock wave, should constant, and that the wave is planar, one finds that the longiproduce a similar effect, a compressive dilatancy. For, each tudinal pressure just behind the shock front, where no reaction element of fluid passing through the shock front is compressed has yet occurred, is twice that at the Chapman-Jouguet point. longitudinally, but not transversely. About a typical molecular The corresponding peak in conventional pressure is called the von center, the radial distribution function for density moves to its Neumann spike.9 Introducing curvature and introducing nonhigh pressure form in the longitudinal direction, while in the steady behavior have little effect on the behind-shock-front transverse directions little change takes place. During the compressure. On the other hand, these perturbations have a conpression, the molecular units do not have time to move around siderable effect on the conditions at the Chapman-Jouguet point, and pack in the transverse directions. But as the element moves as is well-known.8 So if local equilibrium did exist through the supporting zone, through the supporting zone following the shock front, this rethe density would be considerably higher in the early part of the arrangement process proceeds and the transverse packing increases reaction zone. For common explosives, however, measurements to the equilibrium value asymptotically. of density there give results that do not differ much from those at points later in the zone.9 But one cannot repeal Newton’s laws. (1) G. H. Duffey, “Physical Chemistry”, McGraw-Hill, New York, 1962, The only way out of the dilemma is to admit that a compressive pp 185-346. dilatancy analogous to the established shear dilatancy is present. (2) F. W. Sears,‘An Introduction to Thermodynamics, the Kinetic Theory Let us consider an explosive undergoing one-dimensional comof Gases, and Statistical Mechanics”, 2nd ed,Addision-Wesley, Reading, MA, 1953, pp 1-199. pression and subsequent reaction in an apparently steady manner. (3) D. ter Haar and H. Wergeland, “Elements of Thermodynamics”, AdThe fore part of the wave can then be used as a reference structure, dison-Wesley, Reading, MA, 1966, pp 1-132. attaching to it a coordinate system with axis 1 pointing normal (4) K. G. Denbigh, “The Thermodynamics of the Steady State”, Wiley, New York, 1951, pp 1-99. (5) I. Prigogine, “Introduction to Thermodynamics of Irreversible Processes”, 2nd ed, Interscience, New York, 1961, pp 3-115. (6) H. J. M. Hanley and D. J. Evans, J . Chem. Phys., 76, 3225-3232 (1982). (7) D. J. Evans, J . Chem. Phys., 78, 3297-3302 (1983).
0022-3654/86/2090-0205$01.50/0
0 1986 American Chemical Society
