Models and algorithms for the next generation of glass transition studies

Andrea Ninarello; Ludovic Berthier; Daniele Coslovich

Source-linked AI summary

Models and algorithms for the next generation of glass transition studies

Andrea Ninarello, Ludovic Berthier, Daniele Coslovich

arXiv:1704.08864v1Apr 28, 2017cond-mat.stat-mech

TL;DR

Glass-transition simulations must address both structural ordering and sharply increasing relaxation times at low temperatures. This work systematically evaluates model design and swap Monte Carlo to improve glass-forming stability and low-temperature equilibration.

Problem
Simulations provide particle-level observables but cover far less of the glassy slowdown than molecular-liquid experiments, motivating more efficient equilibrium-sampling strategies.
Method
The study systematically varies particle-size distributions and pair interactions while combining ordinary particle displacements with swap Monte Carlo moves.
Results
Particular combinations of model parameters achieve both excellent glass-forming ability and a dramatic reduction in the computer time needed to obtain thermalized low-temperature configurations.
Takeaways & Limitations
Optimizing model parameters together with swap Monte Carlo enables simulations of thermalized configurations in regimes otherwise inaccessible to ordinary dynamics.

Abstract

from arXiv · show

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 non-additivity 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 optimisation of the potential's softness, polydispersity and non-additivity leads to novel computer models with excellent glass-forming ability. For such models, we achieve over ten orders of magnitude gain in the equilibration timescale 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 insights into the performance of the swap algorithm which should help optimizing models and algorithms even further.

I. INTRODUCTION

Glass-transition simulations face a mismatch between particle-level spatial detail and experimentally relevant timescales, compounded by crystallization and demixing. This work addresses both constraints by combining optimized polydisperse models with swap Monte Carlo.

Motivation: Simulations provide particle-level resolution, but typically access only 4–5 decades of slowdown versus 12–13 decades in molecular-liquid experiments.The resulting gap is about eight orders of magnitude between simulation and molecular-glass regimes.
Motivation: Alternative algorithms are needed because brute-force computing improvements have extended the simulation time window by only about 3 orders of magnitude over 30 years.Collective moves and replica exchange provide gains of at most factors of 40 and about two orders of magnitude, respectively, in the cited dense-fluid studies.
Related simulation strategies: Replica exchange also scales poorly with particle number, losing most efficiency for systems containing thousands of particles.This limits its usefulness for bulk glass-transition studies.
Swap Monte Carlo: Swap Monte Carlo has reported a factor-180 equilibration speedup in a binary mixture, but that model crystallizes and has poor glass-forming ability.A later ternary model improved stability, yet its accessible dynamical window increased by only about 2–3 orders of magnitude after detailed analysis.
Paper approach: The paper systematically varies size distributions, softness, attractivity, and additivity across eleven models to optimize both stability against ordering and thermalization efficiency.The study reports particular parameter combinations with excellent glass-forming ability and dramatically reduced equilibration times.
Model design: Appropriate polydispersity must prevent crystallization and phase separation while also maintaining high particle-swap acceptance rates for efficient equilibrium sampling.The authors introduce discrete, continuous, and hybrid distributions; hybrid distributions combine discrete-mixture glass formation with continuous-distribution swap efficiency.
Model design: Non-additive interactions improve structural stability by allowing greater overlap between small and large particles while leaving identical-diameter interactions unchanged.The paper treats additivity and non-additivity as distinct interaction classes.

B. Physical observables

The paper defines structural, dynamical, and equilibration observables to assess stability and thermalization in supercooled models. It also compares relaxation-time fitting forms and notes that swap dynamics is algorithmically useful but microscopically nonphysical.

