Abstract
In contrast with crystallization, there is no noticeable structural change at the glass transition. Characteristic features of glassy dynamics that appear below an onset temperature, T0 (refs 1,2,3), are qualitatively captured by mean field theory4,5,6, which assumes uniform local structure. Studies of more realistic systems have found only weak correlations between structure and dynamics7,8,9,10,11. This raises the question: is structure important to glassy dynamics in three dimensions? We answer this question affirmatively, using machine learning to identify a new field, ‘softness’ which characterizes local structure and is strongly correlated with dynamics. We find that the onset of glassy dynamics at T0 corresponds to the onset of correlations between softness (that is, structure) and dynamics. Moreover, we construct a simple model of relaxation that agrees well with our simulation results, showing that a theory of the evolution of softness in time would constitute a theory of glassy dynamics.
Similar content being viewed by others
Main
To look for correlations between structure and dynamics, one typically tries to find a quantity that encapsulates the important physics, such as free volume, bond orientational order, locally preferred structure, and so on. In contrast to this approach, we use a machine learning method designed to find a structural quantity that is strongly correlated with dynamics. Earlier, we applied this approach to the simpler problem of classifying particles as being ‘soft’ if they are likely to rearrange or ‘hard’ otherwise12. We describe a particle’s local structural environment with M = 166 ‘structure functions’13 that respect the overall isotropic symmetry of the system and include radial density and bond angle information. We then define an M-dimensional space, , with an orthogonal axis for each structure function. The local structural environment of a particle i is thus encoded as a point in M-dimensional space. We assemble a ‘training set’ from molecular dynamics simulations consisting of equal numbers of ‘soft’ particles that are about to rearrange and ‘hard’ particles that have not rearranged in a time τα preceding their structural characterization, and find the best hyperplane separating the two groups using the support vector machines (SVM) method14,15. Finally, we define the softness, Si, of particle i as the shortest distance between its position in and the hyperplane, where Si > 0 if i lies on the soft side of the hyperplane and Si < 0 otherwise.
We study a 10,000-particle 80:20 bidisperse Kob–Andersen Lennard-Jones glass16 in three dimensions at different densities ρ and temperatures T above its dynamical glass transition temperature. All results here are for particles of species A only. However, the results are qualitatively the same for particles of both species. At each density we select a training set of 6,000 particles, taken from a molecular dynamics trajectory at the lowest T studied, to construct a hyperplane in . We then use this hyperplane to calculate Si(t) for each particle i at each time t during an interval of 30,000τ at each ρ and T.
We can deduce the most important structural features contributing to softness either by training on fewer structure functions or by examining the projection of the hyperplane normal onto each orthogonal structure function axis. Both analyses yield a consistent picture (see Supplementary Information): the most important features are the density of neighbours at the first peaks of the radial distribution functions gAA(r) and gAB(r); these two features alone give 77% prediction accuracy for rearrangements. Particles with more neighbours at the first peaks of g(r) have a lower softness, and are thus more stable. These results are reminiscent of the cage picture, in which an increase of population in the first-neighbour shell suppresses rearrangements, or the free-volume picture, in which particles whose surroundings are closely packed are more stable than those with more loosely packed neighbourhoods17. Overall, soft particles typically have a structure that is more similar to a higher-temperature liquid, where there are more rearrangements, whereas hard particles have a structure that is closer to a lower-temperature liquid18.
Figure 1a is a snapshot with particles coloured according to their softness. Evidently, S has strong spatial correlations. Figure 1b shows the distribution of softness, P(S), and the distribution of softness for particles just before they go through a rearrangement, P(S ∣ R). We see that 90% of the particles that undergo rearrangements have S > 0. We have also tested other sets of structure functions (see Supplementary Information) and found nearly identical accuracy. Softness is therefore a highly accurate predictor of rearrangements that is reasonably robust to the set of structure functions chosen.
We next show that the probability that particles rearrange is a function of their softness. This probability is calculated as the fraction of particles of a softness, S, that are rearranging at a given time, PR(S). We plot PR(S) in Fig. 2a in solid lines at temperatures ranging from T = 0.47 (blue) to T = 0.58 (red). At each T we see that PR(S) is a strong function of softness, increasing by several orders of magnitude, especially at the lower temperatures, in the range S = −3 to S = 3. A similar, albeit more modest, relationship was seen in ref. 19. When PR(S) is plotted as a function of 1/T for several values of softness (Fig. 2b), the probability that a particle of softness S will rearrange has Arrhenius behaviour, PR(S) = P0(S)exp(−ΔE(S)/T), where P0(S) and ΔE(S) depend on S. Confirming this observation, PR(S)/P0(S) collapses over many orders of magnitude for all temperatures when plotted against ΔE(S)/T, as shown in the inset of Fig. 2b.
An Arrhenius form emerges when a kinetic process depends on a single energy scale ΔE(S). In Fig. 2c we plot ΔE(S) and Σ(S) ≡ lnP0(S) versus S. Both terms depend nearly linearly on S: ΔE = e0 − e1S and Σ = Σ0 − Σ1S, where all four coefficients are positive and independent of T. Our results are consistent with the interpretation that, at low temperatures, harder regions of the glassy liquid with higher energy barriers are frozen out whereas softer regions are not, leading to heterogeneous dynamics. These heterogeneities smooth out with increasing temperature, and vanish altogether once PR(S) no longer depends on softness. This occurs at the temperature T0 where the softness dependence of Σ exactly cancels that of ΔE/T0, and so T0 = e1/Σ1. This result can also be seen visually in Fig. 2b, where the different Arrhenius predictions for PR(S) all intersect at a single temperature, T0, where the probability of rearrangement will be independent of softness. In Fig. 2d we compare our prediction for T0 to the onset temperature of glassy dynamics measured by Keys et al. 20, T0m, at different densities. The excellent agreement between the predicted T0 and the measured values implies that the onset of glassy dynamics at T = T0 coincides with the onset of correlations between structure (softness) and dynamics.
We explore next the relationship between softness and the non-exponential decay of the overlap function
where N is the number of particles in the system, ri is the position of particle i, and Θ is the Heaviside function. We take a = 0.5 (ref. 21). In Fig. 3a, b we plot the overlap function for different temperatures at ρ = 1.20. Our aim is to understand the form of the decay of q(t) from the behaviour of the rearrangement probability, PR(S). To begin, we define the contribution to the overlap from particles whose softness was initially S at t = 0, q(S, t). The total overlap is q(t) = ∫ dSq(S, t)P(S). Because q(S, t) is the fraction of particles with initial softness S that have not rearranged after a time t, we expect (dq(S, t)/dt)|t=0 = −caPR(S) (see Supplementary Information for details), where ca is the fraction of rearrangements that displace particles by more than a. This is indeed the case, as is evident from the data in Fig. 2a, when (dq(S, t)/dt)|t=0 (points) is overlaid with PR(S) (solid lines).
If we now assume that each particle rearranges with probability PR(S) as an independent Arrhenius process according to Fig. 2, then we can predict the decay of q(S, t) using a simple discrete model: it can be written in terms of the probability that a particle of softness S does not rearrange for t − 1 time steps before finally rearranging at time t, given by (1 − PR(S))t−1PR(S). The resulting prediction for q(t) (dashed) is shown in Fig. 3a for several different temperatures. Although the prediction is not poor, its accuracy decreases at longer times, particularly at lower temperatures.
We now show that the discrepancy between our naive theory and the decay of q(t) primarily results from a crucial neglected feature: even if a given particle does not rearrange, its local structural environment—and therefore its softness—can be altered if nearby particles rearrange. This physics is reminiscent of facilitation.
To take this facilitation into account, we calculate the ‘softness propagator’, G(S, S0, t), defined as the distribution of softness at time t for particles that start with a softness S0 at t = 0 and move less than a distance a after a time t (that is, that do not rearrange in a time t). Figure 4a shows a Gaussian approximation to G(S, S0 = −3, t). We see that G(S, S0, t) is sharply peaked around S0 at small t, but widens and shifts with increasing t, reminiscent of directed diffusion. Figure 4b shows G(S, S0 = −3, t) at several different times, where points are measured probabilities and dashed lines are their Gaussian approximations. In Fig. 4c we plot the mean softness evaluated as = ∫ dSSG(S, S0, t) for several different values of S0. For each S0 the average softness of particles evolves towards the mean of the equilibrium softness distribution over a time period of approximately τα. The softness propagator is evaluated only for particles that have not rearranged, so Fig. 4 shows that rearrangements of nearby particles affect a particle’s softness significantly.
Our first naive prediction based on the assumption that particles rearrange independently corresponds to G(S, S0, t) = δ(S − S0). We refine our theory by using the actual softness propagator in connecting the probability of rearranging, PR(S), with the overlap q(S, t) (see Supplementary Information). For ease of calculation, we approximate G(S, S0, t) as a Gaussian distribution in S and calculate its mean and variance as functions of S0 and t from simulated data. The resulting prediction for the overlap is shown in Fig. 3b. The agreement with the actual q(t) is excellent, suggesting that an understanding of the time evolution of the softness field, or equivalently of the softness propagator, would suffice to understand the non-exponential decay of the overlap function.
Our results show that there is structure hidden in the disorder of glassy liquids. This structure can be quantified by softness, which controls glassy dynamics at temperatures below T0. According to our analysis, simple Arrhenius relaxation for each softness, coupled with the time evolution of softness, leads to the observed slow, non-exponential relaxation dynamics of glassy liquids below T0. Thus, our results suggest that the challenge of understanding glass transition dynamics can be reframed as the challenge of understanding the evolution of softness.
Methods
System information.
We study a 10,000-particle Kob–Andersen model, a 80:20 binary LJ mixture16 with parameters: σAA = 1.0, σAB = 0.8, σBB = 0.88, εAA = 1.0, εAB = 1.5, εBB = 0.5, mA = mB = 1. Time is measured in units of and the Boltzmann constant is kB = 1. We cut off the LJ potential at 2.5σAA and smooth the potential so that force varies continuously. This mixture has been characterized extensively. In particular, we compare our predictions to the measurements of the onset temperature in Keys and colleagues20. Simulations were done using LAMMPS (ref. 22) in an NVT ensemble with a Nosé–Hoover thermostat and a time step of 0.0025τ. We output states every τ and quench them to their nearest inherent structure using a combination of conjugate gradient and FIRE algorithms. Throughout this study we use inherent structure positions. However, qualitatively similar results can be obtained using time-averaged positions. We study this system over the temperatures and number densities listed in Table 1.
Identifying rearrangements.
We adapt a method first proposed by Candelier and colleagues23,24. A timescale tR = 10τ is chosen to be commensurate with the amount of time the system takes to complete a rearrangement. Then two time intervals are defined as A = [t − tR/2, t] and B = [t, t + tR/2]. An indicator function can then be written as
where 〈〉A and 〈〉B are averages over the intervals A and B respectively. phop is large when the mean position of a particle changes appreciably. Otherwise, it is similar in magnitude to the variance in particle positions due to noise from the inherent structure calculation.
To find rearrangements we restrict our attention to events in which phop exceeds a threshold of 0.05, which is large compared to the scale of fluctuations in particle positions but small compared to the typical value of phop during a rearrangement. As discussed in the Supplementary Methods, we define rearrangements to be those events with phop∗ > pc = 0.2. Changing this cutoff affects the results only quantitatively and manifests itself primarily as a shift in the energy scale, ΔE, which is approximately logarithmic in the cutoff. This agrees with the observations of Keys et al. 20, who saw a similar logarithmic shift in the energy scale governing rearrangements with the size of the rearrangements.
Note that rearrangements defined using phop result in particle displacements that follow a distribution that depends on the cutoff pc used. This pc dependence needs to be addressed when comparing the probability of rearrangement to the overlap function and its derivative, which are defined in terms of a length scale a. To do this, we multiply PR by a temperature-independent constant ca, namely the fraction of rearrangements that displace particles by more than a.
Computing softness.
We have made two improvements that greatly increased the prediction accuracy for rearrangements compared to ref. 12. First, we identified rearrangements more carefully, as detailed above. Second, we defined our training sets more carefully. Each training set contains 6,000 particles that rearrange in the next time step, each labelled with yi = 1, as well as 6,000 particles that have not rearranged for a time τα before the structure was calculated, each labelled with yi = 0. These particles were chosen randomly from the set of all particles satisfying these conditions from MD simulations at a low temperature. Then, a training set of N particles can be written as {(F1, y1), …, (FN, yN)}, where Fi = {Fi1, …, FiM} are the M structure functions that describe the local neighbourhood of particle i.
As in ref. 12 we use two classes of functions to generate the set of M structure functions. The first class measures the density of particles a distance r ± δ from a reference particle, i,
where Rij is the distance between particles i and j and X denotes a species whose density we wish to probe. By varying r, δ and X these functions investigate different aspects of the density of particles near the particle i. In this work we keep δ = 0.1σAA fixed. The second class of functions counts the number of small (or large) bond angles between triples of particles within a distance ξ of one another. It is given by
where θijk is the angle between Rij and Rik, λ = ±1 determines whether we consider small or large bond angles, ζ determines the angular resolution, and X and Y are species. By varying X, Y, ζ, λ and ξ we probe different aspects of a particle’s angular neighbourhood.
We then use an SVM to find the hyperplane w ⋅F − b = 0 that separates the points with yi = 1 from those with yi = 0. This hyperplane is used on the rest of the data to reach the results reported. The SVM is trained—that is, the hyperplane is constructed—on the binary variable y using the LIBSVM package14. It is not possible to find a hyperplane that perfectly separates the two different classes. We use a penalty parameter C and find the optimal hyperplane equation by minimizing
with the constraint yi ⋅(wT ⋅Fi + b) ≥ 1 − χi and χi ≥ 0, where χi are the slack variables. The C parameter was chosen through cross-validation12. The hyperplane obtained from this training can be used to classify a new particle neighbourhood, Fn, as soft or hard. Fn is soft if w ⋅Fn − b > 0, and hard otherwise. The continuous variable softness is defined by Sn = w ⋅Fn − b. See Supplementary Information for other classification schemes explored.
References
Angell, C. A. Formation of glasses from liquids and biopolymers. Science 267, 1924–1935 (1995).
Sastry, S., Debenedetti, P. G. & Stillinger, F. H. Signatures of distinct dynamical regimes in the energy landscape of a glass-forming liquid. Nature 393, 554–557 (1998).
Debenedetti, P. G. & Stillinger, F. H. Supercooled liquids and the glass transition. Nature 410, 259–267 (2001).
Parisi, G. & Zamponi, F. Mean-field theory of hard sphere glasses and jamming. Rev. Mod. Phys. 82, 789 (2010).
Berthier, L., Biroli, G., Bouchaud, J.-P., Cipelletti, L. & van Saarloos, W. Dynamical Heterogeneities in Glasses, Colloids, and Granular Media (Oxford Scholarships Online, 2011).
Charbonneau, P., Kurchan, J., Parisi, G., Urbani, P. & Zamponi, F. Exact theory of dense amorphous hard spheres in high dimensions. III. the full RSB solution. J. Stat. Mech. P10009 (2014).
Gilman, J. Metallic glasses. Phys. Today 28, 46–53 (May, 1975).
Berthier, L. & Jack, R. L. Structure and dynamics of glass formers: predictability at large length scales. Phys. Rev. E 76, 041509 (2007).
Royall, C. P., Williams, S. R., Ohtsuka, T. & Tanaka, H. Direct observation of a local structural mechanism for dynamic arrest. Nature Mater. 7, 556–561 (2008).
Manning, M. L. & Liu, A. J. Vibrational modes identify soft spots in a sheared disordered packing. Phys. Rev. Lett. 107, 108302 (2011).
Jack, R. L., Dunleavy, A. J. & Royall, C. P. Information-theoretic measurements of coupling between structure and dynamics in glass-formers. Phys. Rev. Lett. 113, 095703 (2014).
Cubuk, E. D. et al. Identifying structural flow defects in disordered solids using machine learning methods. Phys. Rev. Lett. 114, 108001 (2015).
Behler, J. & Parrinello, M. Generalized neural-network representation of high-dimensional potential-energy surfaces. Phys. Rev. Lett. 98, 146401 (2007).
Chang, C.-C. & Lin, C.-J. LIBSVM: a library for support vector machines. ACM Trans. Intell. Syst. Technol. 2, 27 (2011).
Cortes, C. & Vapnik, V. Support-vector networks. Mach. Learn. 20, 273–297 (1995).
Kob, W. & Andersen, H. C. Scaling behavior in the β-relaxation regime of a supercooled Lennard-Jones mixture. Phys. Rev. Lett. 73, 1376 (1994).
Berthier, L. & Biroli, G. Theoretical perspective on the glass transition and amorphous materials. Rev. Mod. Phys. 83, 587–645 (2011).
Percus, J. K. & Yevick, G. J. Analysis of classical statistical mechanics by means of collective coordinates. Phys. Rev. 110, 1–13 (1958).
Smessaert, A. & Rottler, J. Structural relaxation in glassy polymers predicted by soft modes: a quantitative analysis. Soft Matter 10, 8533–8541 (2014).
Keys, A. S., Hedges, L. O., Garrahan, J. P., Glotzer, S. C. & Chandler, D. Excitations are localized and relaxation is hierarchical in glass-forming liquids. Phys. Rev. X. 1, 021013 (2011).
Keys, A. S., Abate, A. R., Glotzer, S. C. & Durian, D. J. Measurement of growing dynamical length scales and prediction of the jamming transition in a granular material. Nature Phys. 3, 260–264 (2007).
Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117, 1–19 (1995).
Candelier, R. et al. Spatiotemporal hierarchy of relaxation events, dynamical heterogeneities, and structural reorganization in a supercooled liquid. Phys. Rev. Lett. 105, 135702 (2010).
Smessaert, A. & Rottler, J. Distribution of local relaxation events in an aging three-dimensional glass: temporal correlation and dynamical heterogeneity. Phys. Rev. E 88, 022314 (2013).
Acknowledgements
We thank D. Chandler, K. Mandadapu, M. Parrinello, J. Sethna and D. Wales for helpful discussions. This work was supported by the US Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award DE-FG02-05ER46199 (S.S.S. and A.J.L.), the Advanced Materials Fellowship of the American Philosophical Society (D.M.S.), a Harvard IACS Student Scholarship (E.D.C.), and a Simons Investigator award from the Simons Foundation (A.J.L.).
Author information
Authors and Affiliations
Contributions
E.D.C. and S.S.S. contributed equally to this work. S.S.S. and E.D.C. conceived the project, did the simulations, and analysed the data. A.J.L. supervised the project and analysis. All authors contributed to analysis, discussion and writing of the manuscript.
Corresponding authors
Ethics declarations
Competing interests
The authors declare no competing financial interests.
Supplementary information
Supplementary information
Supplementary information (PDF 2918 kb)
Rights and permissions
About this article
Cite this article
Schoenholz, S., Cubuk, E., Sussman, D. et al. A structural approach to relaxation in glassy liquids. Nature Phys 12, 469–471 (2016). https://doi.org/10.1038/nphys3644
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1038/nphys3644
This article is cited by
-
Anomalous temperature dependence of elastic limit in metallic glasses
Nature Communications (2024)
-
Modern computational studies of the glass transition
Nature Reviews Physics (2023)
-
Finding defects in glasses through machine learning
Nature Communications (2023)
-
Relevance of structural defects to the mechanism of mechanical deformation in metallic glasses
Scientific Reports (2023)
-
Entanglement Detection with Complex-Valued Neural Networks
International Journal of Theoretical Physics (2023)