Attractive versus truncated repulsive supercooled liquids: The dynamics is encoded in the pair correlation function

TL;DR

This study shows that subtle differences in pair correlation functions, captured through machine learning, can explain the distinct glassy dynamics in liquids with similar structures but different interaction potentials.

Contribution

It introduces a machine learning-based weighting method to connect pair correlation functions with dynamical differences in supercooled liquids.

Findings

Weighted integral of pair correlation function predicts dynamical differences.

Attractive interactions significantly alter glassy dynamics.

Machine learning effectively captures subtle structural influences.

Abstract

We compare glassy dynamics in two liquids that differ in the form of their interaction potentials. Both systems have the same repulsive interactions but one has also an attractive part in the potential. These two systems exhibit very different dynamics despite having nearly identical pair correlation functions. We demonstrate that a properly weighted integral of the pair correlation function, which amplifies the subtle differences between the two systems, correctly captures their dynamical differences. The weights are obtained from a standard machine learning algorithm.

Equations

$$p_{hop}(i,t)\equiv \sqrt{{⟨{\left({\vec{r}}_i-{⟨{\vec{r}}_i⟩}_U\right)}^2⟩}_V^{1/2}{⟨{\left({\vec{r}}_i-{⟨{\vec{r}}_i⟩}_V\right)}^2⟩}_U^{1/2}},$$

$$G_{Y,\alpha }(i,t)=\frac{1}{4\pi r_{\alpha }^2}\cdot \sum _{j\in \{Y\}}{\mathbb{I}}_{|{\vec{r}}_j-{\vec{r}}_i|\in [r_{\alpha },r_{\alpha +1}]},$$

$$S(i,t)=\vec{w}\cdot \vec{G}(i,t)+b,$$

$$\text{Score}=\frac{\mathbb{P}(S>0|y_{GT}=1)+\mathbb{P}(S<0|y_{GT}=1)}{2},$$

$$⟨S⟩=\vec{w}\cdot ⟨\vec{G}⟩+b,$$

$$P(R|S)=\exp \left[-\frac{\Delta E(S)-T\Delta \Sigma (S)}{T}\right]$$

Full text

The central tenet of liquid state theory is that the structural and thermodynamic properties of liquids are primarily governed by the short-range repulsive part of their interaction potentials Widom (1967); Weeks et al. (1971). The extension of these ideas to dynamical properties, however, is still in question. Studies of glass-forming liquids with purely repulsive potentials show that their relaxation times $\tau_{\alpha}$ are well-described using a perturbative analysis around the hard-sphere limit Haxton et al. (2011); Schmiedeberg et al. (2011). However, a study of the Kob-Andersen Lennard-Jones (LJ) mixture and its repulsive counterpart, the Weeks-Chandler-Andersen (WCA) mixture, found that the two systems have very different dynamics despite their structural similarities Berthier and Tarjus (2009, 2011a). In particular, the relaxation time $\tau_{\alpha}$ is much smaller in the WCA system than in the LJ one at the same supercooled temperatures. While hidden potential scale invariances Pedersen et al. (2010), complex structural indicators such as the point-to-set length Biroli and Bouchaud (2004) and triplet correlation functions markedly differ in these two systems Hocky et al. (2012); Coslovich (2013), the pair correlation functions are instead very similar. These findings indicate that long-ranged attractions strongly affect dynamics while leaving simple structural features relatively unaffected.

This result is a challenge to theories directly based on the exact (numerically computed) thermodynamically-averaged pair correlation functions $\left\langle g(r)\right\rangle$ Berthier and Tarjus (2011b). Mode coupling theory is the most prominent example of such a theory, and indeed it was shown to largely underestimate the difference in dynamics of the WCA and LJ mixtures Berthier and Tarjus (2010). From this negative result, it was argued that $\left\langle g(r)\right\rangle$ does not contain the physical information relevant for predicting the dynamical slowdown, at least not in a way amenable to predictions, and that any theory of the glass transition based on the sole basis of $\left\langle g(r)\right\rangle$ is bound to fail.