Observables: The study monitors structure and dynamics to characterize model stability, equilibration, and thermalization under supercooled conditions.Observables include structure factors, bond-orientational order, potential energy, and self and collective correlation functions.
Structural observables: The structure factor probes long-wavelength density fluctuations, demixing, and instability in size-disperse systems.Continuously polydisperse particles are grouped into comparable-size families for partial structure factors.
Structural observables: The six-fold bond-orientational order parameter and potential energy are inspected across time and temperature to detect crystalline ordering.Neighbors are selected within a radius tied to the minimum of the rescaled radial distribution function.
Interpretation: Swap trajectories are not physical, although their self and collective time correlations still quantify particle diffusion and density-fluctuation relaxation timescales.Swap moves exchange particle diameters rather than positions.
Dynamical observables: The self-relaxation time τα is defined by Fs(k, τα) = e^-1, with k at the first peak of the total structure factor.Because particle diameters change during swaps, the self-correlation sum runs over all particles.
Relaxation-time fits: VFT, parabolic, and Arrhenius fits bracket low-temperature relaxation behavior: VFT likely overestimates growth, while Arrhenius underestimates it.The combined fits provide a confidence interval for extrapolating accessible relaxation times.
Dynamical observables: The collective relaxation time τo is defined by Fo(τo) = e^-1 using a cutoff distance a = 0.3.The overlap function provides information similar to the coherent intermediate scattering function while reducing statistical fluctuations.

C. Efficiency of swap moves

Swap efficiency depends non-monotonically on the attempt probability, so the algorithm must balance diameter exchanges with translational motion. In the tested model, a low swap cost accompanies orders-of-magnitude faster structural relaxation, while efficiency remains insensitive to system size.

Optimization: Swap efficiency has two limiting cases: p = 0 gives standard Monte Carlo, whereas p = 1 prevents structural relaxation because positions never change.The useful regime therefore requires both displacement and swap moves.
Optimization: p ≈0.20 optimizes swap efficiency in the tested soft repulsive non-additive system at T = 0.101.The relaxation-time curve has a broad minimum when normalized by its p = 0 value.
Computational cost: One swap sweep with p = 0.2 takes only 20% longer in CPU time than a standard sweep, despite the much larger relaxation-time gain.A swap attempt evaluates local energy changes for two particles rather than one.
System-size scaling: Swap efficiency and implementation are insensitive to particle number N, unlike replica exchange, which scales poorly with N.This makes swap Monte Carlo suitable for larger systems typically used in bulk glass-transition studies.
Acceptance: Swap acceptance decreases as particle diameter differences increase, making continuously polydisperse systems especially favorable for efficient exchanges.Small diameter differences allow particles to fit more readily into one another’s local environments.

D. Equilibration and metastability

Glass simulations must maintain equilibrium long enough to sample observables while avoiding crystallization, demixing, or other ordering. The paper applies strict equilibration and instability tests and shows that swap dynamics extends reliable equilibrium sampling far beyond standard Monte Carlo.

Motivation: Glass-forming simulations face a compromise between sufficiently long equilibrium sampling and avoiding crystallization or structural ordering.Binary and weakly polydisperse models can crystallize over long simulation times.
Structural stability: Enhanced sampling can expose instabilities hidden by conventional simulations, including crystallization in repulsive mixtures and fractionated crystals at high polydispersity.Previously acceptable glass models may fail stability criteria once swap moves probe lower temperatures.
Protocol: All models are compared using an identical protocol combining standard-Monte-Carlo observables with swap-based equilibration checks.The protocol extracts energy, structure factors, and relaxation times before assessing target-temperature equilibration.
Protocol: Equilibrated swap simulations are extended for 200τα to measure static and dynamic properties over a broad time window.Thermalization therefore requires a simulation window two orders of magnitude longer than the structural relaxation time.
Metastability: A state point is classified as unstable if at least one of five independent simulations shows ordering or other instability within 200τα.The criterion is intentionally strict because ordering fluctuations can interfere with metastable-fluid physics.
Results: Swap dynamics remains at equilibrium to much lower temperatures than standard dynamics under cooling-rate comparisons.Both methods agree at high temperatures, while swap curves retain equilibrium behavior deeper into cooling.
Results: Swap simulations show consistent first- and second-half scattering functions, indicating no aging, while standard dynamics reaches a long-lived plateau at the lowest temperature.The collective overlap also decorrelates in the swap simulations.

