Classification of solid and liquid structures using a deep neural network unveils origin of dynamical heterogeneities in supercooled liquids Open Access

The glass transition occurs in a wide range of substances, including molecules, colloids, granular matter, and biological materials. Over the years, researchers have conducted numerous studies, including experimental, numerical, and theoretical work, to understand the mechanism behind this transition.^1–9 When a liquid is rapidly cooled toward the glass transition, its dynamics slow down significantly, eventually forming a disordered solid state known as glass.^2,3 Even a slight decrease in the temperature leads to a significant increase in relaxation time and a significant decrease in the diffusion coefficient.¹ Interestingly, there is only a minor change in the structure, as indicated by two-point density correlation functions, such as a radial distribution function and a static structure factor.^4,10 This is in contrast to other phase transitions, such as crystallization, where arrested dynamics accompany a significant change in the structure, i.e., the emergence of a long-range periodic ordered structure.

It is a well-known fact that slow dynamics are associated with dynamical heterogeneities.^11–15 As we decrease temperature, we observe an increase in spatial heterogeneities in the dynamics, which results in more pronounced differences between faster and slower regions. Previous studies have attempted to understand the origin of dynamical heterogeneities^5,7,8 and have particularly worked toward establishing connections between structures and dynamics.^16–20 The role of structures in dynamical heterogeneities and slow dynamics has been actively investigated for decades. To this end, we need to characterize the amorphous structures, and several different structural order parameters have been proposed so far. Examples of such order parameters include those that characterize icosahedral bond-orientational orders,^18,21,22 crystal-like structures,^16,17,19,23 and low-energy topological clusters.²⁰ Recently, new order parameters have been proposed to detect entropically stable clusters.^24,25 These structures have often been referred to as locally preferred structures. Furthermore, the free volume has been introduced,^26–28 which has been shown to correlate with dynamical heterogeneities.²⁹ The dynamics are strongly related to the local free volumes, such as the local-average free volume within the neighboring shells, and the quasivoid composed of several neighboring free volumes.^30,31

However, a generalized order parameter adapted to all glass formers has never been found. The main issue is that we need to employ a suitable order parameter depending on the target system, which is usually unknown beforehand and required to be changed heuristically. For some systems, such as the Gaussian-core models,³² we have not discovered an appropriate order parameter yet. Also, some dynamical theories, such as the dynamical facilitation theory³³ based on the kinetically constrained models,³⁴ can also account for glass transitions.^35–39 Therefore, it has not been settled whether the glass transition phenomenon is a thermodynamic or dynamic phase transition.

When encountering difficulties in describing amorphous structures, researchers turn to machine-learning techniques for assistance. Pioneering studies^40,41 have introduced a machine-learning-based order parameter called “softness,” which is strongly related to dynamical behaviors. Since then, researchers have extensively used machine-learning techniques for extracting structural features and understanding their correlations with dynamics in amorphous systems.^42–57 For instance, Refs. 44 and 47 have used support-vector machines to distinguish more mobile particles in glass formers by training data on dynamics and structures. Recently, structural identification and high precisions for predicting dynamical heterogeneity have been realized by supervised approaches,^48,52 which need dynamic data for training. Additionally, works based on unsupervised algorithms without training any dynamic information were also reported to be effective.^49,50

In many of the machine-learning techniques introduced so far, evidence is mounting for the existence of a static structure in glass that leads to dynamic heterogeneity. However, the specific characteristics of these structures are far from clear. In contrast, our recent study⁵⁸ can extract structural features that are easily interpretable. In Ref. 58, we proposed a method that automatically extracts the characteristic structures of amorphous systems based solely on static configurational information. In this method, we first train a convolutional neural network (CNN) to classify configurations of liquid and solid (glass) states. We then utilize the gradient-weighted class activation mapping (Grad-CAM) method^59,60 to identify which part of the system is responsible for the classification. This examination results in Grad-CAM scores, $Γ L$ and $Γ G$⁠, which quantify the degree of liquid-like and solid-like structures in every local part of the system, respectively. Higher values of $Γ L$ indicate more liquid-like structures, while higher values of $Γ G$ indicate more solid-like structures. This approach relies only on structural data and requires no information on dynamics. Additionally, the parameters $Γ L$ and $Γ G$ have clear physical meanings as indicators of liquid-like and solid-like structures, respectively.

