Molecular dynamics

Molecular dynamics moves atoms by integrating the forces a calculator returns. The API follows ASE, and timesteps are femtoseconds at the Python boundary.

Argon trajectory

Solid argon, 32 atoms, 1.0 fs timestep, 500 steps of NVE Verlet. Drag to orbit, and use the playback bar to watch the atoms oscillate about their lattice sites.

Ar 32 atoms, fcc 2x2x2 cubic cell. 101 frames sampled every 5 steps from a 500-step NVE run.

The script that produced it:

from atomli.build import bulk
from atomli.calculators import LennardJones
from atomli.md.velocitydistribution import MaxwellBoltzmannDistribution
from atomli.md.verlet import VelocityVerlet
from atomli.units import fs

atoms = bulk("Ar", "fcc", a=5.26, cubic=True).repeat((2, 2, 2))
atoms.calc = LennardJones(epsilon=0.0103, sigma=3.4, rc=5.0, shift=True)

MaxwellBoltzmannDistribution(atoms, temperature_K=60)

dyn = VelocityVerlet(atoms, timestep=1.0 * fs)
dyn.run(500)

Temperature equilibration

That run starts at 60 K and settles near 34.8 K, roughly half. That is not a bug and not a thermostat: MaxwellBoltzmannDistribution puts 60 K worth of kinetic energy in, and in a solid the atoms immediately trade half of it into potential energy climbing out of their lattice wells. Equipartition takes it back to about half the starting value.

So if you want to equilibrate at 300 K, either initialise at 600 K and let it fall, or use a thermostat. Initialising at 300 K and reporting 300 K is the single most common mistake in a first MD script.

Energy conservation

Velocity Verlet integrates NVE: no thermostat, so total energy is the quantity that must stay put while potential and kinetic energy trade against each other. Over the run above, total energy spans 0.017% of the mean kinetic energy (2.47e-5 eV).

Measure it yourself rather than assuming:

energies = []
dyn.attach(lambda: energies.append(
    atoms.get_potential_energy() + atoms.get_kinetic_energy()
), interval=5)
dyn.run(500)

drift = max(energies) - min(energies)

A drift that grows steadily rather than oscillating means the timestep is too large, the cutoff is wrong, or the calculator’s forces are not the gradient of its energy.

Timestep

Halving the timestep should sharply reduce drift. Here is what it actually does for this system, same seed, same 500 steps:

Timestep	Energy drift
0.5 fs	0.0052%
1.0 fs	0.0171%
2.0 fs	0.7158%
4.0 fs	2.4431%

Between 2.0 fs and 1.0 fs the drift falls by a factor of about 42. Below that the return diminishes and you are paying twice the compute for accuracy the rest of your model cannot justify. Argon is heavy and slow; a system containing hydrogen needs roughly 0.5 fs.

Cutoff radius

This one silently destroys a run. The cell above is 10.52 Å across, so the cutoff is set to 5 Å. Raising it past half the box makes an atom interact with two images of the same neighbour, the forces stop being the gradient of the energy, and “conserved” energy drifts by tens of percent while the trajectory still looks perfectly plausible on screen.

box = atoms.get_cell().lengths().min()
assert calc_rc < box / 2, "cutoff exceeds the minimum image convention"

Langevin thermostat

Langevin adds friction and a matching random force, so the system samples NVT at a temperature you choose instead of conserving energy.

from atomli.md.langevin import Langevin
from atomli.units import fs

dyn = Langevin(atoms, timestep=1.0 * fs, temperature_K=300, friction=0.01)
dyn.run(1000)

friction is a rate in inverse ASE time units. Large values thermostat hard and damp real dynamics; small values equilibrate slowly. For structural sampling 0.01–0.02 is a reasonable starting range. If you care about transport properties such as diffusion or vibrational spectra, equilibrate with Langevin and then measure in NVE with Verlet, because the thermostat’s random force corrupts the dynamics you are trying to observe.

Live visualization

view() returns a viewer you can stream frames into, so a relaxation or an MD run animates in the cell as it progresses rather than after it finishes.

from atomli.visualize import view

v = view(atoms)
dyn = VelocityVerlet(atoms, timestep=1.0 * fs)
dyn.attach(v.update, interval=5)
dyn.run(500)

v.save("run.gif", fps=20)

See Visualization for the full viewer API.