On the Role of Inherent Structures in Glass-Forming Materials: I. The

J. Phys. Chem. B 2008, 112, 10619–10627

10619

On the Role of Inherent Structures in Glass-Forming Materials: I. The Vitrification Process Dimitrios G. 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 Eduard VII, 75009 Paris, France ReceiVed: February 12, 2008; ReVised Manuscript ReceiVed: April 21, 2008

In this work, we investigate the role of inherent structures in the vitrification process of glass-forming materials by using a two-component Lennard-Jones mixture. We start with a simplified model that describes the dynamics of the atomistic system as a Poisson process consisting of a series of transitions from one potential energy minimum (inherent structure) to another, the rate of individual transitions being described by a first-order kinetic law. We investigate the validity of this model by comparing the mean square displacement resulting from atomistic molecular dynamics (MD) trajectories with the corresponding mean square displacement based on inherent structures. Furthermore, in the case of vitrification via stepwise cooling, we identify the role of the potential energy landscape in determining the properties of the resulting glass. Interestingly, the cooling rate is not sufficient to define the resulting glass in a stepwise cooling process, because the time spent by the system at different temperatures (length of the steps) has a highly nonlinear impact on the properties of the resulting glass. In contrast to previous investigations of supercooled liquids, we focus on a range of temperatures close to and below the glass transition temperature, where the use of MD is incapable of producing equilibrated samples of the metastable supercooled state. Our aim is to develop a methodology that enables mapping the dynamics under these conditions to a coarse-grained first-order kinetic model based on the Poisson process approximation. This model can be used in order to extend our dynamical sampling ability to much broader time scales and therefore allow us to create computer glasses with cooling rates closer to those used experimentally. In a continuation to this work,1 we provide the mathematical formulation for lifting the coarsegrained Poisson process model and reproducing the full dynamics of the atomistic system. Introduction The role of glasses in our everyday life is far more important and essential than one would intuitively expect. This happens because their use is not limited to the traditional inorganic glasses, such as wine bottles and window panes, but includes also the technologically most important category of organic (usually polymer) glasses. Today, it is known that any substance is able to form a glass, under the condition that the nucleation of a crystal phase is repressed or limited. The correct question is consequently not whether a substance is able to vitrify but under which conditions this can happen.2 Several attempts have been made to categorize glasses, according to various criteria. A first categorization of glassforming liquids was attempted by Angell.3 The distinction is based on the temperature dependence of dynamic viscosity. A linear change of the logarithm of the relaxation time with reciprocal temperature (usually reduced by Tg) corresponds to an Arrhenius-type dependence. Liquids for which the deviation from the Arrhenius-type behavior is small are called strong glass formers, whereas, on the contrary, those which present large deviations are called fragile glass formers. Another distinction of glasses was proposed by Sciortino,4 who distinguished glasses into attractive and repulsive. The molecular origin of attractive glasses is the existence of strong (in relation to the temperature) attractive forces. In cases where the range of attractive forces * Author to whom correspondence should be addressed. Fax: +30 210 772 3112. E-mail: [email protected]. † National Technical University of Athens. ‡ Scienomics.

is much smaller than the particle size, the glassy state can be found at very low densities.5 On the other hand, repulsive glasses are the result of hard-core interatomic interactions due to the short-range repulsion of the electron clouds and are usually observed at high densities, where repulsive interactions become dominant. Various theories have been proposed for the interpretation of glass formation. One such theory is the mode-coupling theory (MCT).6 MCT attempts to describe the dynamic behavior of strongly supercooled liquids at temperatures above the glass transition temperature Tg. In its simplest version, the so-called ideal MCT, the model predicts the existence of a temperature Tc at which the system goes from ergodic to nonergodic behavior. For temperatures lower than the melting temperature Tm, thermodynamics dictates that a crystallizable material should be in the crystalline solid state; nevertheless, because of kinetic reasons, it may be trapped in the supercooled liquid state for time scales of the range of the experimental observation. In the case of semicrystalline materials, part of the system is crystallized, whereas the rest is kinetically trapped. A lower bound of the supercooled liquid state is the glassy state, where the system has been vitrified. Therefore, the supercooled liquid state describes the condition in which a liquid lies between crystallization and vitrification. For a crystallizable material, the supercooled liquid is metastable with respect to the stable crystalline state. As the supercooled liquid is cooled to increasingly lower temperatures, its density and viscosity increase, and the molecules move more and more slowly. This implies that the molecules, as temperature decreases, need an increasingly

10.1021/jp801296k CCC: $40.75  2008 American Chemical Society Published on Web 08/01/2008

