10628
J. Phys. Chem. B 2008, 112, 10628–10637
On the Role of Inherent Structures in Glass-forming Materials: II. Reconstruction of the Mean Square Displacement by Rigorous Lifting of the Inherent Structure Dynamics Dimitrios Tsalikis,† Nikolaos Lempesis,† Georgios C. Boulougouris,‡ and Doros N. Theodorou*,† School of Chemical Engineering, National Technical UniVersity of Athens, Zografou Campus, GR-15780 Athens, Greece, and Scienomics, 17, Square Edouard VII, 75009 Paris ReceiVed: February 14, 2008; ReVised Manuscript ReceiVed: June 07, 2008
In a previous paper, we investigated the role of inherent structures in the vitrification process of glass-forming materials, showing that the dynamical transitions between inherent structures (states) can be well predicted by a first-order kinetic scheme based on infrequent-event theory at low temperatures. In this work, we utilize and extend that methodology in order to completely reconstruct the system dynamics in the form of the mean square atomic displacement as a function of time at finite temperatures on the basis of the succession of transitions in a network of states, the vibrational contribution being evaluated on the basis of short molecular dynamics runs artificially trapped within each one of the states. In order to do so, we provide the mathematical formulation for lifting the coarse-grained Poisson process model of transitions between states back to the atomistic level and thereby reproducing the full dynamics of the atomistic system within the Poisson approximation. Our result shows excellent agreement for temperatures around and below the glass transition temperature of our model Lennard-Jones two-component mixtures. Clearly, our approach is able to reproduce the full dynamics of the atomistic system at low temperatures, where the Poisson approximation is valid. Introduction Glassy materials have acquired an important role in our life and consequently attract the interest of the scientific community in both applied and basic research. Although theoretical and experimental studies of the glass transition and the glassy state have been a focus of research for many years in the field of materials science, many basic questions still remain unanswered, especially regarding the connection of macroscopically observed properties to atomic-level structure and dynamics. From a thermodynamic point of view, materials in the glassy state have not reached equilibrium, and therefore, their properties depend on their formation history. From the dynamical point of view, the longest relaxation time exceeds macroscopic observation times below the glass transition temperature Tg. The study of glassy systems by means of molecular simulation faces serious challenges, because one needs to bridge time scales spanning some 20 orders of magnitude, from the period of fast atomic vibrations (10-14 s) up to the longest time for structural, volume, and enthalpy relaxation (on the order of years at 20 °C or so below Tg). The transition from the ergodic liquid to the nonergodic glass does not appear to be discontinuous. Initially, the liquid, when cooled below the melting point, reaches a metastable state described as supercooled liquid. The supercooled state is metastable with respect to crystallization for systems where molecular geometry allows packing into regular lattices. Nevertheless, it can be considered as an equilibrium liquid for most practical purposes. One may use traditional sampling techniques to calculate the thermodynamic properties of this state, because the time needed to explore all relevant microscopic states which shape the properties of the supercooled liquid is much shorter * Author to whom correspondence should be addressed. Fax: +30 210 772 3112. E-mail:
[email protected]. † National Technical University of Athens. ‡ Scienomics.
than the observation time in both experiment and simulation far above Tg. As the system is cooled toward Tg, the necessary time for establishment of metastable equilibrium increases dramatically. As we have shown in our previous work,1 the necessary time for an equilibrium sampling of liquid states increases beyond the observation window of simulation, leading to strong nonequilibrium behavior. Over the past decades, molecular simulation has proved to be a reliable tool for predicting the thermodynamic and mechanical properties of matter and has provided new insights into the molecular mechanisms underlying macroscopic properties. For glass-forming materials, considerable progress has been made with the help of molecular simulations in the area of supercooled liquids.2 Unfortunately, the broad range of time scales for molecular motion present in glassy systems poses severe limitations for molecular simulation in the vicinity of and below Tg. Any traditional numerical solution of the time evolution equations of a microscopic particulate model is bound by the fastest process present, thereby limiting the ability to track the time evolution out to the desired longest time scale. As a result, brute force molecular dynamics (MD) simulations are doomed to describe only a very short part of the spectrum of time scales characterizing motion in a glassy system and therefore cannot answer questions concerning the long time relaxations present in glasses. On the other hand, different vitrification conditions lead to glasses with different properties. For example, in the case of vitrification via cooling, the faster the cooling rate, the closer the resulting glass is to the higher temperature conditions of the metastable liquid from which it is created, underlining the importance of the nonequilibrium character of the vitrification process.1 State-of-the-art theories of the supercooled state include mode-coupling theory3,4 and theories for enumerating the number of stationary points,5 which focus on the topography of the multidimensional energy hypersurface of the system.6 These theories try to connect the phenomenon of glass transition
10.1021/jp8013223 CCC: $40.75 2008 American Chemical Society Published on Web 08/01/2008
Reconstruction of the msd by Rigorous with the thermodynamic arrest of the system in the vicinity of some local minima in its energy landscape.7,8 In a more general sense, analyses of the potential energy landscape have been used in the past in both the supercooled region2,9 and in the deep glassy state.10,11 For a comprehensive recent review of landscape ideas the reader is referred to ref 9. One of the most important features of the potential energy landscape is the local minima of the energy, or inherent structures,12 around which the system is expected to spend most of its time trapped, at least at low temperatures. Below Tg, the importance of inherent structures to molecular simulation has been extensively explored.9–11 Far below Tg, Boulougouris and Theodorou10 were able to simulate the dynamics of an atomistic model of glassy atactic polystyrene over more than 10 orders of magnitude on the time scale. To do so, they developed an alternative to kinetic Monte Carlo sampling, the dynamical integration over a Markovian web (DIMW) methodology,10 which samples in parallel all possible paths for evolution in time and enables the direct evaluation of the time-dependent probabilities of occupancy of states for a system undergoing successive transitions with firstorder kinetics in a landscape of infinite extent. On the basis of DIMW, it is possible to construct an ever-expanding network of known, or explored, states, bounded by a set of boundary states, by starting from an initial (small) set of states. For each explored state, all transitions connecting it with its neighboring states have been located, and the corresponding transition rate constants have been computed by atomistic infrequent-event analysis. Boundary states are connected to explored states but are not yet explored themselves. The time-dependent probability distribution among explored states is determined via analytical solution of the master equation for the explored states under absorbing boundary conditions for the boundary states. The set of explored states is expanded systematically whenever necessary through a stochastic scheme that each time selects to include in the set of explored states a boundary state, according to the probability flux to that state. In their work, Boulougouris and Theodorou10 used multidimensional transition state theory within the harmonic approximation in combination with the dimer saddle-point search method13 in order to locate and evaluate transitions and rate constants out of the states being explored on the fly. Note that with DIMW, one is able to follow the time evolution of the probability density itself and therefore effectively track all possible trajectories and not a single one, as with MD or traditional kinetic Monte Carlo schemes. DIMW can be also viewed as the extension to kinetic Monte Carlo of the integration over a Markovian web (IMW) method, introduced by Boulougouris and Frenkel,14 which, via a rigorous statistical mechanics formulation, allows the enrichment of ensemble averages and the combination of multiple integration levels in ensemble averaging, including combination of dynamics and static integration in Monte Carlo. In our previous work,1 we have investigated the significance of the inherent structure concept for the prediction of the dynamical behavior of glass-forming materials. To do so, we observed the mean squared displacement (msd) of atoms in the system, as obtained from full-blown MD simulations and from the corresponding time evolution of the inherent structures. We found that, for temperatures significantly higher than Tg, the dynamics of the system can be captured well by the hopping process between regions (basins) surrounding different energy minima in configuration space, in the diffusive regime (long observation times). In particular, the self-diffusivity of atoms
J. Phys. Chem. B, Vol. 112, No. 34, 2008 10629 in the system can be computed simply through the hopping process among inherent structures. Furthermore, at temperatures close to and lower than Tg, we showed that the dynamics of the system can once again be well described as a sequence of transitions between basins not only for long observation times but also for intermediate observation times, that is, when the system is considered to be trapped (cage effect regime). Finally, we concluded that the sequence of transitions between different basins at low temperatures (T < Tg), where the system is considered to be trapped, can be modeled with an equivalent first-order reaction network, where the rate constants for individual transitions are calculable from atomistic MD by use of hazard plots.15,16 This work aims to provide an investigation of the msd of atoms in a glass-forming system by expressing it rigorously in terms of subcontributions. This analysis has two objectives; one is to reveal the relative importance of interbasin transitions and local intrabasin vibrations for the dynamics, and the second is to provide the tools for lifting the dynamics from the coarsegrained inherent-structure picture to the full atomistic picture. Both objectives are strategically very important in view of the existence of efficient computational approaches, such as DIMW,10 based on multidimensional transition-state analysis and capable of overcoming the long-time limitations of MD, for tracking the evolution of the distribution of the system among its basins. By coupling DIMW with the lifting methodology introduced herein, one can construct the mean square atomic displacement as a function of time for times much longer than those accessible by brute-force MD simulation with a statistical error orders of magnitude lower than that afforded by MD or kinetic Monte Carlo. One should note that, although the inherent structure approach has been used extensively in past simulations, it has never been possible, to our knowledge, to reconstruct the full dynamics of the system on the basis of this approach. It has always been argued that the interbasin dynamics is important for long time evolution, but this idea has never been tested quantitatively. In our work, we numerically validate this idea for a number of temperatures below Tg. Our ability to reconstruct the full dynamics of the system at low temperatures provides strong evidence that the sequence of interbasin transitions can be treated as a Poisson process, which constitutes one of the main approximations invoked in practically all landscape approaches. At higher temperatures, lumping of fast-communicating basins into metabasins7 would retain the Poisson character of the dynamics; the strategy developed herein would remain valid for lifting the dynamics from the level of metabasins to the full configuration space of the system. Theory In our previous work on a model Lennard-Jones (LJ) twocomponent mixture,1 we showed that for temperatures around Tg, it was possible to follow the dynamics by constructing a network of transitions between discrete states, where each state of our model corresponds to the basin around one inherent structure in configuration space. For this system in the supercooled state, Kob,17 Shell, et al.18 have performed extensive studies, on the basis of which the mode coupling critical temperature Tc is reported as 0.435 in reduced units (∼52.2 K).17 For the same system, the glass transition temperature has been predicted18 to be Tg ) 0.32, that is, ∼38.4 K. By assuming Poisson statistics and extending hazard plot analysis15,16 to the case of multiple outcomes, we have been able to use a combination of MD and energy minimization in order to evaluate the rate constants describing transitions from
10630 J. Phys. Chem. B, Vol. 112, No. 34, 2008 one basin to the other.1 In our approach, each phase-space point along an isothermal-isochoric (NVT) MD trajectory is mapped to the local energy minimum resulting from minimization of the potential energy of the phase-space point under examination with respect to the positions of all atoms in the system. The minimization of potential energy proceeded on the basis of the Fletcher-Reeves conjugate gradient method, which was found to produce identical results to steepest descent in all cases tested. Through this mapping process, the same local minimum of the potential energy may be revisited many times within any given trajectory. We group all phase-space points according to the minimum to which they are mapped, and we assign them as belonging to the same inherent structure. In this way, from the initial MD trajectory, we obtain a reduced trajectory consisting of a sequence of inherent structures visited with recorded times of entry, exit, and residence in each inherent structure. The sequence of inherent structures visited gives us a picture of how the system evolves on its potential energy landscape, depending on the conditions. In practice, we do not perform the minimization at every MD step, simply because that would be computationally prohibitive; instead, we set an interval of several MD steps between successive minimizations. By performing the same calculation with smaller and smaller intervals, we make sure that our analysis of the dynamics is insensitive to our choice of interval. In our previous work, we showed that it is possible to follow the coarse-grained inherent dynamics and extract meaningful rate constants for interbasin transitions.1 In this work, we show that it is possible, for conditions close to and below Tg, to reconstruct the full dynamics of the system by accounting for both the interbasin transitions and the local vibrations within the inherent structures. To quantify the vibrations, we propose a procedure that traps the system within each of the basins by applying reflective conditions, that is, reversing the atomic momenta in Newton’s equations of motion every time the minimization procedure shows that a switch to a neighboring inherent structure has occurred. The inversion of the atomic momenta is followed by a randomization of the momentum of the thermostat that smoothly leads to a new trajectory. This procedure has been inspired by the novel methods developed by Voter19 for analyzing the dynamics of infrequent events. Voter has introduced the reflective conditions in order to prevent his biased system from leaving a minimum, because he is interested in biasing the local search for a transition, for example, by increasing the temperature. At the same time, he wants to avoid passing this bias to the selection of the minima to which transitions will occur; in his case, this is accomplished by reweighting the rates of different observed transitions according to an Arrhenius law, based on a harmonic approximation for the activation enthalpies. In Voter’s work,19 the reader can also find an excellent description of the importance of constraining the time evolution of the system in the case of biasing the shorttime dynamics. In this work, we do not attempt to bias the dynamics. Instead, we use reflective conditions in order to sample the local vibration within each inherent structure for the given statistical ensemble and therefore reconstruct the full dynamics of the system within a Poisson process approximation. As described above, our main objective is to start from the first-order kinetic model of transitions between inherent structures and end up recovering the full msd of the atomistic system, within the Poisson process approximation. Consequently, the first step is to start from eq 1, the so-called master equation20 that describes the time evolution of the system as a sequence
Tsalikis et al. of dynamical transitions between inherent structures under the initial condition of starting from a given state at time t ) 0.
∂Pt)0,Rft)τ,β Pt)0,Rft)τ,γ kγfβ ) -Pt)0,Rft)τ,β kβfγ + ∂τ γ*β γ*β
∑
∑
(1) where Pt)0,Rft)τ,β is the probability of observing the system in state β at time t ) τ, given that it was at state R at time t ) 0, and kβfγ is the rate constant describing the transition probability per unit of time from a state β to the state γ. In the case of systems evolving on a finite set of states obeying detailed balance (Pβ(∞) kβfγ ) Pγ(∞) kγfβ, ∀β,γ), the solution of the master equation can always be expressed via spectral decomposition.10,20 Here, Pβ(∞) is the equilibrium probability of visiting state β and is the result that will be reached at infinite time, Pβ(∞) ) limτf∞ Pt)0,Rft)τ,β, irrespectively of the initial state, for an ergodic system. In eq 1, we imply that we are interested in all solutions of the master equation that are consistent with the initial condition of starting from state R, that is, Pt)0,Rft)0,β ) δRβ, ∀β at time t ) 0. In this manner, one can describe all paths that at time t ) 0 start from state R and end up in any possible state after a time t ) τ, including R itself. The next step is to separate the inter- and intrabasin dynamics. To do so, we express the position ri(τ) of each atom i at time t ) τ as the position of atom i in the minimum of the R inherent structure R, ri,min (τ), and the deviation from it, rRi (τ) R R ri,min(τ) ≡ ∆ri (τ), as described by eq 2. R R ri(τ) ≡ riR(τ) ) ri,min (τ) + (riR(τ) - ri,min (τ)) ) min (τ) + ∆riR(τ) (2) ri,R
Given the initial (boundary) condition for the probability density, that is, the condition that the system is prepared at a time t ) 0, the master equation, eq 1, completely describes the dynamics of the system at the coarse-grained level. We now express the probability density for observing the total system as the probability of being in a given inherent structure at a given time in conjunction with the probability density of being in a given configuration within the inherent structure itself. More specifically, the probability density that describes the atomistic system at finite temperature is the product of two probability measures in two distinct spaces: the probability of being in the basin of one of the inherent structures in the space of inherent structures and the conditional probability to be in one of the configuration space points within the given basin. Note that, within the Poisson approximation, the distribution among configuration space points within the basin of an inherent structure is considered to be locally equilibrated, and therefore, the conditional probability density is considered to be dictated by the ensemble in which we work (NVT in our case). A very similar approach has been followed by Boulougouris and Theodorou in order to describe the free energy of an aging atactic polystyrene glass by an inherent structure approach, within the harmonic approximation.10 It is easy to see that the msd can be calculated from eq 3, by summing over all possible paths that start at time t ) 0 from all possible starting states R, weighted according to the probability of being initially at any R state, PR(0),∀R.
Reconstruction of the msd by Rigorous
∆r2(τ) )
〈
N
∑ (ri(τ) - ri(0))2 i)1
N
〉
J. Phys. Chem. B, Vol. 112, No. 34, 2008 10631
N
)
∑
1 〈(r (τ) - ri(0))2 〉 N i)1 i (3)
For simplicity, we will continue our derivations by focusing on a single atom i; subsequent extension to the average over N atoms is trivial. The key point to our derivation is to translate the Poisson approximation into mathematical equations for the probability density of our system. One of the main properties of the Poisson approximation is that there is no memory after a transition is made. In terms of the mathematical formulation, this requires a definition of the probability per unit time of entering any given state β at time τ - δτ and remaining in that state without exiting until time τ. We will call this probability per unit time fβ(δτ,τ). This quantity should be distinguished from Pt)0,Rft)τ,β, which only describes the probability of finding the system at a given time t ) τ in state β, irrespectively of its history. Note that, if one integrates fβ(δτ,τ) over all past times, one will recover Pt)0,Rft)τ,β:
Pt)0,Rft)τ,β )
[∑
∫0 fβ(δτ, τ) dδτ
∫0δτ [ ∑ kβfγ e- Σ k γ*β
γ*β
] } [∑
βfγt
dt )
γ*β
{ [
n
n
τ ∑ PR(0) ∫ F∆r ∑ ∫0 fβ(δτ, τ)[∫ F∆r (ri(τ) R i
R)1
β i
β)1
] }
ri(0))2 d∆riβ dδτ d∆riR
]
(6)
By substituting eq 2 and expanding, one obtains:
〈(ri(τ) - ri(0))2 〉 )
{ [
n
n
τ ∑ PR(0) ∫ F∆r ∑ ∫0 fβ(δτ, τ)[∫ F∆r ((ri,βmin - ri,Rmin)2 + R i
R)1
β i
β)1
] }
(4)
]{
Pt)0,Rft)τ-δτ,γ kγfβ + δRβ δ(τ - δτ) 1 -
γ*β
〈(ri(τ) - ri(0))2 〉 )
min min (∆riβ - ∆riR)2 + 2(ri,β - ri,R )(∆riβ -
τ
In Appendix A1, we show, for the first time to our knowledge, that fβ(δτ,τ) can be directly evaluated from Pt)0,Rft)τ,β by including the initial (boundary) condition of eq 1. There are three main contributions to fβ(δτ,τ) (eq 5). The two terms within the first bracket correspond to the probability of entering state β at time t ) τ - δτ. The last term, multiplying the bracket, corresponds to the conditional probability that once the system has entered state β at t ) τ - δτ, it does not leave this state in the remaining interval [τ - δτ,τ]. Remarkably, the probability of entering state β, as described within the square bracket of eq 5, has two contributions: the flux from neighboring states toward that state at the time t ) τ - δτ, and the probability of actually starting to count time while the system is in that particular state, a direct implication of the initial (boundary) conditions.
fβ(δτ, τ) )
On the basis of the above description, it is straightforward to realize that the msd can be evaluated as an average over all starting states R weighted according to PR(0), summed over all trajectories leading from R to any state β at time t ) τ, where the actual deviation from the minima is described by the probability density F∆riβ and the residence time δτ is integrated from 0 to τ, as formulated in eq 6.
Pt)0,Rft)τ-δτ,γ kγfβ +
]
Σ kβfγδτ δRβ δ(τ - δτ) e-γ*β (5)
In eq 5, we have expressed the probability of not leaving the state β as one minus the probability of leaving it, integrated over the time interval δτ. Note that, if one substitutes eq 5 into eq 4, one gets an integral form of the master equation that includes the initial condition of the initial differential form as a delta function; this integral form is completely equivalent to the differential form. The last quantity that we have to define is the probability, once we are in a given state β, to be distanced from the position of atom i at the minimum by ∆rβi at time τ, given that we reached state β for the first time at t ) τ - δτ. We call this probability F∆riβ d∆rβi . Note that, in the specific example of the msd, this quantity is only needed as an intermediate quantity; at the end, an artificially entrapped MD is used in order to evaluate the resulting ensemble averages. Correspondingly, we call F∆riR d∆rRi the probability of atom i being at position ∆rRi relative to its position at the energy minimum of state R at time t ) 0.
∆riR)) d∆riβ dδτ d∆riR
]
(7)
The first stage is to separate the contribution of the interbasin transitions. This only requires information on which state the system occupies at time t ) τ and is not affected by the length of time that the system has spent in that state. As has been shown,1 this is the most important contribution for the longtime diffusivity. For sufficiently long times, it will dominate the diffusive behavior of the total system.21 min At this point, we set to zero the average of the term (ri,β min β R β ri,R )(∆ri - ∆ri ), under the assumption that 〈∆ri 〉 and 〈 ∆rRi 〉 are zero. Strictly speaking, this requires a specific symmetry around the minima, which is a stronger condition than the Poisson process condition but weaker than the harmonic approximation. Because in our case we are able to describe the results by setting this term equal to zero, this is expected to be a logical assumption. In the more general case that one can think of, this approximation has to be reconsidered The next step is to distinguish two cases for the contribution arising from the (∆rβi - ∆rRi )2 term: the case where states R and β are the same and the case where R * β. By using the normalizations of F∆riR and F∆riβ, along with eq 4, we obtain:
〈(ri(τ) - ri(0))2 〉 ) n
{ [∑ ∫ n
}
∑ PR(0) ∑ Pt)0,Rft)τ,β(ri,βmin - ri,Rmin)2
R)1 n
∑
R)1
β)1 n
PR(0)
β)1
τ
0
+
fβ(δτ, τ)∫∫F∆riβ F∆riR[δRβ +
]
(1 - δRβ)](∆riβ - ∆riR)2 d∆riR d∆riβ dδτ (8) β Now, for the case of R ) β, we put ∆rβi ≡ rβi (τ) - ri,min (τ), R R R ∆ri ≡ ri (0) - ri,min(0), and we eliminate the contribution β R ri,min (τ) - ri,min (0). For the case where R * β, we expand the β term (∆ri - ∆rRi )2, and again, we set the contribution from averaging 2∆rβi ∆rRi equal to zero on the basis of the symmetry of fluctuations around the minimum.
10632 J. Phys. Chem. B, Vol. 112, No. 34, 2008
Tsalikis et al.
〈(ri(τ) - ri(0))2 〉 )
{∑
n
∑
n
n
PR(0)
R)1 n
〈(ri(τ) - ri(0))2 〉term 2 )
∑
min min 2 Pt)0,Rft)τ,β(ri,β - ri,R ) +
R)1
β)1
∑ [∫0 fβ(δτ, τ)∫F∆r τ
β i
β)1
}
n
∑
R)1
)
n
R)1
β)1
∑ Pt)0,Rft)τ,β(ri,βmin(τ) - ri,Rmin(0))2 +
∫0 fβ(δτ, τ)(〈 τ
(∆riβ)2
n
〈(ri(τ) - ri(0))2 〉 )
〉 +〈
(∆riR)2
〉 ) dδτ]
(9)
∑
n
PR(0)
R)1
n
∑ PR(0) ∑
β)1
]
[
∑ Pt)0,Rft)τ,β(ri,βmin - ri,Rmin)2+
β)1
δRβ
∫0τ [ ∑ Pt)0,Rft)τ-δτ,γ kγfβ + γ*β
〉
n
R)1
β)1
∑ PR(0) ∑ Pt)0,Rft)τ,β[(1 - δRβ)(〈(∆riβ)2〉+ 〈(∆riR)2 〉)] (10)
To summarize, there are three terms present in equation 10: (a) The interbasin contribution from the transition from one inherent structure to the other.
〈(ri(τ) - ri0))2 〉term 1 ) n
R)1
n
n
R)1
β)1
n
PR(0)
∑ Pt)0,Rft)τ,β(ri,βmin(τ) - ri,Rmin(0))2
The quantities necessary for the reconstruction of total dynamics of the system at finite temperature are: (a) The rate constants kβfγ that describe the transition between any accessible states of the system. (b) The square of the (continuous) displacement Cartesian min distance for all particles between inherent structures, (ri,β min 2 ri,R ) . (c) The msd within each basin as a function of time since entry to the basin 〈(rRi (δτ) - rRi (0))2〉. (d) The equilibrium in-basin variance of atomic positions relative to the energy minimum, 〈(∆rRi )2〉, for each basin. Several approaches may be used in order to calculate or estimate the necessary quantities. In this work, we use MD trajectories as our simulation tool and the hazard plot analysis as our strategy for extracting the rate constants from the MD trajectories. An alternative would be to use transition state theory in combination with a harmonic or quasiharmonic approximation.10,22 Such approaches have been shown to be very efficient at low temperatures, that is, below the glass transition temperature,10,19 whereas for conditions above the glass transition temperature, the current approach that is based on short MD runs is expected to be more efficient. In addition, it can be generalized very easily and will not suffer from unharmonicities in the motion at high temperatures. Molecular Simulation approach
Σ kβfγδτ δ(τ - δτ) e-γ*β 〈(riR(δτ) - riR(0))2 dδτ + n
]
〈(riR(δτ) - riR(0))2 〉 dδτ (12)
〈(∆riR)2 〉)] (13)
Finally, by substituting eq 5, assuming that 〈(∆rβi )2〉 is independent of time, and by using eq 4, we end up with our main result, eq 10. In the case where R ) β, strictly speaking, one should include two mechanisms of movement within the given basin: (a) the diffusive movement from t ) 0 to the time δτ that the system remained in the given basin and (b) the random reinsertion of the system into basin R after it has visited any other state. According to the Poisson approximation, the reinsertion is completely uncorrelated with the removal. In Appendix A2,we show how the integration over all possible points of reinsertion according to their equilibrium weight reduces the contribution to the first mechanism under the same assumptions invoked up to now.
∑
γ*β
τ
(1 - δRβ)
R)1
γ*β
∑ PR(0) ∑ Pt)0,Rft)τ,β[(1 - δRβ)(〈(∆riβ)2〉 +
∑ PR(0) ∑ [δRβ∫0 fβ(δτ, τ)〈(riR(δτ) - riR(0))2 〉 dδτ +
n
∫0τ [ ∑ Pt)0,Rft)τ-δτ,γ kγfβ +
〈(ri(τ) - ri(0))2 〉term 3 )
β)1
n
δRβ
(c) The intrabasin contribution for cases where a switch in minimum has occurred, which incorporates the variance in position due to oscillations around the original and final minima.
n
PR(0)
β)1
[
δRβδ(τ - δτ)] e
(1 - δRβ)[(∆riβ)2 + (∆riR)2 -
]
∑
- Σ kβfγδτ
F∆riR{δRβ(riβ(τ) - riR(0))2 +
2∆riβ∆riR]} d∆riR d∆riβ dδτ
n
PR(0)
(11)
β)1
(b) The intrabasin contribution for cases where no switch in minimum has occurred, which is a function of the residence time, that is, the time that the system has been spending in the minimum since its last reentry.
Our model system consists of a binary LJ mixture of particles of unequal size and interaction strength, with enhanced miscibility due to the nontraditional cross-interaction term, designed to suppress crystallization, as has been proposed by Kob et al.21,23 More information on the details of the system can be found in our previous work1 and in the original work of Kob et al. All results are reported in reduced units according to the definition used by Kob et al.21,23 In order to measure the vibrational contribution (terms 2 and 3 from equations 12 and 13) for every inherent structure, we need to have a sample of thousands of configurations in every specific basin. The sampling of a series of configurations, all belonging to the same basin (i.e., mapped to the same inherent structure), is based on the following procedure: we select one configuration that corresponds to the desired inherent structure and perform MD simulation, in the NVT ensemble, in the course of which we minimize the potential energy every 50 fs. As long as our system is still in the same basin, we store the configuration and continue our simulation. As soon as the system moves to another basin, we reverse the particle momenta and the
Reconstruction of the msd by Rigorous
J. Phys. Chem. B, Vol. 112, No. 34, 2008 10633
Figure 2. Reconstruction of the msd at T ) 23 K based on equation 10. Symbols and lines are the same as those in Figure 1. Figure 1. Reconstruction of the msd at T ) 9 K based on the firstorder kinetic model (eq 10): term 1, interbasin diffusion at the minima (eq 11); term 2, intrabasin motion (eq 12); term 3, cross term with independent intrabasin vibrations (eq 13). Both distance and time are reported in reduced units according to the definition used by Kob et al.21,23
momentum of the thermostat. In this way, we force the system not to leave the current energy well. In order not to follow the same trajectory (our trajectories are reversible for more than 500 fs in the NVT ensemble) and to sample better our energy well, after the first reversal of momenta and every 500 fs, we modify the system velocity by choosing new velocities randomly for all particles from a Gaussian distribution. The new velocity distribution is generated in such a way that the system total momentum is zero, and the new kinetic energy, via the equipartition theorem, corresponds to the imposed temperature. At the same time (i.e., every 500 fs), we modify the momentum of the thermostat, by choosing a new value randomly from a Gaussian distribution. In our previous MD simulations, we have measured the distribution of the thermostat momentum which is found to be Gaussian, and from this distribution, we choose our new random values. The described procedure has been designed to produce configurations that all belong to the basin but, at the same time, are distributed relative to each other according to the imposed equilibrium ensemble, in our case the canonical (NVT) ensemble. After we have evaluated the rate constants, either by unconstrained MD or by short entrapped MD, we have to solve the master equation. To do so, it is always beneficial if detailed balance holds between the rate constants, as would be expected from a first-order kinetic model. Because by construction our results do not have the detailed balance built in and numerical noise is expected to lead to deviation from it, we imposed the detailed balance condition in the solution of the master equation. To this end, in this work, we used a simple and effective iterative method that is described in our previous work.1 An alternative approach has been followed by Sriraman et al.24 Results and discussion In Figures 1–3, we present the results for the reconstruction of the msd based on eq 10 for our LJ mixture at temperatures of T ) 9, 23, and 38 K. The first term corresponds to the interbasin diffusion between inherent structures (energy minima), as described by eq 11. The second term describes the intrabasin motion within the minima, expressed by eq 12. The third term
Figure 3. Reconstruction of the msd at T ) 38 K based on equation 10. Symbols and lines are the same as those in Figure 1.
TABLE 1: Simulation Details on the Reported Results: Number of MD Steps Used, Number of Configurations Minimized, and Number of Inherent Structures Found T (K)
millions of steps
# configurations
# states
9.0 23.0 38.0
14.5 7.25 0.2
6000 3055 500
10 61 61
corresponds to a cross effect between uncorrelated vibrations about the original and final minima, as described by eq 13. It is evident that eq 10 is able to recover the atomistic dynamics at all three temperatures, providing a rigorous lifting of the coarsegrained model (Poisson process of successive jumps between basins) via the inherent structures and the finite temperature in-basin dynamics produced by classical MD simulation. In Table 1, we report, for the three temperatures studied (9, 23, and 38 K), the number of MD steps used to evaluate the interbasin transition rates, the corresponding time, the number of configurations that have been used for the minimization process every 2 × 10-12 s (2000 MD steps), and finally the number of inherent structures that have been identified within the given MD trajectory. From Table 1, it is clear that the number of states to which transition is attempted during the entrapment of the system in a basin, that is, during the period for which the msd remains at
10634 J. Phys. Chem. B, Vol. 112, No. 34, 2008
Figure 4. Comparison of the msd of atoms and the msd of atoms between inherent structures (energy minima), as it has been produced directly from the MD trajectory and as it has been reconstructed from the first-order kinetic model based on the hazard plot analysis of the same trajectory, at T ) 9 K.
Figure 5. Comparison of the msd of atoms and the msd of atoms between inherent structures (energy minima), as it has been produced directly from the MD trajectory and as it has been reconstructed from the first-order kinetic model based on the hazard plot analysis of the same trajectory, at T ) 23 K.
a plateau value, increases dramatically as temperature changes, as one should expect from the dramatic change in relaxation times and viscosity with temperature. In Figures 4–6, we also report a comparison of the msd of atoms between inherent structures (at the energy minima) as it has been produced directly from the MD trajectory and as it has been reconstructed form the first-order kinetic model based on the hazard plot analysis of the same trajectory (term 1, eq 11). As has been demonstrated in our previous work,1 Figures 46 indicate that our mapping to a first-order Poison process is able to recover the effective observed dynamic process of jumping from one minimum to another to an excellent degree. We also have to stress that, whereas previous studies investigated mostly the supercooled state, we have focused here on the region close to and below (T ) 38, 23, and 9 K) the glass transition temperature Tg (Tg ) 38.4 K),18 where the use of MD is clearly not sufficient to equilibrate the computer specimen.
Tsalikis et al.
Figure 6. Comparison of the msd of atoms and the msd of atoms between inherent structures (at the energy minima), as it has been produced directly from the MD trajectory and as it has been reconstructed from the first-order kinetic model based on the hazard plot analysis of the same trajectory, at T ) 38 K. The nose at the righthand side of the graph is related to the small sample size of the states sampled by the MD trajectory at this temperature.
There is a slow but well distinguishable tendency for the system to move toward the part of the landscape with lower potential energy, as has been shown in our previous work. The system is only temporarily trapped in the local landscape, and there is a slow motion toward the lower potential energies describing an aging process even above the expected Tg, because the aging is clearly dependent on the observation window used, as can be seen in Figure 5 of our previous work.1 Our observation window is set by our use of MD; any relaxation process slower than this will be either completely frozen of will appear as a slow aging process. An important message is that our approach in this temperature range should not be used to calculate equilibrium thermodynamic quantities, as one can do in the case of a supercooled liquid far above Tg. Our aim is simply to demonstrate that we can reproduce the complete dynamics under these conditions on the basis of a first-order kinetic model. This will allow the development of a dynamical sampling scheme that will be able to evolve the system to longer times in a more efficient way and therefore allow the vitrification process to be carried out under conditions much closer to experimental conditions. It should be pointed out that the analysis of the msd into inter-, intra-, and cross-basin contributions reveals that the observed plateau in the atomic msd appears at times shorter than those required for local equilibration of the system entrapped in the vicinity of inherent structures to be achieved. Sometimes, the time for attaining a plateau (caging effect) is shorter than the time for local in-basin equilibration by one, or even two, orders of magnitude. This earlier emergence of a plateau is brought about by the second and third terms of our analysis (eq 10), which, as seen in Figures 13, compensate for the slow rise of the first term. A physical explanation of why term 1 saturates much later than the total msd can be given on the basis of the fact that the plateau value of the first term corresponds to an interbasin local equilibration among basins belonging to the same metabasin. On the other hand, the main contribution to the total MSD under the examined conditions is governed by intrabasin motions (vibrations). As one would expect, the intrabasin dynamics reaches local equilibration on
Reconstruction of the msd by Rigorous a time scale shorter than that of the interbasin dynamics. This explanation of the apparent equilibration times further requires that the intrabasin vibrations be very similar for the set of basins explored, something which agrees with our observations. Summary In our attempt to investigate the role of the inherent structures in the vitrification process of glass-forming materials, we have shown1 that the dynamical transitions between basins surrounding these inherent structures can be well predicted by a firstorder kinetic scheme. In our previous work,1 a simple methodology, based on a combination of MD, potential energy minimization, and an extension of the hazard plot analysis, has been proposed in order to evaluate the rate constants that describe the transitions from each of the basins constructed around sampled inherent structures to any of the neighboring basins. The resulting network was able to reproduce the msd of atoms between inherent structures, as it has been evaluated directly from the MD trajectory. In this work, we utilize that methodology in order to completely reconstruct the total msd at finite temperature on the basis of the first-order kinetic network of interbasin transitions. The contribution to atomic displacement from in-basin vibrational motions is evaluated on the basis of small MD runs artificially trapped within each one of the explored basins. We provide the mathematical formulation for lifting the coarse-grained Poisson process model of a succession of interbasin transitions and reproducing the full dynamics of the atomistic system within the Poisson approximation. Our reconstructed results are in excellent agreement with the full atomistic MD for temperatures around and below the glass transition temperature of our LJ two-component mixture models. This clearly shows that an approach based on infrequent, uncorrelated transitions between basins is able to reproduce the full dynamics of the atomistic glassy system, where the Poisson approximation is valid. Extending the proposed approach to temperatures above Tg requires an additional step, because, as one should expect, the Poisson approximation becomes less accurate with increasing temperature. The reason can be easily understood, because the validity of the Poisson approximation relies on the rare event nature of an activated process. Nevertheless, as shown in our previous work1 and addressed by others,7,21 at higher temperaturesm the dynamics of the inherent structures (not necessarily described by Poisson statistics) is sufficient in order to describe the diffusive regime. An excellent estimate of the diffusivity is obtained from the long-time msd of atoms in the course of successive transitions of the system between single inherent structures. At high temperatures, the rare event element can be retrieved via the introduction of the metabasin picture. In this case, the system by construction spends most of its time in a set of basins that constitute a metabasin; transitions to other metabasins are rare events. In the work of Heuer et al.,7 it is evident that the use of metabasins can facilitate the description, because the transition between metabasins even at high temperatures can be described as a random walk. Applying the proposed approach above Tg has to face yet another obstacle. At those conditions, the number of inherent structures visited by the system for a given observation time grows tremendously with increasing temperature, as we have shown in this work (Table 1) and as has been addressed by others.7 Special techniques, such as interval bisectioning, may be invoked to optimize the efficiency of the energy minimization burden in constructing the inherent structure trajectories. On
J. Phys. Chem. B, Vol. 112, No. 34, 2008 10635 the other hand, it is straightforward to see that the proposed methodology of calculating rates based on the hazard plot analysis can be extended to the case where metabasins are defined, and the whole procedure can be repeated for the evaluation of rate constants between metabasins, instead of basins. Especially for temperatures above Tg, we plan to extend our methodologies by incorporating the idea of metabasins. As has been shown in the work of Heuer et al.,7 the use of metabasins allows resorting to a coarse-grained model of a random walk process on the potential energy landscape. Note that the combination of the metabasin picture with the hazard plot analysis and our proposed methodology for lifting the coarse grained dynamics to the full dynamics of the system will enable us fully to reconstruct the complete dynamics by rigorously attributing inter- and intrametabasin contributions. Finally, by combining our approach with the DIMW method,10 we will be able not only to observe the time evolution of the system but also to sample over a time range far beyond the accessible time of classical MD. Acknowledgment. This paper is part of the 03ED375 research project, implemented within the framework of the Reinforcement Programme Human Research Manpower (PENED) and cofinanced by National and Community Funds (25% from the Greek Ministry of Development-General Secretariat of Research and Technology and 75% from E.U.-European Social Fund). This work has been also funded by Scienomics SARL. Simulations have been performed under the Project HPCEUROPA (RII3-CT-2003-506079), with the support of the European Community–Research Infrastructure Action under the FP6 “Structuring the European Research Area” program. G.C.B. would also like to acknowledge the invitation to the workshop on Metastability and Rare Events in Complex Systems”, 18-22, Feb 2008, Vienna, organized by the Erwin Schro¨dinger Institute; participation in this workshop was extremely stimulating. Appendix A1 In this Appendix, we provide an integral form of the master equation that links, via eq 4, the probability Pt)0,Rft)τ,β of finding the system at a given time t ) τ in state β, irrespectively of its history, with the probability per unit time fβ(δτ,τ) of entering state β at time τ - δτ and remaining in that state without exiting until time τ. More specifically, by using eq 5 as our definition for fβ(δτ,τ), along with the master equation, eq 1, we will show that eq 4 is valid. By substituting eq 5 in the right-hand side of eq 4 and breaking the integral in two parts, we obtain:
∫0τ fβ(δτ, τ) dδτ ) ∫0τ [ ∑ Pt)0,Rft)τ-δτ,γ kγfβ + γ*β
Σ kβfγδτ δRβ δ(τ - δτ)] e-γ*β dδτ )
∫0τ [ ∑ Pt)0,Rft)τ-δτ,γ kγfβ] e- Σ k γ*β
βfγδτ
dδτ +
γ*β
∫0τ δRβ δ(τ - δτ) e- Σ k γ*β
βfγδτ
∫0τ [ ∑ Pt)0,Rft)τ-δτ,γ kγfβ] e- Σ k γ*β
βfγδτ
dδτ ) Σ kβfγτ dδτ + δRβ e-γ*β
γ*β
(A1) The next step requires the addition and subtraction of the flux out of state β within the first integral on the right-hand side of eq A1. In this way, we recreate the right-hand side of
10636 J. Phys. Chem. B, Vol. 112, No. 34, 2008 the master equation (eq 1) at time t ) τ - δτ within the integral. We then substitute eq 1 in the integral:
∫0τ fβ(δτ, τ) dδτ ) ∫0τ [ ∑ Pt)0,Rft)τ-δτ,γ kγfβ γ*β
Pt)0,Rft)τ-δτ,β
∑ kβfγ] e- Σ k γ*β
βfγδτ
dδτ +
γ*β
∫0τ [Pt)0,Rft)τ-δτ,β ∑
Σ kβfγδτ kβfγ] e-γ*β dδτ +
γ*β
Σ kβfγτ ) δRβ e-γ*β
τ ∂Pt)0,Rft)t′,β
∫0
- Σ kβfγδτ dδτ t′)τ-δτ e γ*β ∂t′ τ d - Σ kβfγδτ Pt)0,Rft)τ-δτ,β ) dδτ + (e γ*β 0 dδτ
∫
Σ kβfγτ (A2) δRβ e-γ*β
Tsalikis et al. appear explicitly in eq 10 for the msd. In this appendix, we analyze in detail the second term, that of intrabasin displacement, and we show that the two mechanisms of motion within a basin result in the expression of eq 12. In the case when the system has not left state R from the origin of time (t ) 0) up to the time of the measurement, t ) τ, there is only a single diffusion-like (ballistic for very short times) mechanism of motion within basin R. On the other hand, if within that period of time, the system has left state R and returned to it, we have to integrate and average over all events. We start by expressing the intrabasin displacement by including the intermediate point riRint (τ - δτ) to which the system returned at the last reentry time τ - δτ. Note that the history is considered to be lost at every reentry.
〈(∆riR(τ) - ∆riR(0))2 〉 ) 〈(riR(τ) - riR(0))2 〉 ) 〈(riR(τ) - riRint(τ - δτ) + riRint(τ - δτ) - riR(0))2 〉 )
At this stage, we can finalize our proof by using integration by parts:
∫0 fβ(δτ, τ) dδτ )
〈(riR(δτ) - riRint(0))2 〉 + 〈(riRint(τ - δτ) - riR(0))2 〉 +
τ
[∫ ( τ
0
)
∂Pt)0,Rft)τ-δτ,β - Σ kβfγδτ e γ*β dδτ ∂δτ
|
Σ kβfγδτ [Pt)0,Rft)τ-δτ,β e-γ*β ]
∫0τ
δτ)τ
+
δτ)0
∂Pt)0,Rft)τ-δτ,β - Σ kβfγδτ Σ kβfγτ dδτ + δRβ e-γ*β ) e γ*β ∂δτ
Σ kβfγτ Σ kβfγ0 + Pt)0,Rft)τ-0,β e-γ*β + -Pt)0,Rft)τ-τ,β e-γ*β Σ kβfγτ Σ kβfγτ ) -δRβ e-γ*β + δRβ e-γ*β Σ kβfγ0 Σ kβfγτ + δRβ e-γ*β ) Pt)0,Rft)τ-0,β e-γ*β
Pt)0,Rft)τ,β (A3) which is the left-hand side of eq 4. Hence, we have proven that eq 4 follows from eqs 1 and 5. Eqs 4 and 5 constitute an integral formulation of the master equation. Recently, during the workshop on Metastability and Rare Events in Complex Systems (Vienna 18-22 Feb 2008), a derivation of the milestoning methodology of Faradjian and Elber came to our attention. Milestoning is a novel methodology for sampling rare events without the need of a reaction coordinate by approximating the reaction as a one-dimensional jump process, not necessarily Poisson, between nonintersecting hyperplanes leading from reactants to products. The probability of observing the system in each of the hyperplanes is then evaluated via an integro-differential equation formulated with similar arguments.25,26
2〈(riR(τ) - riRint(τ - δτ))(riRint(τ - δτ) - riR(0))〉 (A4) where we have invoked a time shift in the first term, 〈(rRi (δτ) - riRint (0))2〉 ≡ 〈(rRi (τ) - riRint (τ - δτ))2〉, because between the intermediate point of last reentry and the ending point, only the relative difference in time is of importance. Our analysis requires that we integrate over all possible starts riRint (0), all possible reentry points riRint (τ - δτ), and finally over all possible ending points. Each point has to be weighted according to the probability of observing it; the starting points and reentry points are sampled on the basis of the assumption of local equilibrium within the basin, whereas the ending points result from deterministic trajectories averaged over all possible paths with the given length in time. Therefore, we define the following probabilities, FriR driRstart, FriR driRint and FriR , riR , δτ driRend where the subscripts start int end int indicate the dependence of each of the probabilities. Note that, for clarity reasons, the averaging over all atoms is not included. Nevertheless, it is trivial to see that expressions developed below are valid for every atom i and therefore for the average, too.
〈〈〈(∆riR(τ) - ∆riR(0))2 〉end 〉int 〉start )
∫r
R istart
FriR
start
∫r
R iint
FriR
int
∫r
R iend
,rR ,δτ end iint
(riRint(τ - δτ) - riRstart(0))2 + 2(riRend(τ) - riRint(τ - δτ))(riRint(τ δτ) - riRstart(0))] driRend driRint driRstart (A5)
Appendix A2 The main goal of this work is to express the msd of the atomistic MD in terms of the inter- and intrabasin motion of atoms in a glass-forming system close to and below the glass transition temperature. To do so, we have decomposed the total atomic msd into three contributions: the interbasin contribution, reflecting switches in the inherent (minimum energy) structure; the intrabasin displacement, applicable when the system at the considered observation time finds itself in the same basin as that from which it started; and a contribution from the uncorrelated intrabasin vibrations in starting and final basins, applicable when a basin switch has occurred. These contributions
[(riRend(δτ) - riRint(0))2 +
FriR
It is easy to see that the first two terms within the brackets parenthesis correspond to the ensemble averages over the process of diffusing for time δτ after the reentry point and to the average of reentry points over all starts. 〈〈〈(∆riR(τ) - ∆riR(0))2 〉end 〉int 〉start ) 〈(riRend(δτ) - riRint(0))2 〉 + 〈(riRint(τ - δτ) - riRstart(0))2 〉 + 2 riRstart(0))[
∫
riRend
FriR
∫
,rR ,δτ end iint
riRstart
FriR
start
∫
riRint
FriR (riRint(τ - δτ) int
(riRend(τ) - riRint(τ - δτ)) driRend] driRint driRstart )
Reconstruction of the msd by Rigorous
J. Phys. Chem. B, Vol. 112, No. 34, 2008 10637
)〈(riRend(δτ) - riRint(0))2 〉 + 〈(riRint(τ - δτ) - riRstart(0))2 〉 + 2
∫
riRstart
FriR
start
2
∫
riRint
∫
riRstart
FriR riRint(τ - δτ)[ int
FriR
start
∫
riRint
∫
riRend
FriR
,rR ,δτ end iint
FriR riRint(τ - δτ)[ int
∫
riRend
riRend(τ) driRend] driRint driRstart -
FriR
,rR ,δτ end iint
riRint(τ -
δτ) driRend] driRint driRstart 2
∫
FriR
∫
FriR riRstart(0)[
∫
FriR
∫
FriR riRstart(0)[
2
riRstart
riRstart
start
start
riRint
riRint
int
int
∫
FriR
riRend(τ) driRend] driRint driRstart +
∫
FriR
riRint(τ - δτ) driRend] driRint driRstart )
riRend
riRend
,rR ,δτ end iint
,rR ,δτ end iint
)[〈(riRend(δτ) - riRint(0))2 〉 + 〈(riRint(τ - δτ) - riRstart(0))2 〉 + 2〈riRint(τ - δτ)〈(riRend(τ))〉end 〉int - 2〈(riRint(τ - δτ))2 〉int 2〈riRstart(0)〈〈riRend(τ)〉end 〉int 〉start + 2〈riRstart(0)〈riRint(τ - δτ)〉int 〉start] (A6)
By expanding and reordering the averages we obtain:
In the manipulations above, we have added and subtracted the position at the minimum riRmin and assumed that the average vector difference from both starting and ending positions to the position of the minimum is zero, invoking a local symmetry. In a few tests that we have performed on our MD trajectories, we have seen that for some of the minima, this difference is indeed close to zero and can be neglected in comparison to the rest of the terms, but it is difficult to extrapolate this notion to any system, and a case by case examination may be required. Finally, we have assumed that 〈(riRint(τ - δτ) - riRmin)2〉 ) 〈(riRstart(0) - riRmin)2〉 as
the result of the selection according to the equilibrium probability distribution of both the starting and the reentry points. We have finally shown that the contribution we have to account for in connection with the intrabasin displacement term comes from the diffusive motion since the last reentry of the system into state R, up to the time of the measurement. The only additional piece of information needed is the probability density for entering state R at time τ- δτ and not having left at τ. This probability density is given by eq 5. References and Notes (1) Tsalikis, D.; Lempesis, N.; Boulougouris, G. C.; Theodorou, D. N. J. Phys. Chem. B 2008, 112, 10619. (2) Debenedetti, P. G.; Stillinger, F. H. Nature 2001, 410, 6825. (3) Go¨tze, W.; Sjo¨gren, L. Rep. Prog. Phys. 1992, 55, 241. (4) Go¨tze, W. In Liquids, Freezing and the Glass Transition, Les Houches Session LI, 1989; Hansen, J.-P., Levesque, D., Zinn-Justin, J., Eds.; North-Holland: Amsterdam 1991; pp 287-499. (5) Shell, S. M.; Debenedetti, P. G.; Panagiotopoulos, A. Z. Fluid Phase Equilib. 2006, 241, 147. (6) Goldstein, M. J. Chem. Phys. 1969, 51, 3728. (7) (a) Bu¨chner, S.; Heuer, A. Phys. ReV. L. 2000, 84, 2168. (b) Doliwa, B.; Heuer, A. Phys. ReV. E 2003, 67, 030501. (c) Saksaengwijit, A.; Doliwa, B.; Heuer, A. J. Phys.:Condens. Matter 2003, 15, S1237. (d) Doliwa, B.; Heuer, A. Phys. ReV. Lett. 2003, 91, 235501. (8) (a) Wang, C.; Stratt, M. R. J. Chem. Phys. 2007, 127, 224503. (b) Wang, C.; Stratt, M. R. J. Chem. Phys. 2007, 127, 224504. (9) Sciortino, F. J. Stat. Mech. 2005, P050515. (10) Boulougouris, G. C.; Theodorou, D. N. J. Chem. Phys. 2007, 127, 084903. (11) Jain, T. S.; de Pablo, J. J. Phys. ReV. Lett. 2004, 92, 155505. (12) Stillinger, F. H.; Weber, T. A. Phys. ReV. A 1982, 25, 978. (13) Henkelman, G.; Jo´nsson, H. J. Chem. Phys. 1999, 111, 7010. (14) Boulougouris, G. C.; Frenkel, D. J. Chem. Theory Comput. 2005, 1, 389. (15) (a) Helfand, E. J. Chem. Phys. 1978, 69, 1010–1018. (b) Helfand, E.; Wasserman, Z. R.; Weber, T. A. Macromolecules. 1980, 13, 526. (16) Helfand, E. J. Chem. Phys. 1978, 69, 1010. (17) Kob, W. Phys. ReV. E 1995, 51 (5), 4626. (18) Shell, S. M.; Debenedetti, P, G.; Panagiotopoulos, A. Z. Fluid Phase Equilib. 2006, 241, 147. (19) (a) Voter, A. F. Phys. ReV. Lett. 1997, 78, 3908. (b) Voter, A. F.; Montalenti, F.; Germann, T. C. Annu. ReV. Mater. Res. 2002, 32, 321. (20) Reichl, L. E. A Modern Course in Statistical Physics, 2nd ed.; Wiley: New York, 1998; Chap. 5, pp 229-261. (21) Shell, S. M.; Debenedetti, P. G.; Stillinger, F. H. J. Phys. Chem. B 2004, 108, 6772. (22) Kopsias, N. P.; Theodorou, D. N. J. Chem. Phys. 1998, 109, 8573. (23) Kob, W. J. Phys.: Condens. Matter 1999, 11, R85. (24) Sriraman, S.; Kevrekidis, I. G.; Hummer, G. J. Phys. Chem. B 2005, 109, 6479. (25) Faradjian, K. A.; Elber, R. J. Chem. Phys. 2004, 120, 10880. (26) West, A. M.; Elber, R.; Shalloway, D. J. Chem. Phys. 2007, 126, 145104.
JP8013223