A. Binary mixtures

The historical 50:50 binary soft-sphere mixture benefits substantially from swaps but remains structurally unstable near its mode-coupling crossover. Its low swap acceptance creates a stability–efficiency trade-off that limits further improvement within simple binary mixtures.

Swap performance: Swap moves accelerate relaxation by about 2 orders of magnitude at the lowest equilibrated temperature despite acceptance rates of order a ∼10^-2.The speedup is measured over the temperature range that satisfies the equilibration criteria.
Stability: The system crystallizes at the lowest studied temperature and becomes unstable below T = 0.202, only marginally above TMCT ≈0.199.Swap moves improve efficiency but do not substantially extend the accessible temperature range.
Stability: At T = 0.2, rapid crystallization is observed when swap dynamics is employed, making the historical model too poor a glass-former for novel low-temperature regimes.Detecting and filtering crystallized configurations would require complex strategies near the mode-coupling crossover.
Limitations: Within simple binary mixtures, increasing the size ratio could improve stability but would make the already low swap acceptance vanishingly small.This stability–efficiency trade-off leaves little room for drastic algorithmic improvement.

B. A ternary mixture

The ternary mixture improves access below the mode-coupling crossover, but detailed relaxation and structural tests reveal low-temperature demixing, crystallization, and failed thermalization. Swap moves substantially accelerate dynamics, yet the stable dynamic-range extension is about two decades rather than the previously claimed ten.

Model: The ternary model uses soft spheres with size ratio σA/σC = 1.25 and compositions xA = 0.55, xB = 0.30, and xc = 0.15.These compositions give approximately equal volume fractions and δ ≈17% polydispersity.
Dynamics: Stable and equilibrated states reach T ≈0.26 below TMCT = 0.288, extending the accessible dynamic regime by about two decades versus standard simulations.The model is nevertheless unstable below T = 0.26, where it demixes and crystallizes.
Thermalization: Swap simulations fail to thermalize the model for T ≤0.24, and energy-distribution reweighting can miss this failure.Structural relaxation and relaxation dynamics provide a more discriminating thermalization test than global static observables alone.
Assessment: Compared with the binary mixture, the ternary mixture preserves swap efficiency and extends thermalization below TMCT, but less dramatically than previously claimed.The authors’ detailed characterization rejects the earlier ten-decade efficiency claim for this model.
Stability: For T ≤0.26, the fluid becomes structurally unstable: composition fluctuations strengthen at low wavevector, followed by demixing and partial crystallization.Longer simulations also trigger demixing at nominally stable T = 0.267 state points.
Dynamics: At T = 0.26, the estimated standard relaxation time is τα ≈10^9, whereas swap simulations give τα ≈5 · 10^5.The corresponding thermalization speedup is about three orders of magnitude.

C. Five-component mixtures

Five-component mixtures increase swap acceptance through reduced size ratios, but both studied models strongly demix during swap simulations and cannot be equilibrated well below TMCT. The authors therefore leave broader discrete-model optimization open.

Model: Two five-component models use linearly spaced diameters with polydispersities δ = 16% and δ = 23%.Their diameter ranges are σmin = 0.847 to σmax = 1.333 and σmin = 0.826 to σmax = 1.771, respectively.
Swap efficiency: Swap acceptance increases to approximately 10%–20% because individual component size ratios are reduced.The acceptance range depends on temperature.
Stability: Both five-component models strongly demix during swap Monte Carlo, making equilibration well below TMCT impossible.Thus, higher component count alone does not resolve the stability–efficiency trade-off.
Scope: The discrete-model parameter space remains largely unexplored as the number of components increases.The authors identify systematic exploration of additional parameter combinations as necessary.