10620 J. Phys. Chem. B, Vol. 112, No. 34, 2008 longer adjustment timesor relaxation timesin the new state. A comprehensive review of supercooled liquids and their glasses can be found elsewhere.7 Extensive studies of atomic systems via molecular simulation resulted in a cage model as the underlying mechanism responsible for the dynamical entrapment of repulsive soft glass-forming materials.8 Every particle is trapped in the cage formed by the neighboring particles that surround it, and a long time is needed for the particle to escape from its cage. Note that the particles forming the cage are, of course, trapped themselves in cages as well, and thus, the motion of all particles is slowed down. With decreasing temperature, the cages become more and more rigid, and thus, the time needed to break them up increases. The above-mentioned MCT9 is an attempt to describe this phenomenon in a self-consistent way and hence to rationalize the dramatic increase of the relaxation time with decreasing temperature on a microscopic level.8 On the other hand, there are theories which focus on the topography of the multidimensional energy hypersurface of the system. These theories try to connect the phenomenon of glass transition with the thermodynamic arrest of the system in the vicinity of some local minima in its energy landscape.10–12 Likewise, they try to connect properties of glasses, such as the faster β-relaxation time and the slower R-relaxation time, with motions within and transitions between the various energy minima. The most important features of the potential energy landscape are local minima or inherent structures,13 where the force vanishes. Associated with an inherent structure is the idea of a basin around the potential energy minimum, formed by the locus of configurational points for which a steepest descent minimization will result in the same local minimum, as it has been introduced by Stillinger and Weber.13 The importance of the inherent structure approach has been linked with the growing use of computer simulation in material science. Since the early work of Theodorou and Suter,14 the inherent structure approach has been used to predict structure and elastic constants in polymer glasses by molecular simulation. Another imporant step in the description of the potential energy landscape is the introduction of the metabasins:15 each metabasin consists of a number of adjacent basins, separated by barriers that are low in relation to the thermal energy kBT. Although there are several approaches for assigning the basins to metabasins, the basic idea is that the transitions between basins within the same metabasin should be much faster than the transitions between basins of different metabasins. Interestingly, Doliwa and Heuer have shown that the dynamics of a system above the mode-coupling temperature can then be interpreted as a random walk between metabasins.16 This work aims at investigating the significance of the inherent structure model for the prediction of the dynamical behavior of glass-forming materials. In order to do so, the atomic mean square displacement (msd) as a function of time is used as a property to describe the dynamics of the system. We find that, for temperatures significantly higher than Tg, the dynamics of the system can be well described by the process of hopping between different energy minima in the diffusive regime (long time of observation). In particular, the self-diffusivity of atoms can be captured well simply by the hopping process among the inherent structures. Furthermore, at temperatures close to and lower than Tg, we show that the dynamics of the system can once again be well described by the transitions between the energy minima not only for long observation time but also for intermediate observation time, that is, when the system is considered to be trapped (cage-effect regime). The last conclu-

Tsalikis et al. sion is further confirmed by an analysis of the msd into various contributions in a continuation of this work.1 Finally, we conclude that the transitions between inherent structures at low temperatures, where the system is considered to be trapped, can be modeled with an equivalent network of first-order reactions, where the reaction (elementary transition) rate constants can be obtained by infrequent event analysis through the use of hazard plots.17 Theory Within the classical mechanics approximation, molecular dynamics (MD) trajectories are used to simulate the dynamical behavior of a system. Each realization of the system is depicted as a point of the trajectory in phase space, that is, the space spanned by the positions and momenta of all particles. The phase-space trajectory gives a complete description of the dynamics, provided it is sufficiently long. For Hamiltonian systems, the trajectory can never cross itself at any time. On the other hand, the projection from phase space to configuration space, that is, the space spanned by the positions of all particles in the system, is not one-to-one, and by all means, the system may revisit several times the same configuration (point in configuration space). In practice, the generation of atomistic MD trajectories is a computationally very intensive process, the times it can follow not exceeding a few microseconds on today’s fastest machines. A coarse-grained description of the dynamics, cast in terms of a smaller number of variables and capable of addressing longer time scales, is highly desirable. Probably, the most rigorous approach for developing such a coarse-grained description is the Mori-Zwanzig projection operator formalism.18,19 Provided that one is able to separate the observables into slow and fast (in this case, one works in the Hilbert space spanned by the observables), this formalism makes possible a rigorous derivation of the equation of motion in the projected space of the slow degrees of freedom, whereas the fast observables provide a bath of random thermal forces and frictional forces. The resulting projected equation of motion is solvable, provided that there is a significant decoupling between the slow and the fast degrees of freedom, where the decoupling is usually based on a separation of time scales. Unfortunately, to our knowledge, there is no rigorous way to identify a priori what is slow and what is fast, and the separation always relies on intuition based on knowledge of the physical system. In this work, we attempt the opposite approach. We start with a simplified model that describes the dynamics of the atomistic system as a Poisson process consisting of elementary transitions with first-order kinetics from the vicinity of one potential energy minimum (inherent structure) to that of another, and we investigate the validity of this model by comparing the msd resulting from the original atomistic trajectory with the msd based on the inherent structure model. To evaluate rate constants for the elementary transitions between inherent structures, we invoke a hazard plot analysis17,20 of the atomistic MD trajectory. Hazard plot analysis is based on the evaluation of the cumulative hazard. The hazard rate, h(t), is defined such that h(t)dt is the probability that a system that has survived a time t since its last transition will undergo a transition at the time between t and dt. We define the cumulative hazard as

H(t) )

∫0t h(t ′ ) dt′

The probability that a transition occurs in a time less than t since the last transition is P(t) ) 1 - exp[-H(t)]. The