The previous work⁵⁸ focused on the non-equilibrium aging dynamics at very low temperatures; here, we conduct a comprehensive study on the equilibrium dynamics for $T L≥T> T G ∗$⁠, where $T G ∗$ is the empirical glass transition temperature. In the present work, we build upon our recent work⁵⁸ and demonstrate the effectiveness of using the CNN and Grad-CAM methods for analyzing the range of temperatures between liquid and glass states used for the CNN training. We train the CNN model using data from two temperatures, namely, the liquid-state temperature, $T L$⁠, and the glass-state temperature, $T G$⁠. We then show that the resulting Grad-CAM scores can accurately extract liquid-like and solid-like structures at any given temperature $T$ within the range of $T L≥T> T G ∗$⁠. Our findings show that the system becomes increasingly solid-like as the temperature decreases, with stronger correlations between the scores and dynamical heterogeneities. This suggests that structures play a more significant role in the dynamics of lower-temperature systems.

This paper is organized as follows. In Sec. II, we introduce the system studied and its properties. In Sec. III, we explain the structural analysis methods based on CNN and Grad-CAM. In Sec. IV, we discuss the static structures characterized by Grad-CAM scores and Voronoi volume $Υ$⁠, and their spatial correlations with the dynamical heterogeneity. Finally, a summary and concluding remarks are given in Sec. V.

In this section, we describe our numerical system. We also present data on the static structure, slow dynamics, and dynamical heterogeneities as the temperature decreases. These data show that our simulations reproduce structural and dynamical properties, which have been observed in many previous works.

We performed molecular-dynamics (MD) simulations using the open-source software LAMMPS,^63,64 with periodic boundary conditions in the $x$ and $y$ directions. We first generated random configurations (shown in the inset of Fig. 1) in the liquid state at a very high temperature $T=4.0$ and equilibrated the system. We then cooled the system toward $T= T G=0.05$ using a constant cooling rate of $dT/dt=8.33× 10 − 5$⁠. Black open circles in Fig. 1 present the temperature dependence of the potential energy per particle in the inherent structure during this cooling process.⁵⁸ It shows that the glass transition occurred at a temperature of $T G ∗≃0.37$⁠. In the present work, we generate the equilibrium liquid states for temperatures in the range of $T L=0.8≥T> T G ∗≃0.37$⁠, as indicated by green closed circles in Fig. 1. At each temperature, we fully equilibrated the system by performing an $NVT$ simulation with the Nos $e ´$–Hoover thermostat to maintain the temperature.⁶⁵ We prepared $50$ independent configurations to perform the ensemble average of physical quantities.

In Fig. 3, we provide data on $⟨Δ r 2(t)⟩$ in (a) and $F s(q,t)$ in (b). The MSD displays ballistic behavior at short times, $⟨ Δ r 2 ( t ) ⟩∝ t 2$⁠, and diffusive behavior at long times, $⟨ Δ r 2 ( t ) ⟩≃4Dt$⁠, where $D$ is the diffusion constant. In between these behaviors, there is an intermediate-time plateau regime, which becomes significantly extended as the temperature is lowered, indicating slow dynamics in the supercooled liquids. $F s(q,t)$ also shows the two-step relaxation behavior with the plateau regime corresponding to the behavior of MSD. We measure the $α$ relaxation time $τ α$ from $F s(q=7.28,t= τ α)=1/e$⁠.

