6.4 Analyze trajectories with units and uncertainty
These are unexecuted teaching templates, with no ABACUS results. Baseline: ABACUS 3.9.0. Verify executable and external-tool versions.
6.4.1 Scientific question
Which conclusions can be justified by a finite atomic trajectory? A simulation contains many frames, but adjacent frames are correlated and may not constitute independent samples. This lesson turns the validated silicon AIMD workflow into a reproducible analysis pipeline. It separates coordinate handling, physical observables and statistical uncertainty. A short solid-state trajectory can diagnose stability and vibrations; it cannot automatically measure diffusion or a reaction rate. The displayed definitions and analysis template do not contain simulated data or claimed results.
6.4.2 Model and prerequisites
Archive the complete INPUT, initial and checkpoint STRU, KPT, executable version, data hashes, restart boundaries and raw OUT directory. Record physical units and output intervals from the files themselves. A parser must identify lattice vectors, species, atom ordering and time indices, not merely count lines. Test it on the initial frame by reconstructing known distances and on a selected later frame by comparison with the raw output. Keep a conversion log when using an external package. The teaching CSV is an intermediate format you must explicitly generate and inspect; the Python example does not parse ABACUS MD_dump directly.
6.4.3 Worked template and interpretation
The output delta requests coordinate-resolved force and velocity information at each saved step. For large production trajectories choose an interval based on storage and the shortest time scale of interest, then document it. The CSV example converts step indices to ps using the actual archived time step. The output stride belongs in the step difference and must not be multiplied in twice. Detect duplicate or reversed steps where restarted segments were concatenated. Verify per-atom normalization before comparing supercells of different sizes. Reject missing, nonfinite or differently defined energy columns. Keep the processed CSV and parser version alongside the raw input; a plot alone is not an adequate archive.
# External analysis of an explicitly prepared CSV, not raw MD_dump:
# Required columns: step, energy_eV_per_atom, temperature_K
import csv
import numpy as np
with open("analysis.csv", newline="") as f:
rows = list(csv.DictReader(f))
step = np.array([int(r["step"]) for r in rows])
dt_fs = 1.0 # replace with the actual archived INPUT value
time_ps = (step - step[0]) * dt_fs / 1000.0
energy = np.array([float(r["energy_eV_per_atom"]) for r in rows])
assert np.all(np.diff(step) > 0), "Duplicate/reversed steps"
assert np.isfinite(energy).all(), "Missing/nonfinite energy"
print("Energy range, eV/atom:", float(np.ptp(energy)))
6.4.4 Physics and units
Wrapped coordinates keep particles inside the simulation cell and are convenient for viewing structure. Unwrapped coordinates preserve boundary crossings and are needed for displacement statistics. For a fixed orthorhombic box, a nearest-image displacement can be tracked between sufficiently close frames; large jumps, changing cells and triclinic geometry require a general cell-aware method. Do not blindly apply Cartesian modulo operations to a skewed cell. Remove center-of-mass motion consistently if the observable requires it, and state that choice. The MSD has units of Ų when positions are in Å; a slope divided by six gives a three-dimensional isotropic diffusion coefficient in Ų/ps if lag time is in ps.
6.4.5 Outputs and analysis
Begin with time traces of energy, temperature, maximum force and selected bond distances. Mark equilibration and restart boundaries. For a solid, bounded displacements around lattice sites and persistent structure are expected diagnostics; a nonzero short-time MSD slope can reflect ballistic motion or vibrations rather than diffusion. A radial distribution function requires pair counting, shell-volume normalization, species populations and a valid density. For same-species pairs avoid self-counts and handle double counting consistently. Limit radii to a range supported by the cell geometry. Report bin width and compare a finer histogram before interpreting a narrow peak.
6.4.6 Convergence and acceptance
Estimate averages after a declared equilibration exclusion. Increase block lengths and inspect the stability of the uncertainty estimate as blocks become less correlated. Repeat independent initial conditions to test reproducibility. A large frame count from one short trajectory cannot replace independent physical time. Compare output cadence for fast observables and supercell size for long-wavelength or finite-size effects. If fitting diffusion, demonstrate a sustained linear diffusive regime over a defensible lag window, enough time origins and finite-size control. Silicon below melting is a useful example of why an absent diffusive regime should lead to withholding a diffusion coefficient rather than forcing a fit.
6.4.7 Pitfalls and exercises
Cross-correlations, thermostat-altered dynamics and restart discontinuities can invalidate naive error bars. Avoid presenting equilibration frames as equilibrium samples. Always label illustrative figures as schematics and distinguish proposed acceptance thresholds from measured values. Exercises: (1) construct a particle crossing a one-dimensional periodic boundary and show why wrapped MSD becomes wrong; then recover the displacement by unwrapping. (2) derive the conversion from Ų/ps to m²/s and show that the factor is 10⁻⁸. (3) explain why doubling saved-frame frequency at fixed duration does not halve statistical uncertainty. (4) choose a block-length test for a structural average and state which remaining finite-time limitation prevents a transport claim.