Vitrification Process probability density for such a transition to occur is dP/dt ) h(t) exp[-H(t)]. For a Poisson process with rate constant λ, the cumulative hazard is H(t) ) λt, and the probability is P(t) ) 1 - exp[-λt]. In our case, the Poisson rate λ can be extracted as the slope of a plot of the cumulative hazard versus the residence time. In the literature, there are several approaches to the estimation of rate constants that are far more rigorous and exact compared to the proposed approach.21–23 Nevertheless, we will show that the hazard-plot-based methodology for evaluation of the rate constants, proposed in this work, is sufficient to reproduce the dynamics of the inherent structures. We target on a variety of conditions, from high temperature, where traditional works on the simulation of supercooled liquids have been conducted, down to low temperatures. We are not interested in the range where MD is able to provide a locally equilibrated sample of the metastable supercooled liquid but rather focus on conditions close to and below Tg, where the use of MD results in an inefficient sampling of the energy landscape. Our interest stems from the importance of this region in shaping the properties of the resulting glass, commonly manifested as a cooling rate effect. For example, it is known that the glass transition temperature is a function of the cooling rate, and the faster the cooling, the higher Tg will be. Our approach indicates that this effect is created by the limited time that the system has at its disposal to explore its energy landscape. The actual part of the landscape where the system will be trapped as the temperature is lowered below the glass transition temperature depends on the time that the system was allowed to relax. From the temperature dependence of the relaxation time (or, equivalently, of viscosity) in glass-forming materials, we know that the time needed to relax rises extremely steeply close to and below Tg. It is easy to imagine that from the temperature range where MD is able to sample the landscape efficiently down to the point where the relaxation time reaches the order of seconds, there will be a region where MD is not able to produce anything but a small window of the approach to the local equilibrium state of the corresponding supercooled liquid. By providing a way of mapping the dynamics of such systems onto the coarse-grained level of inherent structures, we aim at developing a method that will allow the observation of a much broader window of time scales, allowing computer vitrification under much slower cooling rates, such as those used experimentally. In our approach, each trajectory phase-space point of MD is mapped to the minimum resulting from the minimization of the potential energy with respect to the positions of all atoms of the system at constant volume. As a consequence, the same minimum of the potential energy is revisited (mapped) many times within any given trajectory. We group all phase-space points according to the mapped minima, and we assign them as belonging to the same inherent structure. In this way, we obtain a view of the distribution among basins in the potential energy landscape, and we examine how this distribution changes, depending on the adopted conditions. Molecular Simulation Approach In the present work, a two-component Lennard-Jones (LJ) mixture was studied. The total number of atoms was set to 641. The selection of the specific system is based on its simplicity and on the existence of extensive prior research results on it. One of its main advantages is the suppression of the tendency for crystallization present in pure LJ systems. The mixture, initially proposed by Kob et al.,24 consists of two different types of atoms, A and B. The mixture has a concentration of 80% in

J. Phys. Chem. B, Vol. 112, No. 34, 2008 10621 A atoms and 20% in B atoms. The parameters of the model have been selected24 in such a way that, although the two species have different sizes and interaction strengths, demixing is suppressed in order to suppress nucleation. Although A atoms are larger than B atoms, they are assumed to have the same mass mA ) mB ) 6.634 × 10-26 kg. The LJ interaction parameters are εAA ) 1.65678 × 10-21 J, σAA ) 3.4 × 10-10 m, εBB ) 0.82839 × 10-21 J, σBB ) 2.992 × 10-10 m, εAB ) 2.48517 × 10-21 J, and σAB ) 2.72 × 10-10 m. The unit for reducing time is selected24 as [mAσ2AA/(48εAA)]1/2 ) 3.10 × 10-13 s, and the unit for temperature is εAA/kB ) 120 K. If the above LJ interaction parameters are reduced25 by the values of the A-A interaction parameters, they read24 εAA ) 1.0, σAA ) 1.0, εBB ) 0.5, σBB ) 0.88, εAB ) 1.5, and σAB ) 0.8. For this system in the supercooled state, Kob26 and Shell et al.27 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).26 For the same system, the glass transition temperature has been predicted27 to be Tg ) 0.32, that is, ∼38.4 K. In order to have continuous first and second derivatives, necessary for MD but even more vital for the energy minimization process, the LJ intermolecular potential has been truncated at 1.71σAA, and a fifth-order spline has been added up to rspl ) 2.0σAA, bringing potential and first and second derivatives to zero.14 The system has been studied for a series of temperatures under constant density 1.1908σAA-3 in the canonical (NVT) ensemble. In order to validate our code, we have performed extended comparisons between our results and those originally reported by Kob et al.24,26 for both thermodynamic and dynamical properties. Data from these comparisons are not reported here,but showed excellent agreement. The standard Nose´-Hover thermostat, in combination with the Verlet algorithm as implemented by Melchionna et al.,28 has been used for the solution of Newton’s equations of motion in the NVT ensemble. The time step has been set to 1 fs, or 0.32 × 10-2 in reduced units. In the following, we summarize the procedure adopted for the evaluation of the rate constants for the transitions from one inherent structure to a neighboring one. Along our atomistic MD simulations, at regular time intervals, the potential energy was minimized with the method of conjugate gradients.29 Then, in a separate step, the coordinates of the center of mass for each configuration were subtracted from all atomic coordinates of the energy-minimized configuration (inherent structure). When calculating the center of mass, we used the images of atoms originally in the simulation box evolved along continuous paths via the MD; that is, we did not allow image switches via the periodic boundary conditions (cf. calculation of self-diffusivities via MD). This was done to facilitate comparisons of inherent structures with each other and also with the MD configurations, the center of mass of which remained motionless during the integration of the equations of motion. A hazard plot17 analysis was conducted to extract the transition rate constants based on the residence times between successive transition events. In our case, the residence time is the time interval during which the system stays in the neighborhood of a certain minimum, that is, the time interval for which the minimization procedure along the trajectory leads to the same minimum; this interval terminates as soon as a switch in minimum occurs. For any given minimum β, we have to count (a) the number n of times that we leave this minimum to any other minima and (b) the time elapsed from entry into β until

10622 J. Phys. Chem. B, Vol. 112, No. 34, 2008 exit to any minimum neighboring γ, for each of the n visits to minimum β. Hazard analysis requires that we order the n transitions according to the value of the residence time between transitions that preceded them: t1 e t2 e t3 e tk e... e tn. By ordering the transition times, the cumulative probability of observing a transition up to time tk is approximated in the hazard analysis as17 k/(n + 1). The transition rate out of minimum β ˆ (k) versus tk, where H ˆ (k) ) will equal the slope λ of a plot of H ∑l)k-1 l)0 1/(n - l) is an estimate of the cumulative hazard function, n is the total number of transitions, k is the current transition, and tk is the residence time of this transition. A linear ˆ (k) on the residence time implies that the dependence of H process is of Poisson’s type. In a Poisson process, each step of the process does not have memory of the system’s previous condition; in our case, the residence time in a minimum is assumed independent of the system’s previous history; that is, transitions are uncorrelated events. In a Poisson process, with the assumption that the cumulative probability is of the form