IV. CONTINUOUSLY POLYDISPERSE SYSTEMS

Continuous polydispersity is introduced to improve swap acceptance while suppressing crystallization and fractionation. Within this family, steeper repulsions provide better efficiency and structural stability, with n = 24 stable and efficiently thermalized down to about 0.6TMCT.

Motivation: Reducing diameter differences improves swap acceptance but lowers polydispersity, increasing crystallization risk; multicomponent mixtures also tend to demix at low temperature.These competing effects motivate continuous particle-size distributions.
Design principle: Continuous size distributions provide many pairs of particles with similar diameters, enabling successive accepted swaps that facilitate thermalization.Higher polydispersity can also stabilize the liquid against crystallization and fractionation.
Model: The studied systems contain N = 1500 particles at ρ = 1, with particle-size polydispersity δ varied through σmax/σmin.The interactions use a cutoff distance rcut = 1.25σij.
Influence of the particle softness: Swap moves significantly accelerate thermalization for all softness exponents n = 8, 12, 18, and 24.Temperatures are scaled by each model’s TMCT for direct comparison.
Influence of the particle softness: Larger n values permit equilibration at increasingly lower temperatures relative to TMCT and yield better structural stability.The n = 24 system remains stable and efficiently thermalized down to T ≈0.6TMCT.
Stability: Even continuously polydisperse models can demix at sufficiently low temperatures during swap simulations.Unconnected symbols indicate structurally unstable state points where relaxation times are only roughly estimated from short simulations.

C. Non-additive interactions

Non-additivity suppresses phase separation and can improve both structural stability and swap-Monte-Carlo efficiency, with an optimum near ǫ = 0.1. The resulting models remain equilibrated below conventional glass-transition temperature estimates, although timescale gains become extrapolation-sensitive.

Swap efficiency: Non-additivity improves swap efficiency because lower temperatures relative to TMCT can be studied, although its detailed physical basis remains unclear.The improvement is reported relative to the additive model with the same soft-repulsion exponent n = 12.
Structural stability: ǫ = 0.1 uniquely remained stable and equilibrated across the entire accessible temperature range among the eleven studied models.At ǫ = 0.2 and 0.3, crystallization events were frequently observed at low temperature.
Structural stability: Non-additivity strongly suppresses phase separation for ǫ = 0.1 and 0.2, while larger values can permit crystallization.The authors associate the optimal parameter with competition between demixing at ǫ = 0 and crystallization at large ǫ.
Accessible temperatures: T ≈ 0.5TMCT became accessible after optimizing additivity and the pair potential, compared with T ≈TMCT for the corresponding binary mixture.This decrease in the thermalization and stability limits is described as a major methodological improvement.
Experimental timescales: 10 (±1) orders of magnitude is the estimated thermalization speedup at Tg for the illustrated non-additive model.The maximal speedup is larger because temperatures below Tg can also be thermalized, but estimating that gain depends sensitively on extrapolation.
Experimental timescales: Selected models can access a temperature regime unavailable to experiments, opening a previously unexplored observational window on glass physics.The paper identifies efficient glass-forming models and algorithms as its main achievement.

V. MICROSCOPIC INSIGHTS INTO THE SWAP DYNAMICS

Swap dynamics accelerates relaxation through efficient motion in particle-diameter space rather than by simply allowing particles to escape positional cages. Diameter, positional, and collective-density relaxation share the same temperature-dependent timescale, while local correlations remain imperfect.

