6.2 NVE trajectories and energy drift

These are unexecuted teaching templates, with no ABACUS results. Baseline: ABACUS 3.9.0. Verify executable and external-tool versions.

Open figure at full size

6.2.1 Scientific question

Does the numerical trajectory approximately conserve the energy of an isolated model system? NVE fixes particle count, volume and total energy. Potential and kinetic energy exchange as atoms vibrate; neither part must remain individually constant. The diagnostic is the sum. This lesson uses a relaxed silicon supercell with all relevant atoms mobile, fixed cell vectors and the electronic model established earlier. It is a short integrator diagnostic rather than a production thermal-property calculation. No trajectory or drift value has been computed for this course.

6.2.2 Model and prerequisites

Start from a structure whose maximum residual force is already acceptable and whose electronic force convergence has been tested. Copy its complete INPUT, STRU and KPT into a fresh directory. Retain the pseudopotential/orbital pairing, atom ordering and basis settings. A minimal primitive cell forces every periodic image to move coherently, so a supercell is preferable for thermal diagnostics; test its size for eventual observables. Initial velocities matter as much as coordinates. Record the seed and initialization policy, and inspect the actual initial temperature and center-of-mass motion rather than assuming the requested temperature uniquely defines a state.

6.2.3 Worked template and interpretation

Replace the listed keys rather than adding a second copy. md_dt is in femtoseconds; md_tfirst is in kelvin. The proposed 200 steps at 1 fs cover only 0.2 ps. They are a commissioning test, not proof of equilibrium or adequate sampling. Do not change the cell during NVE. Dumping every step preserves enough time resolution to inspect numerical oscillations. The restart interval is distinct from the trajectory interval: an output frame is not automatically a restart state. Use the restart lesson to preserve velocities, electronic initialization and the actual program checkpoint together. No scheduler submission is required by this teaching text.

# Delta from a converged, relaxed silicon supercell INPUT:
calculation md
suffix si_nve
cal_force 1
symmetry 0
scf_thr 1e-10
md_type nve
md_nstep 200
md_dt 1.0
md_tfirst 300
md_seed 17
md_dumpfreq 1
md_restartfreq 10
dump_force 1
dump_vel 1

6.2.4 Physics and units

The total-energy definition should match the electronic force model. For a gapped silicon baseline, use a consistent fixed-occupation electronic treatment. Metals with electronic smearing require careful distinction between reported free-energy-like quantities and extrapolated zero-temperature energies; do not combine a force-consistent potential with an unrelated energy column. Express energy in eV per atom and time in ps when reporting drift, giving units of eV atom⁻¹ ps⁻¹. ABACUS input time units and an analysis plot's time units can differ by a factor of one thousand. The equation is a definition, not an assertion that a finite-step algorithm conserves it exactly.

\[ E_{\mathrm{tot}}(t)=E_{\mathrm{BO}}(\mathbf R(t))+\sum_i\frac{m_i|\mathbf v_i(t)|^2}{2},\qquad b=\frac{d(E_{\mathrm{tot}}/N)}{dt} \]

6.2.5 Outputs and analysis

Inspect running_md.log and MD_dump in OUT.si_nve, including the printed units and SCF status for each step. Read coordinates, velocities and forces with a version-compatible parser. Verify the number of frames, time origin and atom count before making a plot. Plot potential, kinetic and their consistent total against physical time. Fit a drift slope over a declared window and report its uncertainty; also report the bounded oscillation amplitude. A slope close to zero over a very short window does not prove long-time stability. Check for abrupt energy jumps aligned with an SCF failure or a restart boundary.

6.2.6 Convergence and acceptance

Repeat the same initial state using 0.5 fs and 1 fs with identical physical duration, doubling the number of steps in the smaller-step branch. Compare drift and oscillation amplitude. Then tighten SCF at fixed step size to separate integration error from force noise. Repeat selected forces with richer orbitals or higher cutoff. Long trajectories diverge microscopically because of nonlinear dynamics, so compare statistical and conservation diagnostics rather than requiring identical positions at late times. Choose acceptance tolerances appropriate to the smallest energy or dynamical signal of interest. Keep rejected runs and their failure reasons in the evidence table.

6.2.7 Pitfalls and exercises

A thermostat can conceal energy drift by exchanging energy with the system; it cannot validate an NVE integrator. A temperature rise can reflect conversion of strain energy from an unrelaxed starting structure rather than simple numerical heating. Frozen atoms alter the mobile degrees of freedom used to infer temperature. Exercises: (1) convert a trial slope from eV per cell per fs to eV per atom per ps without inserting a fabricated measured value. (2) Explain why halving the time step requires doubling steps for a fair-duration test. (3) Design a crossed time-step/SCF test and explain which pattern implicates force noise. Archive the plot recipe and raw columns so the conclusion can be reproduced.