Temperature ramping simulations

Ouro

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:

https://matgl.ai/tutorials%2FRelaxations%20and%20Simulations%20using%20the%20M3GNet%20Universal%20Potential.html

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.

Computation

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 $O(N^2)$ . This is a unique property of using an MLIP model, whereas with DFT traditionally is $O(N^3)$ .

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.

Dynamics environment

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

python

# NVE
from ase.md.verlet import VelocityVerlet
dyn = VelocityVerlet(atoms, timestep=1.0 * units.fs)

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

python

# NVT (Langevin)
from ase.md.langevin import Langevin
dyn = Langevin(atoms, timestep=1.0 * units.fs, temperature_K=300, friction=0.01/units.fs)

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

python

# NPT
from ase.md.npt import NPT
from ase.md.nptberendsen import NPTBerendsen

dyn = NPT(atoms, timestep=1.0 * units.fs, temperature_K=300, externalstress=0.0, ttime=25*units.fs, pfactor=75*units.fs)
# OR
dyn = NPTBerendsen(atoms=atoms, timestep=time_step_fs * units.fs, temperature_K=t_start, pressure_au=pressure * units.bar)

Supercell

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:

python

# Fetch a structure from Materials Project
with MPRester(MP_API_KEY) as mpr:
    struct = mpr.get_structure_by_material_id("mp-22851")

# Convert Pymatgen structure to ASE Atoms
atoms = AseAtomsAdaptor().get_atoms(struct)
atoms = atoms.repeat((3, 3, 3))  # 3x3x3 supercell

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.

Simulation parameters

Friction (relevant for NVT, Langevin)

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

Simulation design

Time step

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.