Mechanism: Swap moves do not explain their speedup simply by letting particles escape positional cages, because exchanged particles still occupy caged positions.The authors instead analyze swaps as diameter exchanges without positional changes.
Mechanism: Diameter fluctuations wander through continuous size space and provide the efficient dynamics that drives structural relaxation in position space.At low temperature, diameter relaxation appears as intermittent jumps after an initial caged regime.
Collectivity: Diameter and positional jumps sometimes coincide but also occur independently, indicating collective relaxation rather than a sufficient single-particle trigger.Changes in a particle’s neighborhood may trigger displacement even without a matching diameter change.
Timescales: τα, τσ, and τo have the same temperature dependence, showing that positional, diameter, and collective-density relaxation occur on a common timescale.Here τσ is defined by Cσ(τσ) = e^-1.
Collectivity: The spatial regions of fastest positional and diameter dynamics correlate clearly, despite weak correlation at the individual-particle level.The comparison uses the 10% of particles with the largest displacements in each space.

B. Spatially heterogeneous dynamics

Swap dynamics preserves equilibrium static properties while changing dynamic correlations. Near TMCT, its weak spatial dynamic correlations and the matching temperature dependence of diameter and positional fluctuations indicate that the acceleration is largely dynamic and cooperative.

Cooperativity: The correspondence of diameter and positional timescales, combined with weak single-particle correlation, supports a cooperative mechanism for diffusion in both spaces.The proposed spatial picture is that efficient diameter dynamics facilitates diffusion in corresponding regions.
Dynamic heterogeneity: χd4(t) and χσ4(t) behave similarly and follow the temperature dependence of the swap relaxation time.The susceptibilities peak around τα and τσ, respectively, and quantify spatial correlations of local observables.
Dynamic heterogeneity: Swap dynamics near TMCT displays essentially no spatial dynamic correlations, unlike standard dynamics whose susceptibility grows toward χd*4 ≈ 12.Because swaps preserve equilibrium static properties, the contrast identifies the standard-dynamics correlations as primarily dynamic rather than structural.
Dynamic heterogeneity: χd*4 and χσ*4 are quantitatively close under swap dynamics, confirming strong correlation between real-space and diameter-space fluctuations.Their shared temperature evolution follows the swap relaxation time.

VI. IDEAS FOR THE FUTURE DESIGN OF GLASS-FORMING MODELS

The paper uses its microscopic understanding of swap dynamics to propose new directions for designing glass-forming models in which swap Monte Carlo remains efficient.

Future design: The microscopic analysis motivates new design strategies for glass-forming models optimized for efficient swap Monte Carlo.This section builds directly on the preceding understanding of the swap mechanism.

A. “Hybrid” models for binary mixtures

Hybrid size distributions combine binary-mixture-like structural stability with swap-compatible intermediate particle sizes. The resulting model equilibrates efficiently below T_MCT, remains resistant to ordering down to about 0.75T_MCT, but develops low-temperature phase-separation instabilities.

Hybrid design: Hybrid distributions retain two main particle-size peaks while adding an interpolating third species to enable swap-compatible exchanges.Particles can move between the two main species through many intermediate swaps with high acceptance rates.
Hybrid design: A size ratio of σB/σA = 1.6 separates the main peaks, while σC = (σA + σB)/2 connects them continuously.The reported composition is xA = 0.33, xB = 0.34, and xC = 0.33, with 20% polydispersity.
Performance: For the hybrid repulsive model, equilibration is easily attainable below T_MCT and resistance to ordering persists to approximately 0.75T_MCT.The study uses N = 1000 particles at density ρ = 1.3.
Performance: At lower temperatures, the hybrid model develops instabilities associated with a tendency toward phase separation.The authors suggest studying variants with a smaller concentration of the interpolating species.
Lennard-Jones variant: Adding Lennard-Jones attraction changes the dynamics little beyond rescaling temperature, shifting T_MCT from 0.680 to 0.543.The comparison uses standard and swap dynamics for corresponding repulsive and attractive systems.
Overall conclusion: Across the studied models, carefully chosen particle-size distributions and interactions support efficient swap Monte Carlo simulation while preserving strong glass-forming ability.The paper considers variations in pair interactions, size distributions, and non-additivity.