In Fig. 4, we display the values of $D$ in (a) and $τ α$ in (b) as a function of the inverse of temperature $1/T$⁠. As the temperature decreases, $D$ drastically decreases, while $τ α$ drastically increases. Specifically, as $T$ decreases from $T=0.8$ to $0.38$⁠, $D$ decreases from approximately $10 − 3$ to $10 − 6$⁠, while $τ α$ increases from $10 1$ to $10 4$⁠: both of $D$ and $τ α$ change by almost four orders of magnitude. Our simulations, thus, reproduce the slow dynamics in supercooled liquids as in many previous simulations.^1,61,70,71

Figure 5 plots $χ 4(t)$ as a function of time $t$ for various temperatures. We confirm that $χ 4(t)$ reaches a maximum value around the relaxation time $t= τ α$⁠. The height of this peak evaluates the size of clusters of mobile particles. We pick up this peak value as $χ 4 p$ and plot $χ 4 p$ as a function of $1/T$ in the inset of Fig. 5. As the temperature is lowered, $χ 4 p$ increases, which indicates that dynamical heterogeneity is developed. At a temperature of $T=0.39$⁠, we observe that the growth of $χ 4 p$ stops and becomes saturated. This behavior of $χ 4 p$ is in agreement with previous studies.^78,79 Our simulations, thus, reproduce the dynamical heterogeneities in supercooled liquids.

In this section, we outline our machine-learning approach, which combines the CNN and Grad-CAM methods.⁵⁸ This approach generates spatial distributions of Grad-CAM scores, $Γ L$ and $Γ G$⁠, which assess the presence of liquid-like and solid-like structures, respectively. Higher $Γ L$ values indicate a higher proportion of liquid-like structures, while higher $Γ G$ values indicate a higher proportion of solid-like structures. Additionally, we examine the Voronoi volume $Υ$⁠. By utilizing $Γ L$⁠, $Γ G$⁠, and $Υ$⁠, we characterize the structural heterogeneities in the system.

The aim of this work is to identify structural heterogeneities responsible for dynamical heterogeneities. Previous works have reported that the inherent structure underlying the equilibrium liquid state is more correlated with dynamics than the equilibrium configuration (or an instantaneous structure).^41,80 Therefore, the present work focuses on the inherent structure and investigates its spatial heterogeneities and correlations with dynamical heterogeneities.

The inherent structure represents the local minimum in the potential energy landscape of the system. Starting from the equilibrium configuration, we derive the inherent structure by minimizing the total potential of the system using the steepest descent algorithm.⁸¹ We obtain the inherent structure underlying the equilibrium liquid state at each temperature in the range of $0.8≥T> T G ∗≃0.37$⁠, as indicated by the green closed circles in Fig. 1. We then analyze the obtained inherent structure using the CNN and Grad-CAM methods described below.

In our previous work,⁵⁸ we constructed a CNN using two datasets of liquid configurations at $T L=0.8$ (represented by the red cross in Fig. 1) and glass solid configurations at $T G=0.05$ (yellow cross). Note that these configurations are instantaneous ones, not inherent structures. To achieve this, we prepared $5000$ independent configurations, which include $4000$ for training, $400$ for validation, and $600$ for test data, for each state at $T L=0.8$ and $T G=0.05$ (10 000 samples in total). Please see also the caption of Fig. 1.

We utilized a CNN with an architecture identical to the one discussed in Ref. 58. The network does not contain any pooling layers, and it consists of three stacked convolutional layers. Each convolutional layer incorporates a Rectified Linear Unit (ReLU) as the activation function. Following the three convolutional layers and subsequent ReLU activations, a fully connected layer and a dropout layer are introduced before the output layer. For further details and specific choices of hyperparameters, please refer to the supplementary material of Ref. 58. Note that when utilizing the particle configuration data generated by MD simulations as input to the CNN, a projection operator $M$ composed of Gaussian kernels was applied to map the data onto the grids: $x ~= M(x)$⁠, where $x= ∑ i = 1 N x iδ( r− r i)$ and $δ( r)$ represents Dirac’s delta function. Hereafter, we express the data on grids using tilde symbols.

