Models and Algorithms for the Next Generation of Glass Transition Studies

Successful computer studies of glass-forming materials need to overcome both the natural tendency to structural ordering and the dramatic increase of relaxation times at low temperatures. We present a comprehensive analysis of eleven glass-forming models to demonstrate that both challenges can be efficiently tackled using carefully designed models of size polydisperse supercooled liquids together with an efficient Monte Carlo algorithm where translational particle displacements are complemented by swaps of particle pairs. We study a broad range of size polydispersities, using both discrete and continuous mixtures, and we systematically investigate the role of particle softness, attractivity, and nonadditivity of the interactions. Each system is characterized by its robustness against structural ordering and by the efficiency of the swap Monte Carlo algorithm. We show that the combined optimization of the potential’s softness, polydispersity, and nonadditivity leads to novel computer models with excellent glass-forming ability. For such models, we achieve over 10 orders of magnitude gain in the equilibration time scale using the swap Monte Carlo algorithm, thus paving the way to computational studies of static and thermodynamic properties under experimental conditions. In addition, we provide microscopic insight into the performance of the swap algorithm, which should help optimize models and algorithms even further.

Popular Summary

If a liquid is cooled below its freezing temperature fast enough, it avoids freezing into a crystalline solid and becomes what is known as a supercooled liquid. As the temperature is decreased further, the motion of molecules within the liquid slows down enormously. This variation can be experimentally analyzed over about 12 orders of magnitude—roughly a factor of one trillion—down to the glass temperature, below which the liquid transforms to a glasslike solid. Computer simulations can track this change at the particle scale, but over short time scales only. Supercooled liquids in the laboratory are about one billion times slower than those currently simulated on a computer. In this work, we develop computational models and efficient simulation algorithms that close this colossal gap for the first time and even go beyond the dynamical regime currently accessible to experiments.

To accelerate thermalization, we perform simulations in which randomly chosen pairs of particles can exchange their positions. Such “swap moves” were introduced 30 years ago to study dense liquids, but previous attempts produced a speedup of at most 2 orders of magnitude or led to partly crystalline solids. We achieve speedups of more than 10 orders of magnitude without incurring crystallization. We systematically explore different classes of glass-forming liquids by varying distributions of particle size and interactions to simulate liquids at unprecedented degrees of supercooling, thus paving the way for the next generation of computer studies of the glass transition.

Our work shifts theoretical studies of disordered materials to an entirely new territory, which should close the gap between theory and experiments and get us close to solving the long-standing enigma of how the glass transition works.

Figure 1

Summary of the results obtained for the eleven models studied in this work. The top row shows the particle size distribution for each model; the next two rows specify the pair potential and its additivity. Below, we use a temperature axis rescaled by the location of the mode-coupling crossover, where blue points indicate equilibrium disordered fluid configurations and red diamonds indicate instability towards crystalline or demixed states. The bottom row indicates the estimated range of equilibrium relaxation times ${\tau }_{\alpha }$ that can be studied in stable equilibrium conditions for each model, using ${\tau }_0$ as the relaxation time at the onset temperature [53]. We have constructed several models that remain stable and can be equilibrated deep in the temperature regime $T/T_{\mathrm{MCT}}<1$, which conventional simulation studies are unable to penetrate, three of which allow us to reach temperatures below the experimental glass transition, conventionally defined as ${\tau }_{\alpha }/{\tau }_0=10^{12}$.

Figure 2

The relaxation time ${\tau }_{\alpha }$ as a function of the swap attempt probability $p$, normalized by ${\tau }_{\alpha }^{*}$, its value for standard Monte Carlo simulations when $p=0$. A broad minimum indicates that $p\approx 0.20$ optimizes the efficiency of the swap Monte Carlo algorithm. The inset shows the probability distribution of swap acceptance as a function of the diameter difference $\mathrm{\Delta }\sigma =|{\sigma }_1-{\sigma }_2|$ between the particles for which the swap move is attempted. The system is a soft repulsive system of nonadditive particles studied in Sec. 4c, with $\epsilon =0.2$ and $T=0.101$.

Figure 3

Potential energy per particle in simulations with different cooling rates $\gamma$ using both standard and swap Monte Carlo dynamics. Both dynamics produce consistent results at high temperatures, but the swap dynamics remains at equilibrium down to much lower temperatures than the standard one. The system is the soft repulsive system of nonadditive particles studied in Sec. 4c, with $\epsilon =0.2$.

Figure 4

Swap dynamics in a soft repulsive system of nonadditive particles studied in Sec. 4c with $\epsilon =0.2$. (a) Self-incoherent scattering $F_s(k,t)$ computed, respectively, on the first (solid lines) and second (dashed lines) halves of the simulation run. Results for standard dynamics without swap for the lowest temperature are shown with a dotted line. (b) Collective overlap correlation function $F_o(t)$. In both panels, temperatures are $T=0.25$, 0.175, 0.125, 0.092, 0.075, 0.065, 0.062, 0.058, and 0.0555, from left to right. Swap Monte Carlo simulations fully decorrelate both single-particle and collective density fluctuations in a regime where standard Monte Carlo simulations may be fully arrested.

Figure 5

Relaxation times ${\tau }_{\alpha }$ of standard (black empty points) and swap (red solid squares) simulations for the binary mixture of soft repulsive spheres studied in Sec. 3a. The speedup offered by the swap moves is obvious, but the system is unstable below $T=0.202$ where it crystallizes, and temperatures below the mode-coupling temperature $T_{\mathrm{MCT}}=0.199$ cannot be studied. The inset shows a time series of the potential energy for standard (black) and swap (red) simulations at $T=0.2$, and we see that crystallization is easily observed when swap moves are introduced.

Figure 6