P(t) )

{

1 - e-λ(t-t0) 0

t > t0 t e t0

}

the gradient λ will be

^ λ)

n - rl n

∑

k)rl + l

(tk - tl)

with n being the total number of transitions and tk the residence time of transition k. The quantity rl stands for the number beyond which transitions are taken into account for the calculation of the slope. It is designed in order to exclude some (rl) initial transitions at short times, where the Poisson approximation is not expected to hold. The time tl is the total time for this first rl cutoff. If no cutoff is used, one recovers that the total rate constant (sum of rate constants to other minima) leading out of minimum β is the inverse mean first-passage time for the event of leaving the minimum to reach any other minimum. We should note that the original hazard analysis has been developed for a single event, whereas in this work, we extend its application to the case of multiple outcomes. To do so, we identify as a measured event (as described by the original scheme of hazard analysis),17 any transition to any neighboring minimum. That is to say, an event is any transition to a new minimum. In this way, the resulting rate corresponds to the sum of rates out of the given minimum. In order to recover the rate constants to each of the neighboring minima, we also evaluate the conditional probability of leaving minimum β and ending up in each of the neighboring minima γ, ωβfγ. The rate constant for the transition from minimum β to any minimum γ is equal to the product of ωβfγ times the total rate constant for leaving minimum β: kβfγ ) ωβfγ∑γ′kβfγ′ (irrespectively of the time this took). This is the exact inverse procedure of what is used in an event-driven kinetic Monte Carlo simulation,29 where one first evaluates the time that any event will take place on the basis of the sum of the rate constants and then selects one of the events according to their individual rate constants. In Figure 1, we present a characteristic hazard plot for an inherent structure studied in this work at T ) 9 K. As described above, our methodology performs minimization steps at fixed intervals. In Figure 1, we report the hazard plots for two calculations: one where the time between minimizations was set to 0.8 ps, and one where it was set to 2 ps. As reported in the work of Helfand,17 a linear dependence of the cumulative hazard on the transition time indicates that the

Tsalikis et al.

Figure 1. Hazard plot of the time between transitions to a neighboring basin for one inherent structure at T ) 9,K. Two sets of calculations are presented, with a 0.8 ps interval between minimizations (squares) and a 2 ps interval between minimizations (crosses).

assumption of the Poisson process is valid. The slope of the cumulative hazard corresponds to the transition rate constant. By applying the same approach to each of the inherent structures which resulted from the minimization procedure of the MD trajectory, a network of first-order transitions is created. The time evolution of the probability of observing the system in inherent structure β at time τ results from the solution of the master equation

∂Pt)0,Rft)τ,β ) -Pt)0,Rft)τ,β ∂τ

∑ kβfγ + ∑ Pt)0,Rft)τ,γkγfβ γ

γ

(1) Equation 1 actually describes a set of linear first-order differential equations, where the first term on the right-hand side corresponds to the flow of probability from state β to all possible other states γ and the second term corresponds to the flow from other states γ to state β. State R determines the initial boundary condition for a time equal to zero. The integration of eq 1 is always feasible, provided that detailed balance holds, with the use of a spectral method, which leads to an analytic expression of the evolution in time, based on the eigenvectors and the eigenvalues of the transition matrix K.30 The elements of the transition matrix K are defined as Kβγ ) kγfβ ∀γ*β, Kββ ) -∑γkβfγ. Here, ∑γkβfγ are the rate constants calculated from the slopes of the hazard plots, and the individual rate constants kβfγ are evaluated with the additional knowledge of the conditional probabilities ωβfγ for starting from minimum β and ending up in minimum γ, as described above. To solve the master equation, eq 1, 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 detailed balance built-in and numerical noise is expected to lead to deviations from it, we imposed the detailed balance condition in the solution of the master equation. An alternative approach has been followed by Sriraman et al.;31 distinct parts of phase space were identified, and effective rate constants were calculated by connecting those parts, by imposing the condition of detailed balance. In this work, we used a simple and effective self-consistent method that modifies the transition matrix so as to obey detailed balance, without changing the equilibrium probabilities that result directly from the balance condition, eq 2 below. The equilibrium probabilities Pβ,eq are calculated directly from (0) the raw estimates of the rate constants kβfγ using eq 2.

Vitrification Process

J. Phys. Chem. B, Vol. 112, No. 34, 2008 10623

0 0 0 0 ) ∑ Pγ,eqkγfβ or Pβ,eq∑ kβfγ ) ∑ Pγ,eqkγfβ ∑ Pβ,eqkβfγ ( )

γ

( )

( )

γ

γ

( )

γ

(2) To solve eq 2, an iterative scheme is employed. A converged value forPβ,eq is obtained by an iterative solution of eq 3, by 0 starting from any initial guess for the probabilities, Pβ,eq , under the constraint that they are normalized to 1:

(l+1) Pβ,eq )

(l) (0) kγfβ ∑ Pγ,eq γ

(3)

(0) ∑ kβfγ γ

If the rate constants obey detailed balance, the transition matrix K can be transformed to a symmetric one, Π, with the same eigenvalues. The diagonal elements of the new, symmetrized matrix are (Pβ/Pγ)1/2kβfγ ) (Pγ/Pβ)1/2kγfβ. In our case, we build this symmetric matrix as the average of the two estimates from the forward and the backward transition, Πβγ ) [(Pβ,eq/ (0) (0) Pγ,eq)1/2kβfγ +(Pγ,eq/Pβ,eq)1/2kγfβ ]/2, by using the raw estimates (0) of the rate constants, kβfγ. It is now trivial to invert the transformation from Π to K and evaluate a set of rate constants (0) kβfγ that result in the same Pβ,eq as the kβfγ but also obey detailed balance exactly, by using kβfγ ) Πβγ(Pγ,eq/Pβ,eq)1/2.