One unique feature of this study is the extraction of the characteristic structures of glasses and liquids, which are crucial for the classification task. This extraction is enabled by the Grad-CAM.^59,60 In this method, the structures that significantly contribute to the classification are extracted by averaging all feature maps in a weighted manner according to the magnitude of their contribution.

Hereafter, the term “Grad-CAM scores” refers to the particle-based values that are re-projected onto each particle: $Γ= ∑ i = 1 N Γ iδ( r− r i)$⁠, where $Γ i$ represents the Grad-CAM score of particle $i$⁠. The calculation of $Γ$ involves the use of the inverse projection operator $M − 1$ from grids to particles, defined as $Γ≡ M − 1( L ~ C)$⁠. This method can identify glass-like structures independent of the classification results or temperatures. In other words, it is capable of extracting glass-like local structures that gradually develop as the temperature decreases in the supercooled liquid state. We refer to the Grad-CAM scores of glass-like and liquid-like structures as $Γ G$ and $Γ L$⁠, respectively. These can be calculated by specifying $C= glass$ or $C= liquid$ in Eqs. (6) and (7).

In this study, we investigate dynamical heterogeneities, specifically examining spatial correlations with the static structure. As previously indicated in Sec. III, we stress that our approach focuses on the underlying inherent structure and its spatial correlations with dynamics. The inherent structure is characterized by spatial distributions of Grad-CAM scores $Γ L$ and $Γ G$ and Voronoi volume $Υ$⁠.

We display the spatial distribution of the dynamic propensity $U( τ α)$ in Fig. 6(d). Note that we normalize $U i$ within the range of $0$– $1$⁠. We clearly observe spatial heterogeneities in $U( τ α)$⁠. In the following, we will examine the correlations of $U( τ α)$ with the structural order parameters of $Γ L$⁠, $Γ G$⁠, and $Υ$⁠.

We now examine the Grad-CAM scores. Figure 6 displays the spatial distributions of $Γ G$ and $Γ L$ for three representative temperatures in panels (a) and (b). It is important to remember that $Γ G$ and $Γ L$ measure solid-like and liquid-like structures: higher values of $Γ G$ indicate more solid-like structures, while higher values of $Γ L$ indicate more liquid-like structures. We also note that $Γ G$ and $Γ L$ are particle-based quantities assigned to each particle, while the figure shows the CG values obtained through the method explained in Sec. E. Here, we employed the same CG length $ξ=4$ for all plots of structure indicators. From Fig. 6, we observe that $Γ G$ and $Γ L$ present the heterogeneities in the static structure, where solid-like and liquid-like structures coexist.

When comparing the values of $Γ G$ and $Γ L$⁠, we observe that they exhibit anticorrelations in their spatial distributions [notice that the direction of the color bar is reversed only in (a)]. Specifically, structures that are more solid-like correspond to less liquid-like structures and vice versa. To quantify this relationship, Fig. 7(a) presents Pearson’s correlation coefficient between $Γ G$ and $Γ L$ as a function of the CG length $ξ$⁠. The figure shows that the coefficient is approximately $−0.7$ for $ξ≥2$ at all temperatures examined. This suggests that $Γ G$ and $Γ L$ convey essentially the same information about the static structure (underlying an inherent structure).

Additionally, we plot the average values of Grad-CAM scores, denoted as $⟨ Γ G⟩$ and $⟨ Γ L⟩$⁠, in Fig. 8, where the average is taken over all samples and particles. The CG operation was not applied here. Because we analyze the inherent structures, the value of $⟨ Γ L⟩$ is much smaller than $⟨ Γ G⟩$⁠. We observe that $⟨ Γ G⟩$ and $⟨ Γ L⟩$ exhibit opposite dependencies on the temperature. Specifically, as the temperature decreases, $⟨ Γ G⟩$ increases monotonically, while $⟨ Γ L⟩$ decreases monotonically. This observation indicates that the system becomes more solid-like and less liquid-like at lower temperatures. It provides evidence that the Grad-CAM scores are a reasonable way to characterize the static structures of supercooled liquids.