The view that $\left\langle g(r)\right\rangle$ alone cannot explain the dynamical difference between the WCA and LJ systems has recently been questioned. Combining it with the force-force correlation function and treating the attractive and slowly varying forces separately, Schweizer et. al. Dell and Schweizer (2015) proposed an alternative microscopic theory for the dynamical differences between LJ and WCA systems. Modeling dynamics as that of a single-particle in an effective caging potential and using the pair-correlations based configurational entropy, Nandi et al. (2017) points toward $\left\langle g(r)\right\rangle$ containing enough information to predict the transition temperature. A comparison of that theory with Schweizer’s has recently shown they produce comparable results, although the latter is more precise Saha et al. (2019). Independently, approximate measures of the configurational entropy based on $\left\langle g(r)\right\rangle$ can be used to estimate the dynamics via the Adam-Gibbs relation, and it was argued that the resulting difference of configurational entropy can explain the dynamical difference between the two systems  Banerjee et al. (2014, 2016, 2017). Although it is unclear why the configurational entropy should play a role at the temperatures studied, which are above the mode-coupling temperature, these results Banerjee et al. (2014, 2016, 2017) suggest that the relationship between $\left\langle g(r)\right\rangle$ and dynamics should be revisited.

In this RRapid Communication we demonstrate that a properly weighted integral of the pair correlation function obtained by machine-learning (ML) techniques Cubuk et al. (2015); Schoenholz et al. (2016); Cubuk et al. (2016); Schoenholz et al. (2015); Cubuk et al. (2017); Sussman et al. (2016), the so-called ”softness” designed to correlate strongly with instantaneous mobility, fully captures the dynamical differences between the LJ and the WCA supercooled liquids, despite their minute differences in average structure. Our results offer a simple explanation of the relationship between $g(r)$ and relaxation time, and show that softness is an important structural fingerprint of dynamics. We emphasize that in contrast to previous works on Softness, here we are able to extract dynamics in systems with very similar $\left\langle g(r)\right\rangle$, showing that fluctuations in $g(r)$ are relevant to the dynamics. We conclude that, at least in principle, one should be able to compute the relaxation time of these two liquids on the sole basis of the pair correlation.

Following Berthier and Tarjus Berthier and Tarjus (2011a), we compare two Kob-Andersen Kob and Andersen (1995a, b) binary mixtures of A and B particles, one with Lennard-Jones interactions and the other with WCA interactions Weeks et al. (1971). The latter system shares the same repulsive short-range interaction with the LJ system, but possesses no attractive forces, i.e. the potential is zero at the minimum $r=r_{cut}=2^{1/6}\sigma_{XY}$ and beyond. The potentials are smoothed at the cutoff distance, so that their second derivatives are continuous everywhere. The two liquids are first equilibrated in the NVT ensemble; we then record particle trajectories $\vec{r}_{i}(t)$ for $\sim 200\tau_{\alpha}$ in the NVE ensemble Anderson et al. (2008); Glaser et al. (2015) and extract the relaxation time $\tau_{\alpha}$ from the decay of the self-intermediate scattering function. Following Ref. Berthier and Tarjus (2011a), we show in Fig. 1b that $\tau_{\alpha}$ increases much more rapidly in the LJ system than in the WCA one when temperature is lowered, while concomitantly the differences in the average pair correlation function $\left\langle g(r)\right\rangle_{AA}$ between systems remain small relative to the fluctuations within each system, as shown in Fig. 1a.

In the supercooled regime particles essentially fluctuate around a given average position, effectively trapped in a cage of surrounding particles. From time to time a rearrangement, involving a small number of particles, occurs. At the level of a single particle $i$, a rearrangement is characterized by a cage-jump, corresponding to a relatively rapid change in the average position of the particle that is on the scale of the particle size. Cage-jumps taking place during a time interval $W=[t_{1},t_{2}]$, for a given particle $i$, may be captured by computing the quantity $p_{hop}(i,t)$ Candelier et al. (2009); Candelier et al. (2010a, b); Smessaert and Rottler (2013):