Relaxation times ${\tau }_{\alpha }$ of standard (black empty points) and swap (red solid squares) simulations for the ternary mixture of soft repulsive spheres studied in Sec. 3b. The speedup offered by the swap moves is obvious, but the system is unstable below $T=0.26\approx 0.9T_{\mathrm{MCT}}$ where it demixes and crystallizes. Disconnected squares are a rough estimate of ${\tau }_{\alpha }$ obtained using short simulations in the unstable region. Stable and equilibrated states can be accessed down to $T\approx 0.26<T_{\mathrm{MCT}}=0.288$, extending the dynamic range by about 2 orders of magnitude as compared to standard simulations.

Figure 7

Partial structure factor $S_{CC}(q)$ for small particles in the ternary mixture of Sec. 3b. It is featureless at high enough temperatures, $T=0.350$, and displays strong composition fluctuations at low $k$ in the fluid at $T=0.267$, which may eventually lead to a demixed state at long times. For $T\le 0.26$, the system is demixed, as shown for $T=0.256$. The inset shows a representative snapshot of a demixed and partially crystallized system at $T=0.256$.

Figure 8

Relaxation times for continuously polydisperse systems of repulsive soft spheres with different softness exponents $n=8$, 12, 18, and 24. Temperatures are scaled by $T_{\mathrm{MCT}}$ to allow direct comparison between models, with $T_{\mathrm{MCT}}=0.143$, 0.267, 0.468, and 0.662, respectively. Open symbols represent the standard Monte Carlo dynamics, and closed symbols show the swap algorithm, for which unconnected symbols represent structurally unstable state points where only a rough estimate of ${\tau }_{\alpha }$ is obtained in short simulations. A larger $n$ yields better efficiency and better structural stability.

Figure 9

Relaxation times for systems with continuous polydispersity, $n=12$, and different nonadditivity parameters $\epsilon$. Temperatures are scaled by $T_{\mathrm{MCT}}$ to allow direct comparison between models, with $T_{\mathrm{MCT}}=0.267$, 0.176, 0.104, and 0.0534, respectively. Open symbols represent the standard Monte Carlo dynamics, and closed symbols show the swap algorithm, for which unconnected symbols represent structurally unstable state points where only a rough estimate of ${\tau }_{\alpha }$ is obtained in short simulations. A well-chosen amount of nonadditivity, $\epsilon \approx 0.1–0.2$, considerably improves the efficiency of thermalization and the structural stability.

Figure 10

(a) Relaxation times for the nonadditive model with $n=12$ and $\epsilon =0.2$ for standard and swap Monte Carlo dynamics. The standard dynamics is fitted with the VFT, parabolic, and Arrhenius laws, as shown with lines, which are used to estimate the location of the experimental glass temperature $T_g$, as shown with vertical dashed lines. For this system, the swap dynamics is able to provide stable and thermalized configurations at temperatures below $T_g$. (b) Relaxation times obtained from standard (open symbols) and swap (filled symbols) dynamics for various size polydisperse models of various softness ($n$) and nonadditivity ($\epsilon$) are shown in an Arrhenius form with rescaled temperature $T_g/T$, where $T_g$ is estimated as in (a). For all models, the thermalization speedup near $T_g$ is of about 10 orders of magnitude, some models being structurally stable down to temperatures below $T_g$.

Figure 11

Time series of individual displacement $\mathrm{\Delta }r(t)$ and diameter value $\sigma (t)$ for two tagged particles in the nonadditive model of Sec. 4c with $\epsilon =0.2$ at $T=0.0555$. Intermittent diffusion in real and diameter space is observed, with strong correlations between $\mathrm{\Delta }r(t)$ and $\sigma (t)$ highlighted with dashed lines, but we also observe many events in one observable that have no counterpart in the other, indicating that the correlation between the two observables is nonlocal.

Figure 12

Nonadditive model of Sec. 4c with $\epsilon =0.2$. (a) Time autocorrelation of particle diameters $C_{\sigma }(t)$ measured during the swap dynamics for temperatures, as in Fig. 4. (b) Relaxation times ${\tau }_{\alpha }$, ${\tau }_{\sigma }$, and ${\tau }_o$ as a function of temperature, with ${\tau }_{\sigma }$ and ${\tau }_o$ rescaled to coincide with ${\tau }_{\alpha }$ at $T=0.175$ (shown with the horizontal bar). The three time scales obviously have the same temperature dependence.

Figure 13

Snapshots of the 10% of particles with (a) the largest displacements and (b) the largest diameter change, computed between times $t=0$ and $t={\tau }_{\alpha }/2$ for $T=0.0555$. There is a clear correlation between the spatial regions where dynamics in real and diameter spaces are fast, but the correlation is weak at the single-particle level.

Figure 14

Temperature evolution of dynamic susceptibilities (right axis) and relaxation times (left axis). The vertical dashed line is $T_{\mathrm{MCT}}=0.104$. Whereas ${\chi }_4^d(t)$ grows together with ${\tau }_{\alpha }$ in standard simulations, ${\chi }_4^d(t)$ and ${\chi }_4^{\sigma }(t)$ behave similarly and have a very different temperature dependence which instead mirrors the evolution of the swap relaxation time.

Figure 15

Relaxation times for repulsive and Lennard-Jones potentials with a hybrid particle size distribution. Temperatures are scaled by $T_{\mathrm{MCT}}$ to allow direct comparison between models, with $T_{\mathrm{MCT}}=0.680$ and 0.543 for repulsive and LJ potentials, respectively. Open symbols represent the standard Monte Carlo dynamics, closed symbols the swap algorithm, for which unconnected symbols represent structurally unstable state points where only a rough estimate of ${\tau }_{\alpha }$ is obtained in short simulations.