We now analyze the spatial correlations between Grad-CAM scores and dynamical heterogeneities. We compare the spatial distributions of the dynamic propensity $U$ and the Grad-CAM scores, $Γ G$ and $Γ L$⁠, as shown in Fig. 6. Especially at the lowest $T=0.38$⁠, we can visually observe spatial correlations between $U(t= τ α)$ and $Γ G, Γ L$⁠. Specifically, we observe that smaller values of $Γ G$ and larger values of $Γ L$ correlate with larger values of $U$⁠, while larger values of $Γ G$ and smaller values of $Γ L$ correlate with smaller values of $U$ [note again that the color bar direction is different only in (a)]. This result demonstrates that particles are more mobile in more liquid-like and less solid-like regimes, and less mobile in more solid-like and less liquid-like regimes.

We conduct a quantitative investigation into the spatial correlations observed above by measuring Pearson’s correlation coefficients, $C Γ G , U$ between $Γ G$ and $U$⁠, and $C Γ L , U$ between $Γ L$ and $U$⁠, with changing the CG length $ξ$⁠. Figure 9 shows the plot of $C Γ G , U$ in (a) and $C Γ L , U$ in (b). We note that $C Γ G , U$ indicates negative correlations, while $C Γ L , U$ shows positive correlations. Maximum values of $− C Γ G , U$ and $C Γ L , U$ are observed at specific CG length $ξ m a x$⁠, which is consistent with the suggestion in Refs. 24 and 25. As suggested by Refs. 24 and 25, this result indicates the nonlocal link between static structural order and dynamical heterogeneity: the correlations emerge only when we coarse-grain the static structure using an appropriate CG length.⁸² Also, the correlations are evident with local structural quantity, such as free volumes over some length scale, rather than with a single-particle structure.^30,31

Based on the data from Fig. 9, we determine the maximum length $ξ α m a x$ at which $C α$ is maximized and its corresponding maximum value $C α m a x$ at each temperature, where $α$ represents $Γ G,U$⁠, or $Γ L,U$⁠. We then plot these values of $ξ α m a x$ in (a) and $C α m a x$ in (b) as a function of $1/T$ in Fig. 10. From panel (a), we observe that $ξ Γ G , U m a x$ and $ξ Γ L , U m a x$ are in agreement, indicating the same length in the static structure. Furthermore, the values of $ξ α m a x$ range from $ξ α m a x=3$ to $5$ and tend to increase monotonically as the temperature decreases. As a result, we deduce that the Grad-CAM scores measure lengthscales of $3$– $5$ in the current temperature range.

In panel (b) of Fig. 10, we observe that both $− C Γ G , U m a x$ and $C Γ L , U m a x$ consistently increase as the temperature decreases. This suggests that the structure becomes more significant to the dynamics as temperature decreases and thermal fluctuations are reduced. At the lowest temperature we investigated, $T=0.38$⁠, the correlation reaches approximately $0.5$⁠, indicating strong correlations with dynamical heterogeneities. Therefore, we can conclude that our approach, which combines CNN and Grad-CAM, appropriately characterizes amorphous structures and predicts dynamical heterogeneities.

At the end of this section, we mention that the Grad-CAM scores for the current MKAM system seem to pick up equivalent information to the Voronoi volume. Figure 6(c) shows the spatial distribution of the Voronoi volume $Υ$⁠, which correlates with the Grad-CAM scores $Γ G$ and $Γ L$⁠. We observe that higher values of $Υ$ correspond to higher values of $Γ L$ and lower values of $Γ G$⁠. Additionally, we calculate Pearson’s correlation coefficients between $Υ$ and $Γ G$ (or $Γ L$⁠) in Fig. 7 and find high correlation values of $0.7$– $0.9$ for $ξ≥2$ at all temperatures studied.