[TABLE]

for all $t\in W$, where averages are performed over the time intervals surrounding time $t$, i.e. $U=[t_{1},t]$, $V=[t,t_{2}]$. $t^{*}=\text{argmax}_{t\in W}(p_{hop}(i,t))$ is the time that best separates the trajectory $\vec{r}_{i}([t_{1},t_{2}])$ into two sub-trajectories $\vec{r}_{i}([t_{1},t^{*}])$ and $\vec{r}_{i}([t^{*},t_{2}])$. We record $p_{hop}^{*}(i)=p_{hop}(i,t^{*})$ and the process is iterated on each sub-trajectory. If the trajectory during the time interval $W$ is contained within a cage then $p_{hop}^{*}(i)$ is small and corresponds to the cage’s size. Conversely, if $W$ contains two distinct cages then $p_{hop}^{*}$ is large and corresponds to the jump’s amplitude. We thus interrupt the iteration when $p_{hop}^{*}<p_{c}$, where $p_{c}$ is chosen as the root-mean-squared displacement (RMSD) $\Delta r(t)$, computed at the time where the non-Gaussian parameter $\alpha_{2}=\frac{3\left\langle r^{4}\right\rangle(t)}{5\left\langle r^{2}\right\rangle^{2}(t)}-1$ has a maximum. This choice is rather stringent: only displacements several times larger than the plateau of the RMSD are considered as jumps. Finally, as jumps are never instantaneous, we assign $p_{hop}(i,t)=p_{hop}^{*}(i)$ to all times $t\in[t^{*}-t_{f},t^{*}+t_{f}]$, where $t_{f}=5$ in LJ-time units, i.e. about five ballistic times. This procedure generates a label that classifies every particle at every instant as caged or jumping. This constitutes what is called in ML jargon the ground truth: $y_{GT}(i,t)\equiv\Theta(p_{hop}(i,t)-p_{c})\in\{0,1\}$.

If dynamics is indeed related to local structure, then there should be small but significant structural differences between the local neighborhoods of caged and cage-jumping particles, respectively. A number of approaches have attempted to identify these subtle structural signatures in supercooled liquids Coslovich and Pastore (2007); Chen et al. (2010); Zylberg et al. (2017). They all have in common the fact that they are based on a deep knowledge of the physics of liquids. Here, in contrast, we use a ML approach as introduced in Refs. Cubuk et al. (2015); Schoenholz et al. (2016, 2015). The first step is to identify the structural observables that we will try to correlate with dynamics. We focus on $\vec{G}(i,t)$, a vector of $160$ descriptors of the local structure around particle $i$ at time $t$. Considering particles of type $A$, it reads:

[TABLE]

where $Y$ denotes the type of particles $(A,B)$, and $r_{\alpha}~{}=~{}\{0.70,0.75,...,4.70\}$ is an $80$-point binning of the distance separating pairs of particles. Note that $G_{Y,\alpha}(i,t)$ and the instantaneous pair correlation, $g_{X_{i}Y}(r_{\alpha},t)$, are proportional up to a constant that depends on the binning size.

The second step, called training in the context of ML, consists in fitting a function $f(\vec{G}(i,t)),sign(f(G(i,t)))=y_{GT}(i,t)$, that allows to predict unobserved occurences of $y_{GT}$ from the corresponding $\vec{G}(i,t)$. To do so, we use a linear support vector machine Guyon et al. (1996); Pedregosa et al. (2011), i.e. the function $f$ is an affine function of the input vector’s coefficients. An intuitive picture is to consider $f$ as encoding the hyperplane which best separates the two populations $y_{GT}(i,t)=-1$ or $1$, in the 160-dimensional descriptor space. In our context, the generic fitting function $f$ is called softness $S$:

[TABLE]

