Molecular dynamics¶

This tutorial uses liquid water to illustrate molecular dynamics with a trained FlashCart potential in LAMMPS. It examines energy conservation and density in a periodic system of 256 water molecules. The simulations use Kokkos on one GPU and a time step of 0.5 fs. The two calculations probe different properties:

Energy conservation: 50 ps of NVT equilibration at 300 K, followed by 1.2 ns of NVE dynamics. The analysis reports the drift and fluctuation of total energy over the final 1 ns.
Density: 1.2 ns of NPT dynamics at 1 bar, at temperatures from 273.15 K to 373.15 K in 20 K steps. The analysis reports the mean density and its block-averaged standard error over the final 1 ns.

Run the simulations¶

The input files and scripts are in the repository’s examples/md/water directory. Use a LAMMPS build with ML-IAP, Python support, and CUDA-enabled Kokkos, available as lmp on PATH. See Build LAMMPS for build instructions.

The model must support hydrogen and oxygen and use eV for energy and Å for distance. It should describe water over the temperatures and densities studied.

From the repository root, select a GPU and run:

CUDA_VISIBLE_DEVICES=0 bash examples/md/water/run.sh models/my_run

Replace 0 with the GPU index for your run. This command runs the complete example: NVT equilibration followed by one NVE trajectory, then six NPT trajectories. The simulations run one after another on the GPU you specify.

The script creates the initial box and exports the model with flashcart-lammps models/my_run --compile-mode default if those files do not already exist. Existing model exports are reused. If the checkpoint or compilation setting has changed, run the export again before starting:

flashcart-lammps models/my_run --compile-mode default

With compilation enabled, the model is compiled when first evaluated in each LAMMPS process. See Compiled evaluation in LAMMPS interface.

Logs, trajectories, and saved structures are written to examples/md/water/results/my_run. The analysis prints its summary to the terminal. The directory name comes from the final component of the training directory path. Runs with the same name reuse this output directory and start the simulations again.

Energy conservation tests the force evaluation and time integration. It does not by itself establish the accuracy of the potential.

Build the box¶

The builder places 256 water molecules (768 atoms) in a periodic cube of side 19.706 Å, with random molecular orientations and a fixed seed. It writes a LAMMPS data file and an extended XYZ file. The initial box is unrelaxed. Both LAMMPS inputs minimize it before assigning velocities and starting dynamics.

LAMMPS inputs¶

in.nve

# Energy conservation: minimize, 50 ps NVT (300 K by default), then 1.2 ns NVE (0.5 fs steps).
#
#   lmp -k on g 1 -sf kk -pk kokkos newton on neigh half -in in.nve \
#       -var model /abs/path/flashcart-mliap.pt -var out results/<name>

variable model   index flashcart-mliap.pt
variable out     index results
variable temp    index 300.0        # NVT pre-equilibration
variable eqsteps index 100000       # 50 ps NVT before the NVE stage
variable steps   index 2400000      # 1.2 ns NVE at 0.5 fs
variable dt      index 0.0005       # ps
variable seed    index 4928459
variable tdamp   index 0.2          # ps
variable thermoevery index 200      # 0.1 ps
variable dumpevery   index 4000     # 2 ps

shell           mkdir -p ${out}

units           metal
atom_style      atomic
boundary        p p p
newton          on

read_data       water_256.data
pair_style      mliap unified ${model} 0
pair_coeff      * * H O

neighbor        1.0 bin
neigh_modify    every 1 delay 0 check yes

timestep        ${dt}
thermo          ${thermoevery}
thermo_style    custom step time temp pe ke etotal press
thermo_modify   format float %18.8f flush yes

minimize        1.0e-4 1.0e-6 300 3000
reset_timestep  0

velocity        all create ${temp} ${seed} mom yes rot yes dist gaussian
fix             eq all nvt temp ${temp} ${temp} ${tdamp}
run             ${eqsteps}
unfix           eq
write_data      ${out}/nve_start.data nocoeff

reset_timestep  0
fix             nve all nve
dump            traj all custom ${dumpevery} ${out}/nve.lammpstrj id type element x y z ix iy iz
dump_modify     traj sort id element H O
run             ${steps}

in.npt

# Density: minimize, then 1.2 ns NPT at 1 bar (300 K by default, set -var temp) with 0.5 fs steps.
#
#   lmp -k on g 1 -sf kk -pk kokkos newton on neigh half -in in.npt \
#       -var model /abs/path/flashcart-mliap.pt -var out results/<name> -var temp 320

variable model  index flashcart-mliap.pt
variable out    index results
variable temp   index 300.0        # K
variable press  index 1.0          # bar
variable steps  index 2400000      # 1.2 ns at 0.5 fs
variable dt     index 0.0005       # ps
variable seed   index 4928459
variable tdamp  index 0.2          # ps
variable pdamp  index 1.0          # ps
variable thermoevery index 200     # 0.1 ps
variable dumpevery   index 4000    # 2 ps

shell           mkdir -p ${out}

units           metal
atom_style      atomic
boundary        p p p
newton          on

read_data       water_256.data
pair_style      mliap unified ${model} 0
pair_coeff      * * H O