When looking at the spatial correlations of $Υ$ with the dynamic propensity $U(t= τ α)$ in Fig. 6, we observe clear correlations, as in $Γ G$ and $Γ L$ with $U(t= τ α)$⁠. Additionally, as illustrated in Figs. 9 and 10, Pearson’s correlation coefficients between $Υ$ and $U(t= τ α)$ exhibit similar behaviors as do those between $Γ G$ (or $Γ L$⁠) and $U(t= τ α)$⁠. In particular, $C Υ , U m a x$ increases monotonically as the temperature decreases and reaches the value of $C Υ , U m a x≃0.52$ at $T=0.38$⁠. Based on these observations, we conclude that the Grad-CAM scores pick up essentially the same information as the Voronoi volume for the present MKAM system.

We used a combination of CNN and Grad-CAM to study the inherent structures underlying in equilibrium supercooled states. As the temperature of the supercooled state decreases, the Grad-CAM scores, $Γ G$ and $Γ L$⁠, show that the system becomes more solid-like and less liquid-like, respectively. The Grad-CAM scores then reveal the spatial distributions of solid-like and liquid-like structures, which are correlated with the dynamical heterogeneities. Particles in the liquid-like areas are more mobile, while those in the solid-like areas are less mobile. Pearson’s correlation coefficients between structural and dynamical heterogeneities increase as the temperature decreases, reaching approximately 0.5 at the lowest temperature studied in this work. This result indicates that structural heterogeneities have a greater impact on the dynamical heterogeneities at lower temperatures. Based on these findings, we conclude that our method accurately quantifies amorphous structures and predicts heterogeneous dynamics.

In the present work, we focused on the MKAM. Our findings indicate a strong correlation between the Grad-CAM scores in the MKAM and the Voronoi volume. Therefore, our machine-learning-based approach suggests the Voronoi volume as a suitable structural order parameter for the MKAM. However, we highlight that our method is designed to automatically identify the most appropriate structural order parameter depending on the specific target system. In the future, it will be interesting to investigate systems in which suitable structural order parameters have not yet been identified, such as the harmonic potential system⁸⁷ and the Gaussian-core model.⁸⁸ 

Our study conveys that the dynamical heterogeneity may be governed by the static structures. Also, a work by Monte Carlo simulations established a link between the local thermodynamic fluctuations and dynamic fluctuations.⁸⁹ However, recently, it is reported that dynamic facilitation has become more significant in controlling dynamical heterogeneity, especially at very low temperatures.³⁹ The dynamic origin in the dynamical facilitation theory seems to contradict our observations, while Ref. 90 revealed a structural origin of dynamic facilitation, showing that more mobile particles emerge where more free volumes interact, consistent with our results. Indeed, the conundrum of whether the glass transition is a dynamic transition or a thermodynamic phase transition has not been solved yet, and further research is needed in the future.

We thank Atsushi Ikeda and Kang Kim for their useful discussions and comments. H.M. was supported by JSPS KAKENHI (Nos. 22K03543 and 23H04495). T.K. was supported by the JST FOREST Program (No. JPMJFR212T), the AMED Moonshot Program (No. JP22zf0127009), and the JSPS KAKENHI (Nos. 24H02203, 22H04472, 20H05157, and 20H00128). M.L. was supported by the Collaborative Program from the Chinese Scholarship Council (CSC) and the Japanese Government (MEXT) Scholarship (No. 202106370007).

The authors have no conflicts to disclose.

Min Liu: Data curation (equal); Writing – original draft (equal). Norihiro Oyama: Data curation (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Takeshi Kawasaki: Supervision (equal); Writing – review & editing (equal). Hideyuki Mizuno: Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.