Figure 2. Comparison of msd of A atoms between MD and inherent structure trajectories for T ) 67 K. The x-axis is reduced time, and the y-axis describes the msd, also in reduced units.

Results and Discussion In Table 1, we present the temperatures, in Kelvin and in reduced units, at which we studied sequences of transitions between minima of the potential energy, along with the length of each run in time steps and in absolute time (ns). As described above, for each of the temperatures reported in Table 1, minimizations of the potential energy have been performed at fixed time intervals. In this way, each point along the MD trajectory produced a minimum energy point in the potential energy landscape. After subtraction of the center of mass of the continuously displaced initially selected images of atoms in each configuration, the msd of atoms in the inherent structures is calculated. The results are compared with msd values of the MD for the same temperatures (Figures 2–4). This msd comparison appears in the following figures for the three temperatures studied. From Figures 2–4, one can draw the following conclusions: (i) Initially, the msds along inherent structure trajectories are significantly lower than those of MD. This is evident at the two lowest temperatures but also at 67 K up to t* ) 600. This observation is justified. Local vibrational motions around the minima have been removed in the inherent structure trajectories. If the system remains trapped in the vicinity of a few minima, such vibrational motions constitute the predominant contribution to displacement; therefore, displacement is severely underestimated by the inherent structure trajectories. The situation is dramatically different at high temperatures and long times (T ) 67 K, t* > 600), where the system executes a large number of transitions between different minima. Minima which differ TABLE 1: Temperature Values for Which Potential Energy Minimizations Have Been Performed, Along With the Length of Each Run in Time Steps and in ns. NVT-runs T (K) 1 2 3

67.0 38.0 9.0

T* 0.558 0.317 0.075

millions time timestep minimization of steps (ns) δt (fs) every (ps) 14.5 14.5 14.5

14.5 14.5 14.5

1.0 1.0 1.0

2.0 2.0 2.0

Figure 3. Comparison of msd of A atoms between MD and inherent structure trajectories for T ) 38 K. Axes and symbols are the same as those in Figure 2.

Figure 4. Comparison of msd of A atoms between MD and inherent structure trajectories for T ) 9 K. Axes and symbols are the same as those in Figure 2.

a lot from each other are visited, resulting in large changes in minimum-to-minimum displacements. The displacement due to local vibrations becomes insignificant relative to the diffusive displacement between basins, and the inherent structure trajectory provides an excellent approximation to the MD trajectory. (ii) A second observation concerns the gradient of the two curves at each temperature. It is generally expected that the gradient of the two curves should be similar at long times. The long-time functional form at 67 K in the log-log plot of the MD in Figure 2 is linear with a slope of unity, indicating diffusive motion, and the self-diffusivity of atoms is excellently captured by the corresponding inherent structure curve. At 38 K, the cage effect persists for most of the observation time (Figure 3), whereas at 9 K, complete entrapment is expected

10624 J. Phys. Chem. B, Vol. 112, No. 34, 2008

Tsalikis et al.

Figure 5. Representation of the potential energy per atom at the minimum during the cooling process of the system. The potential energy per atom and time are given in reduced units. In the inset diagram, the evolutionary course of stepwise cooling, involving the three studied temperature levels, is presented.

for times much longer than the observation time (Figure 4). We can conclude that for long times, the diffusivity can be determined equally well by considering transitions of the system from an inherent structure into another one, instead of the complete trajectory on the potential energy hypersurface. This has been pointed out by Shell et al.32 The last conclusion can be mathematically proven as follows: The self-diffusivity corresponds to the long-time gradient of 〈r2(t)〉, that is to say:

6Ds ) lim

t*f∞

∂〈r2 〉 ∂〈(r(t) - r(0))2 〉 ) lim * ∂t ∂t t f∞

(4)

However,

r(t) ) rm(t) + δr(t) ) rm(t) + (r(t) - rm(t))

(5)

and, analogously,

r(0) ) rm(0) + δr(0) ) rm(0) + (r(0) - rm(0))

(6)

The quantity δr(t) is a distance vector from the atom’s position rm in the inherent structure at time t to its instantaneous position r along the system trajectory at the same time. By substituting the quantities, we obtain 2 ∂〈(r(t) - r(0))2 〉 ∂〈(rm(t) + δr(t) - rm(0) + δr(0)) 〉 ) ∂t ∂t

(7)

The quantities δr(t) and δr(0) are instantaneous excursions in the course of oscillations about the inherent structures in which the system finds itself at times t and 0, respectively. They are bounded and comparable in magnitude. On the contrary, rm(t) grows with t in an unbounded fashion, if the motion is diffusive. Thus, |δr(t) - δr(0)| becomes insignificant in relation to |rm(t) - rm(0)| at long times in a diffusive system. Consequently,

∂〈(rm(t) + δr(t) - rm(0) + δr(0))2 〉 6Ds ) lim ) ∂t t*f∞ ∂〈(rm(t) - rm(0))2 〉 (8) lim ∂t t*f∞ From Figure 2, it is evident that the long-time diffusive constant can be calculated either by the full dynamics or by the interbasin transition process mapped to the inherent structure minima. Figure 5 presents the instantaneous values of potential energy at the minima, for three temperatures, in five computational experiments that were started with different initial conditions.