namely the signed distance between the point $\vec{G}(i,t)$ and the hyperplane defined by $b$ and $\vec{w}$. Concretely, for a given training set, an iterative optimization routine adapts the coefficients $\vec{w},b$ so as to minimize the number of errors on the training examples (an error is when $sign(f(G(i,t))\neq y_{GT}(i,t)$. One says the model generalizes well when the error rate on a test set is as low as the training error rate (the test set it typically drawn from the same distribution as the training set, but completely uncorrelated from it). The softness characterizes the local, instantaneous structure of any particle for the two systems (LJ, WCA), at any temperature. Note that $\vec{G}(i,t)$ and $p_{hop}(i,t)$ are computed using the actual finite-temperature structures, similar to the Method of Sussman et al. (2016); unlike that in Schoenholz et al. (2015, 2016), we do not use inherent structures (quenches to the local $T=0$ configuration). The training set is built by sampling uniformly at random $N_{train}=10^{4}$ spatio-temporal instants $(i,t)$ with $p_{hop}^{*}(i,t)>p_{c}$, and the same number from the other class. Finally, the third step of the procedure tests whether the algorithm classifies the data correctly by checking the performance on data that have not been used for training. For any choice of training set (WCA or LJ, $T_{train}$), one can apply the hyperplane $\vec{w}$ to all other test data (WCA or LJ, $T_{test}$) and measure the average classification accuracy (rate of correct answers) on balanced test sets ($50/50$ mix of $y_{GT}=-1$ and $y_{GT}=1$):

[TABLE]

where balance is obtained by downsampling until $\mathbb{P}(y_{GT}=-1)=1/2=\mathbb{P}(y_{GT}=1)$. Fig. 2 shows how this score depends on the choice of training set for various test sets. As expected the score decreases with increasing $T_{test}$: hotter is harder to predict because mobility is less correlated with structure. More interesting is that $T_{train}$ is not crucial: curves corresponding to different $T_{train}$ largely superimpose for the same training system. The choice of training system, however, is important: hyperplanes trained on a given system perform significantly better on that system. Overall, the high scores we achieve demonstrate a strong structure $-\vec{G}(i,t)-$ to mobility $-p_{hop}(i,t)-$ relationship encoded in the softness.

From now on, we shall set a unique hyperplane obtained by training at $T=0.43$, and using a balanced mix of WCA and LJ data. Other choices lead to similar (but slightly less good) results. We compute the local and instantaneous softness for each liquid using Eq. (3) and average it over all trajectories. Fig. 3(a) shows that the average softness $\left\langle S\right\rangle(1/T)$, a purely structural quantity, is quantitatively different for the two systems. However, when plotting $\tau_{\alpha}$ as a function of $\left\langle S\right\rangle$, the differences between the dynamics essentially disappear as shown in Fig. 3(b): not only does $\left\langle S\right\rangle$ amplify minute differences in the structure, it does so in such a way that the softness, on average, is an excellent structural indicator of the dynamics. Since the softness is linearly related to the structural descriptors, one simply has

[TABLE]

where $\vec{\left\langle G\right\rangle}\approx\left\langle g\right\rangle(r)$, in the limit of vanishingly small binning width. We thus conclude that the difference in dynamics between the LJ and the WCA systems is, to a large extent, encoded in $\left\langle g\right\rangle(r)$, but in a subtle way. The two populations of mobile and immobile particles are characterized by different instantaneous local $g(r)$, which our method allows to collapse into a scalar variable, the softness. Its average value is therefore related to the average mobility and, hence, to the relaxation time.

Fig. 3(c) shows that the width of the distribution of softness values is broad compared to the shift of $\langle S\rangle$ with temperature that is shown in Fig. 3(a). Fig. 3 compares the distribution $P(S)$ obtained for the LJ and WCA systems at low temperature ($T=0.43$), and that for the WCA system at a different temperature, $(T=0.14)$, such that they share the relaxation time $\tau_{\alpha}\simeq 10^{4}$. As expected, the distribution for the WCA system is shifted towards higher softness values as compared to the LJ ones at the same temperature. However the distributions have a strong overlap, suggesting that the softness field is highly fluctuating. One should thus not expect a deterministic relationship between the local, instantaneous value of the softness of a particle and the local mobility.

The distributions for the WCA and LJ systems, compared at different temperatures but at the same relaxation time, are quite similar and share a common average. Still, a systematic difference persists on the side of large softness, indicating again system-specific structure-mobility relationships. On the other hand, the negative-$S$ tails are very much alike. This last observation suggests that non-mobile, negative-$S$ particles play a system-independent role in controlling $\tau_{\alpha}$.

From a theoretical viewpoint, the weights $w$ in Eq. 3, or equivalently the softness field, provide a promising route towards a better understanding of the relation between the structure and dynamics of supercooled liquids. These weights encode the properties of the structural fluctuations, and more precisely the fluctuations of $g_{XY}(r)$. Moreover, Eq. (5) shows that the average softness $\left\langle S\right\rangle$ is a weighted integral of the average pair correlation function, where the weights are the components of $\vec{w}$, obtained from the training phase in the ML process, which appears directly related to the relaxation time (Fig. 3(b)). Fig. 4 displays these weights $w_{AX}(r_{\alpha})$, together with $\Delta\left\langle g_{AX}\right\rangle=\left\langle g_{AX}\right\rangle_{WCA}-\left\langle g_{AX}\right\rangle_{LJ}$, the minute differences between the average pair correlation functions for $X=A$ and $B$. High $w$ correspond to features most relevant to the dynamics. Apart from the general observation that the short distance structure is crucial in triggering dynamics, the weighting is non-trivial; it is not a mere amplification of modulations in the pair correlation function. The first minima and maxima of $\Delta\left\langle g_{AX}\right\rangle$ are amplified, but less than the tiny differences observed in the structure of the second shell. Conversely the second minimum in $\Delta\left\langle g_{AY}\right\rangle$ is strongly suppressed by the weighting.

From the machine-learning point of view, we proceed in an unusual way since we trained on a binary classification task (mobile, not mobile) but used the result as in a regression problem, i.e. focusing not only on the sign of $S$ but also on its value. In order to determine whether the magnitude of $S$ has physical meaning, we follow  Schoenholz et al. (2016, 2015) and compute $P(R|S)$ (left panel of Fig. 5, the probability of observing a local re-arrangement $(P(p_{hop}(i,t)>p_{hop}^{*}))$ for a given range of softness $\in[S,S+0.025]$, for different $1/T$. As expected intuitively, large $S$ lead to a tighter relationship with mobility, thus justifying the use of $S$ and not just its sign. As in Schoenholz et al. (2016, 2015) we find that $\log P(R|S)$ decreases linearly with $1/T$ and that

[TABLE]

Note that $\Delta E$ and $\Delta\Sigma$ depend far more strongly on $S$ for LJ than for WCA, as shown in Fig. 5(b,c). Following arguments of Ref. Schoenholz et al. (2016), the observed $S$ dependence leads to a significantly higher value of the onset temperature, $T_{0}$, for LJ than for WCA, consistent with prior calculations Berthier and Tarjus (2011a), as well as a stronger $T$-dependence as in Fig. 1(b).

In summary, we have shown that the softness field provides a quantitative, structural characterization of the dynamical slowdown in the mode-coupling regime, which is precise enough to distinguish liquids with very similar structures, namely LJ and WCA mixtures of equal composition. From the theoretical viewpoint many questions remain open with respect to the nature of the softness field. Perhaps the most important of these questions is whether the kind of fluctuations picked up by the softness are those envisioned in dynamical facilitation theory Keys et al. (2011), connected to mobile defects, or instead are related to the self-induced disorder that facilitates motion within thermodynamic theories such as random first-order transition theory Bhattacharyya et al. (2008); Franz et al. (2011); Biroli et al. (2018). This question can be answered by studying the dynamical evolution of the softness field, and by comparing it to the avalanche dynamics of supercooled liquids. We will address this question in future work.