neighbor        1.0 bin
neigh_modify    every 1 delay 0 check yes

timestep        ${dt}
thermo          ${thermoevery}
thermo_style    custom step time temp pe ke etotal press vol density
thermo_modify   format float %18.8f flush yes

minimize        1.0e-4 1.0e-6 300 3000
reset_timestep  0

velocity        all create ${temp} ${seed} mom yes rot yes dist gaussian
fix             npt all npt temp ${temp} ${temp} ${tdamp} iso ${press} ${press} ${pdamp}

dump            traj all custom ${dumpevery} ${out}/npt_${temp}K.lammpstrj id type element x y z ix iy iz
dump_modify     traj sort id element H O

run             ${steps}

write_data      ${out}/npt_${temp}K_final.data nocoeff

Analysis¶

The analysis reports the slope of a linear fit to total energy as the energy drift in meV/atom/ns. Energy fluctuations are the standard deviation of total energy per atom in meV/atom, without subtracting the fitted drift.

For NPT trajectories, the analysis reports the mean density in g/cm³ and estimates its standard error from 50 ps block averages. The standard error quantifies uncertainty in the mean density due to the finite trajectory length. Compare densities with experimental values at matching temperatures and pressure. Simulation speed in ns/day is taken from the LAMMPS performance summary for the production run.

To repeat the analysis without rerunning the simulations:

python examples/md/water/analyze.py examples/md/water/results/my_run \
    --analysis-ps 1000 --block-ps 50

The default atom count is 768. If you change the box size, pass --n-atoms so that the energy drift and fluctuation are divided by the correct number of atoms.

Analysis script

"""Analyze energy conservation, density, and speed from the water simulations.

    python analyze.py results/<name>

Read log.nve and log.npt_<T>K from the run directory. Over the final
--analysis-ps of each trajectory, report energy drift in meV/atom/ns,
energy fluctuations in meV/atom, and density with a block-averaged standard
error. Simulation speed comes from the LAMMPS production-run summary.
"""

import argparse
import re
from pathlib import Path

import numpy as np

PERF_RE = re.compile(r"Performance: ([\d.eE+-]+) ns/day")
NVE_BLOCK, NPT_BLOCK = 2, 1


def read_log(path, block):
    """Thermo columns of one run block of a LAMMPS log and its ns/day (None if not reached)."""
    blocks = []  # [columns, rows, ns_per_day]
    for line in path.read_text().splitlines():
        parts = line.split()
        perf = PERF_RE.search(line)
        if parts and parts[0] == "Step":
            blocks.append([parts, [], None])
        elif blocks and perf:
            blocks[-1][2] = float(perf[1])
        elif blocks and len(parts) == len(blocks[-1][0]):
            try:
                blocks[-1][1].append([float(p) for p in parts])
            except ValueError:
                pass
    if len(blocks) <= block or len(blocks[block][1]) < 2:
        return None
    cols, rows, ns_per_day = blocks[block]
    table = np.array(rows)
    return {c: table[:, i] for i, c in enumerate(cols)}, ns_per_day


def block_sem(x, t, block_ps):
    """Standard error of the mean from block averages of length block_ps."""
    edges = np.arange(t[0], t[-1] + block_ps, block_ps)
    means = [x[(t >= a) & (t < b)].mean() for a, b in zip(edges[:-1], edges[1:]) if ((t >= a) & (t < b)).any()]
    return float(np.std(means, ddof=1) / np.sqrt(len(means))) if len(means) > 1 else float("nan")


def main():
    parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
    parser.add_argument("run_dir", type=Path)
    parser.add_argument("--n-atoms", type=int, default=768)
    parser.add_argument("--analysis-ps", type=float, default=1000.0, help="final window analyzed, ps")
    parser.add_argument("--block-ps", type=float, default=50.0, help="block length for standard errors, ps")
    args = parser.parse_args()

    log = args.run_dir / "log.nve"
    if log.exists() and (parsed := read_log(log, NVE_BLOCK)) is not None:
        thermo, ns_per_day = parsed
        t = thermo["Time"]
        w = t >= t[-1] - args.analysis_ps
        e = 1e3 * thermo["TotEng"][w] / args.n_atoms  # meV/atom
        slope = np.polyfit(t[w], e, 1)[0]  # meV/atom/ps
        print(f"NVE: drift {1e3 * slope:.4f} meV/atom/ns, std {e.std():.4f} meV/atom, "
              f"T {thermo['Temp'][w].mean():.1f} K, {ns_per_day} ns/day")

    for log in sorted(args.run_dir.glob("log.npt_*K")):
        parsed = read_log(log, NPT_BLOCK)
        if parsed is None:
            continue
        thermo, ns_per_day = parsed
        t = thermo["Time"]
        w = t >= t[-1] - args.analysis_ps
        rho = thermo["Density"][w]
        print(f"NPT {log.name[8:]}: density {rho.mean():.4f} ± {block_sem(rho, t[w], args.block_ps):.4f} g/cm^3, "
              f"T {thermo['Temp'][w].mean():.1f} K, {ns_per_day} ns/day")


if __name__ == "__main__":
    main()

For molecular dynamics with ASE, see Molecular dynamics in ASE calculator.