In the stepwise cooling process, the first temperature is 67 K, the second is 38 K, and the last is 9 K. The system spends 14.5 ns at each temperature level. The cooling rate is thus 2 × 109 K/s. For each of the five sets of runs, the final configuration at each temperature was used as starting configuration for the lower temperature. For the highest temperature of 67 K, one can see that all five computational experiments present the same behavior. Transitions from one minimum to another are very fast, and there is quite a variance among the energies of the minima. It can be observed that the number of visited minima of the potential energy is very large, leading to an almost continuous picture. The simulations at 67 K belong in the first characteristic region of our system. For this temperature, the observation time of the system is longer than the time needed in order to explore the local landscape of the supercooled state, resulting in a local equilibrium sampling between nearby basins. At this temperature, many transitions can occur from minimum to minimum. Furthermore, we are in a position to monitor them and accumulate statistics on them, which is independent of the starting configuration over the time of observation. For any temperatures higher than 67 K, we expect a similar behavior. The second region in Figure 5 (14.5 ns at 38 K, close to the expected value for the glass transition temperature Tg, 38.4 K, or 0.32 in reduced units)27 constitutes a passage from the temperature of 67 K, higher that the expected value of Tc (∼52.2 K, or 0.435 in reduced units),26 to the low temperature of 9 K, where our system is deep in the glassy state. In this second region, the system goes from the supercooled liquid state to the glassy state. The observation time of the simulation is similar to the time that the system needs in order to travel in the local landscape of the supercooled state. The mean inherent structure energy along each of the five trajectories is evidently timedependent, and differences between the different trajectories are clearly significant in comparison to the fluctuations along each trajectory. In other words, the trace of inherent structure energy as a function of time is different, depending on the initial configuration. We observe that, even for temperatures lower than Tc, there exists mobility from minimum to minimum. Moreover, the distribution among different minima relaxes over times similar to, or even longer than, the time scale of observation in the computer experiments. This second region is of great importance, because it determines the quality, and hence the properties, of the glass that we will obtain in the third region. This happens for the simple reason that in the third region of 9 K, the system practically freezes at the point where it was found just before

Vitrification Process the temperature drop to 9 K and remains there for a long time. In this third region, the time of observation of the system obtained by using MD is much smaller than the characteristic relaxation time. The five inherent structure trajectories exhibit minor fluctuations around what appears to be constant energy values and show no tendency to approach each other over the time of observation. As has been shown in the work of Kob and Andersen,24 the msd of atoms in such a system should have a specific form when plotted against simulation time. For temperatures significantly higher than Tc, there are two different regimes that can be observed. In particular, for small times, the ballistic regime is observable, followed by the diffusive regime, in which all particles diffuse freely in the melt. For lower temperatures, in the neighborhood of Tc, there is an extra regime which describes what happens at intermediate times, as the system goes from the ballistic to the diffusive regime. The extra regime describes the so-called cage effect, which is an intermediate state between the ballistic and the diffusive regime. The entrapment of the system for these times results in a leveling off of the msd, before it starts rising again by diffusion. As temperature falls, the cage effect regime becomes more pronounced. From Figure 2, one can see that the diffusive dynamics of the system, for the high temperature and long times, can be described completely via the dynamics of inherent structures. The next step was to model the dynamics of transitions between inherent structures on the basis of a first-order kinetic model. According to this model, the succession of transitions from one state (basin surrounding an inherent structure) to any of its neighboring states constitutes a Poisson process. Each transition from an inherent structure β into another one γ is characterized by a transition rate constant kβfγ, whereas the total of inherent structures and transitions form a network governed by firstorder kinetics. This network is analogous to a network of firstorder, reversible chemical reactions, with basins or states around the inherent structures playing the role of chemical species. When the system remains trapped in a small number of inherent structures, as happens at temperatures smaller than the glass transition temperature, it is possible to derive an analytical solution for the evolution of the probability distribution of the system among states, given the transition rate constants kβfγ. The ability to identify new states to be accessed and calculate kβfγ on the fly can be invoked in order to build an ever expanding network of states and use it to evaluate analytical averages over all possible dynamic outcomes of the Poisson process up to a given observation time. This has been achieved by Boulougouris and Theodorou33 via the development of Dynamic Integration over a Markovian Web (DIMW), which relies upon a distinction between explored states and boundary states and invokes absorbing boundary conditions to evaluate analytically the mean first passage time to reach any of the boundary states. The augmentation of the network (by appending a boundary state to the set of explored states) is accomplished by a stochastic selection scheme based on the flux to each of the boundary states. DIMW can be also viewed as an extension to KMC of the Integration over a Markovian Web (IMW) method, introduced by Boulougouris and Frenkel,34 which allows, via rigorous statistical mechanics formulation, the enrichment of ensemble averages and the combination of multiple integration levels in Monte Carlo simulations. In this work, we do not use an on the fly calculation of the rates, as in ref 33. On the contrary, our purpose is to follow the time evolution of the system via full-blown MD and then evaluate the transition rate constants kβfγ on the basis of a

J. Phys. Chem. B, Vol. 112, No. 34, 2008 10625

Figure 6. Dynamic sequence of local minima, as resulting from MD simulation of the LJ mixture (T ) 9 K), where each inherent structure is characterized by its minimum potential energy.

TABLE 2: Transition Rate Constants kβfγ between the States Visited by One Trajectory at 9 K (ps-1) kβfγ (ps-1)

1

1 2 3 4

0.438 0.212 0.222

2

3

4

0.083

0.003 0.000

0.001 0.000 0.106

0.000 0.000

0.444

