Sharing some things I'm learning as I work on temperature ramping simulations. The goal of these simulations is to learn how a material's lattice changes with temperature, as thermal expansion, decomposition, phase changes, etc.
I'm using ASE and different MLIP (CHGNet, M3GNet) models to model the molecular dynamics.
There are a few different pretrained models we can use immediately from https://matgl.ai/ (which has a nice interface to ASE):
CHGNet-MPtrj-2024.2.13-PES-11M
M3GNet-MP-2021.2.8-PES
See MatGL's documentation for some sample code:
I haven't seen many references to CHGNet in their docs yet. Instead, it all talks about M3GNet but I take it to be because CHGNet is just very new.
Even with a MLIP model taking the place of DFT, there's still a lot of computation happening. N atoms for M time steps adds up.
Thankfully, computation scales as . This is a unique property of using an MLIP model, whereas with DFT traditionally is .
We're absolutely going to need to be CUDA GPU accelerated. Despite cutting computation significantly thanks to MLIP, there is still a lot going on. Step sizes are so small, and you need to run sims slow enough for long enough to get accurate data.
Still, for a single material, we've cut weeks or months of CPU time into hours or days of GPU time.
We're using Langevin for general temperature ramping and NPT for lattice structural studies. Here's a quick look at some of the options we have.
NVE (Microcanonical): Models an isolated system with no energy exchange
Constant: Energy, Volume, Number of particles
Uses Velocity Verlet integration
Good for testing energy conservation
No temperature control
NVT (Canonical) via Langevin: Models a system in contact with a heat bath
Constant: Temperature, Volume, Number of particles
Adds friction and random forces
Good for temperature control and equilibration
NPT (Isothermal-Isobaric): Models a system at constant pressure and temperature. ASE does not recommend their NPT system, but instead their NPTBerendsen.
NPTBerendsen: The size of the unit cell is rescaled after each time step, so the pressure / stress approaches the desired pressure. Disadvantage: Fluctuations in both total energy and pressure are suppressed compared to the correct NPT ensemble. For large systems, this is not expected to be serious.
Constant: Temperature, Pressure, Number of particles
Allows volume fluctuations
Important for thermal expansion studies
Having a supercell is important. You really can't do any real simulation on just the unit cell, especially when it comes to temperature effects.
It's easy to tile your unit cell in ASE:
Minimum Size:
For basic structural transitions or thermal stability studies, use at least a 4×4×4 supercell.
This size ensures proper representation of periodic boundary conditions and allows for structural distortions.
Larger Sizes:
For defects, doping, or surface simulations, use 6×6×6 or larger. Larger systems reduce artificial periodic interactions and allow localized effects like defect clustering.
The friction coefficient in ASE is in frequency units (inverse time):
0.01 / units.fs = 0.01 fs⁻¹ = 10 ps⁻¹
This means the characteristic damping time is:
τ = 1/(friction) = 100 fs
So for crystal superconductor simulations, we can think about it in these timescales:
0.01 / units.fs: Damping time ~100 fs - gentle coupling
0.05 / units.fs: Damping time ~20 fs - medium coupling
0.1 / units.fs: Damping time ~10 fs - strong coupling
This aligns with typical phonon frequencies in crystals (0.1-10 THz or periods of 100-1000 fs). You want your damping time to be comparable to or longer than your fastest relevant dynamics.
For superconductors:
Use gentler coupling (0.01-0.02/fs) when studying phonon dynamics
Use stronger coupling (0.05-0.1/fs) for temperature equilibration
Avoid very strong coupling (>0.2/fs) as it could overdamp important lattice dynamics
Use a time step of 1–2 fs:
For systems with heavier atoms, this is sufficient to capture atomic vibrations without introducing instabilities.
Avoid larger time steps to prevent numerical inaccuracies, especially for high-temperature dynamics.
Cooling Ramps (Phase Transitions):
Ramp down temperature at a rate of 10–50 K/ps for observing phase transitions. Slower ramps (e.g., 1 K/ps) may give more accurate transition data but are computationally expensive.
Ensure the simulation time is long enough for the material to equilibrate at each temperature step.
Heating Ramps (Decomposition):
Ramp up temperature at a rate of 20–100 K/ps. Decomposition usually requires higher temperatures and longer simulation times.
Analyze bond-breaking events and structural distortions as the temperature increases.
Equilibration:
Equilibrate your system for 10–20 ps at your starting temperature (e.g., high temperature for phase transition studies).
Dynamic Simulations:
Run simulations for 50–100 ps for observing transitions or decomposition. Increase the time if you aim to study slower processes or diffusive phenomena.
Use longer simulations (~200 ps or more) for systems with defects or doping, where equilibration is slower.
Why is temperature fluctuating so much, why is it sometimes greater than the simulation environment?
Fixed with larger supercell, slower temperature ramping