postprocess hazard plot analysis. An extension to an on-the-fly dynamical sampling on a Markovian web belongs to our future plans. In our model, the most important dynamical quantity is the conditional probability that we begin from a state R at time zero and end up in a state β at time t (Pt ) 0,Rft ) τ,β). This is given via the time integration of the master eq 1. Figure 6 presents the sampling of local minima at T ) 9 K, as obtained from MD simulation after minimization of the potential energy. Each inherent structure is characterized by its minimum potential energy. It is seen that the dynamics at 9 K involves four states. In Table 2, the transition rate constants kβfγ are listed, as they have been evaluated on the basis of the proposed method: Some of the conditional probabilities Pt ) 0,Rft ) τ,β are compared in Figure 7 as calculated on the basis of the solution of eq 1 and as computed directly from an analysis of the MD trajectory.

Figure 7. Conditional probabilities Pt ) 0,Rft ) τ,β, calculated on the basis of eq 1 and of direct analysis of an MD trajectory at T ) 9 K. In the diagram, from top to bottom, we report Pt ) 0,1ft ) τ,1, Pt ) 0,1ft ) τ,2, Pt ) 0,1ft ) τ,3, and Pt ) 0,1ft ) τ,4. The results from the MD (noisy lines) are compared with the results from the model (straight lines) for every one of the four transitions (four levels). As in the case of msd, the sampling of MD at longer times is extremely poor (e.g., for the longest time, there is only one measurement), and this leads to great noise.

10626 J. Phys. Chem. B, Vol. 112, No. 34, 2008

Figure 8. Mean square atomic displacement (in reduced units) based on inherent structures at 23 K. Comparison of MD and first-order kinetic model (eq 9) obtained by using two time intervals between minimization: 1.2 ps and 2 ps. The blue line corresponds to the kinetic model for a time interval of 2 ps, the purple line corresponds the kinetic model for a time interval of 1.2 ps, the red line corresponds the MD simulation for a time interval of 2 ps, and the green line corresponds to the MD simulation for a time interval of 1.2 ps.

Tsalikis et al.

Figure 9. Mean square atomic displacement (in reduced units) based on inherent structures at 9 K. Comparison of MD and first-order kinetic model (eq 9) obtained by using a time interval of 2 ps between minimizations. The purple line corresponds to the kinetic model, and the blue line corresponds the MD simulation.

The knowledge of probabilities Pt ) 0,Rft ) τ,β, in combination with the Cartesian position vectors of atoms in inherent structures (rRi ), leads to a direct calculation of the msd of atoms via eq 9 N

(∆r)2(τ) )

∑ (ri(τ) - ri(0))2 i

N

∑∑ β

R

|

)

N

∑ (riβ(τ) - riR(0))2 i

N

|

Pt)0,Rft)τ,β PR (9)

where i labels the atoms, R and β label the inherent structures, and PR is the probability of being in state R at the beginning of the measurement. Figure 8 presents the coincidence between the calculation of msd of atoms in the inherent structures, as this results from the MD simulation and from eq 9, by using different time intervals between the minimizations: 1.2 ps and 2 ps. Similar comparisons between MD and the first-order kinetic model at T ) 9 K and T ) 38 K, the latter being very close to the glass temperature, are give in Figures 9 and 10, respectively. Figures 7–10 show a very good agreement between the firstorder kinetic model and the results of MD for the entire temperature range from 9 to 38 K, over the time scales investigated. Furthermore, in Figure 8, we show that the resulting msd is insensitive to the choice of the time interval between minimizations within the tested range. One is led to the conclusion that the transitions from inherent structure to inherent structure can be described well with Poisson statistics at sub-Tg temperatures. The hazard plot method was used to map the dynamical behavior obtained from MD simulations onto a first-order kinetic model. Results have shown that this method of analysis is particularly satisfactory. Conclusion In this paper, we presented a numerical investigation of the role of inherent structures in the vitrification process of glassforming materials in the vicinity of and below the glass transition

Figure 10. Mean square atomic displacement (in reduced units) based on inherent structures at 38 K. Comparison of MD and first-order kinetic model (eq 9) obtained by using a time interval of 2 ps between minimizations. The purple line corresponds to the kinetic model, and the blue line corresponds the MD simulation.

temperature. In most previous studies, the inherent structure picture has been used as a tool in order to understand molecular motion in supercooled liquids far above the glass transition temperature Tg. In contrast to these studies, this work has focused on temperatures in the vicinity of and below Tg. By using the inherent structure picture, we have shown that in the vicinity of the glass temperature, the history of the glass-forming material cannot be uniquely defined by the cooling rate. As temperature decreases, approaching Tg, our fixed time window of observation (dictated by the time scales accessible to MD) becomes comparable to the time necessary for equilibrium sampling. Under these conditions, the glassy state in which the system will eventually be trapped depends on the length of time it has spent in the range of temperatures where the simulation time is smaller but comparable to the equilibration time. For temperatures below Tg, we have envisioned a model that describes the dynamics of the atomistic system as a series of transitions, following first-order kinetics and Poisson statistics, from one potential energy minimum (inherent structure) to another. Our validation of this model is based on its ability to reproduce the mean square atomic displacement between potential energy minima resulting from the original atomistic MD trajectory. The projection of the trajectory onto the Poisson process model has been made possible thanks to the use of hazard plot analysis, which allowed the evaluation of the rate constants for the transitions between the inherent structures. Our

Vitrification Process results for a model two-component LJ mixture clearly show that our Poisson process approach is able to reproduce the inherent structure dynamics obtained directly from MD over a wide range of sub-Tg temperatures, from 9 to 38 K. At higher temperatures, the long time diffusion of the system can be described simply by the hopping rate between inherent structures. Furthermore, in the case of vitrification via step cooling, we identify three temperature ranges around the expected glass transition temperature, Tg, of our model system. In the first range, high above Tg, MD simulation is capable of providing an ergodic sampling of the potential energy landscape with relatively low computational cost. In the second range, in the vicinity of the expected Tg, MD is not able to sample the potential energy landscape efficiently. The properties of the amorphous solid obtained upon further reduction of the temperature depend on the time the system is allowed to spend in the second range. In the third range of temperatures, far below Tg, the system is trapped (as expected) in the vicinity of the states where it found itself at the beginning of this range, for times far beyond the reach of MD simulations. Interestingly, the cooling rate alone is not sufficient for defining the glass resulting from a stepwise cooling process; a physically very relevant aspect of the glass formation history is the length of time that the system is left in the intermediate temperature range, because the relaxation time is a strongly nonlinear function of temperature; the longer the system is left at high temperatures (within the second range), the deeper the resulting glass. On the other hand, a continuous decrease of temperature would result in a less equilibrated glass for the same cooling rate. Finally, in the continuation of this work, we will show that the proposed methodology can be extended in order to account for the vibrational contribution of the visited inherent structures;1 we will rigorously lift, within the Poisson approximation, the coarse-grained inherent structure dynamics to restore mean square atomic displacement in the full atomistic system. 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.

J. Phys. Chem. B, Vol. 112, No. 34, 2008 10627 References and Notes (1) Tsalikis, D. G.; Lempesis, N.; Boulougouris, G. C.; Theodorou, D. N. J. Phys. Chem. B 2008, 112, 10628. (2) Gebremichael, Y. Spatially heterogeneous dynamics and string like correlated motion in supercooled liquids and polymers. Ph. D. Dissertation, University of Maryland, MD, 2004. (3) Angell, C. A. J. Non-Cryst. Solids 1988, 102, 205. (4) Dawson, K. A.; Foffi, G.; Sciortino, F.; Tartaglia, P.; Zaccarelli, E. J. Phys.: Condens. Matter 2001, 13, 9113. (5) Boulougouris, G. C.; Frenkel, D. J. Chem. Phys. 2005, 122, 244106. (6) Go¨tze, W.; Sjo¨gren, L. Rep. Prog. Phys. 1992, 55, 241. (7) Debenedetti, P. G.; Stillinger, F. H. Nature 2001, 410, 6825. (8) Kob, W. J. Phys.: Condens. Matter 1999, 11, R85. (9) (a) Go¨tze, W. Aspects of Structural Glass Transitions. 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. (10) Goldstein, M. J. Chem. Phys. 1969, 51, 3728. (11) Sciortino, F. J. Stat. Mech. 2005, P05015. (12) (a) Wang, C.; Stratt, M. R. J. Chem. Phys. 2007, 127, 224503. (b) Wang, C.; Stratt, M. R. J. Chem. Phys. 2007, 127, 224504. (13) Stillinger, F. H.; Weber, T. A. Phys. ReV. A 1982, 25, 978. (14) (a) Theodorou, D. N.; Suter, U. W. Macromolecules. 1985, 18, 1467. (b) Theodorou, D. N.; Suter, U. W. Macromolecules 1986, 19, 379. (c) Kopsias, N. P.; Theodorou, D. N. J. Chem. Phys. 1998, 109, 8573. (15) Stillinger, F. H. Phys. ReV. B 1990, 41, 2409. (16) (a) Bu¨chner, S.; Heuer, A. Phys. ReV. Lett. 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. (17) Helfand, E. J. Chem. Phys. 1978, 69, 1010. (18) Zwanzig, R. J. Chem. Phys. 1960, 33, 1338. (19) Mori, H. Prog. Theor. Phys. 1965, 33, 423. (20) Helfand, E.; Wasserman, Z. R.; Weber, T. A. Macromolecules 1980, 13, 526. (21) Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geisler, P. L. Annu. ReV. Phys. Chem. 2002, 53, 291. (22) Voter, A. F. Phys. ReV. Lett. 1997, 78, 3908. (23) Voter, A. F.; Doll, J. D. J. Chem. Phys. 1985, 82, 80. (24) Kob, W.; Andersen, H. C. Phys. ReV. Lett. 1994, 73, 1376. (25) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Oxford University Press, 1987. (26) Kob, W. Phys. ReV. E 1995, 51, 4626. (27) Shell, S. M.; Debenedetti, P, G.; Panagiotopoulos, A. Z. Fluid Phase Equilib. 2006, 241, 147. (28) Melchionna, S.; Ciccotti, G.; Holian, B. L. Mol. Phys. 1993, 78 (3), 533. (29) (a) Press, W. H.; Teukolsky, S. A.; Vettering, W. T.; Flannery, B. P. Numerical Recipes: The Art of Scientific Computing, 3rd ed.; Cambridge University Press, 2007. (b) Theodorou, D. N. Principles of molecular simulation of gas transport in polymers. In Materials Science of Membranes for Gas and Vapor Separation; Yampolskii, Y., Pinnau, B. D., Eds.; John Wiley: Hoboken, NJ, 2006, pp 47-92. (c) Kaznessis, Y. N. Chem. Eng. Sci. 2006, 61, 940. (30) Reichl, L. E. A Modern Course in Statistical Physics, 2nd ed.; Wiley: New York, 1998; p 229, Chapter 5. (31) Sriraman, S.; Kevrekidis, I. G.; Hummer, G. J. Phys. Chem. B 2005, 109, 6479–6484. (32) Shell, S. M.; Debenedetti, P. G.; Stillinger, F. H. J. Phys. Chem. B 2004, 108, 6772–6777. (33) Boulougouris, G. C.; Theodorou, D. N. J. Chem. Phys. 2007, 127, 084903. (34) Boulougouris, G. C.; Frenkel, D. J. Chem. Theory Comput. 2005, 1, 389–393.

JP801296K