Molecular Dynamics: Statistics and Validation

A molecular-dynamics trajectory becomes a scientific result only through analysis. This chapter covers the path from trajectory to defensible claim: designing the run around an observable, demonstrating equilibration, measuring statistical convergence and effective sample size, validating the force model along the visited configurations, and reporting with language proportional to the evidence.

The equations of motion, thermostats, barostats, and the exact Tako runtime behavior are developed in Molecular Dynamics; current implementation defects are registered in Known Limitations.

Designing a study from a scientific question

The most reliable MD protocol begins with an observable and a claim, not with a thermostat menu. “Estimate the room-temperature distribution of a torsion within one conformational basin,” “test whether a proposed crystal remains mechanically intact for the accessible simulation time,” and “compare short-time vibrational decorrelation under two force models” are questions with different initialization, ensemble, length, and validation requirements.

Translate the question into an estimand. A structural mean might be

whereas a population is an indicator average, a diffusion coefficient is a long-time slope, and a rate is a probability flux per unit time. State the ensemble in the mathematical definition only after confirming that Tako’s effective route actually samples or approximates it. If Berendsen coupling is used, do not place an exact canonical or NPT expectation symbol on the resulting average without qualification.

Next identify the slow coordinates. A local bond distribution can converge while a torsional state remains trapped. Lattice volume can equilibrate before defect rearrangement. Temperature may reach its target in a few coupling times while solvation structure requires far longer. The production duration should be justified against these slow coordinates, not merely against the integration timestep.

Choose the force model based on the visited states. If the question permits bond breaking, charge redistribution, unusual coordination, or phase change, a model validated only around equilibrium is insufficient. If NPT is required, stress validation is mandatory. If a charged molecular system uses GFN2, charge, spin, and solvation conventions belong to the scientific definition.

Choose controls in stages. A weak-coupling thermostat can remove initialization transients, followed by a validated sampler for structural averages. Real-time correlations motivate minimal coupling or NVE production. Tako cannot preserve velocities or extended state across calls, so a staged design must presently occur within one fixed-control call by discarding an initial interval, or accept a discontinuous reinitialization. This limitation may make some otherwise standard protocols impossible to execute rigorously through the public surface.

Predeclare acceptance conditions. Examples include a maximum NVE drift, agreement of an observable under timestep halving, a minimum number of independent state transitions, a minimum effective sample size, agreement among replicas, and model-reference force errors below a context-specific tolerance. Acceptance conditions prevent a visually attractive trajectory from being promoted after the fact.

Finally align saved resolution with analysis. High-frequency spectra need dense frames; slow structural means do not benefit from saving every step. Because Tako does not save velocities, sparse positions cannot later recover fast dynamics. Because a nondividing frame interval can omit the final coordinates, select a divisible interval when the final state matters or archive the final structure separately with its step label.

Equilibration and production design

An MD study should distinguish preparation, equilibration, and production. Preparation removes gross structural problems: overlapping atoms, unreasonable cell strain, incorrect protonation, missing charge/spin information, and forces outside the model’s reliable regime. Geometry optimization can help, but overoptimization is not always appropriate when the target is a high-temperature phase or liquid.

Equilibration brings selected observables toward a stationary distribution. Temperature, potential energy, volume, pressure, density, structural order parameters, and composition-specific measures may relax on very different timescales. A flat temperature trace alone does not prove equilibrium. A barostat can make volume appear stationary while slow phase reorganization remains incomplete. A short moving average can hide drift.

Production begins only after equilibration criteria chosen in advance are satisfied. The production thermostat and coupling strengths should reflect the observables. Berendsen coupling can prepare temperature or volume but should not support fluctuation-derived properties. Bussi or NHC may be used for canonical structural averages after validation. Langevin or Andersen can aid sampling but alter kinetics. An NVE segment can preserve unthermostatted time correlations after canonical preparation, provided energy drift is controlled.

Because Tako lacks exact restart, public MD calls cannot cleanly carry final velocities and extended variables from equilibration into production. Starting a second call reinitializes velocities. This is a major methodological limitation. A single call can contain only one fixed set of controls, so it cannot represent a rigorous staged protocol with continuous phase-space state. One may discard an initial portion of a single longer run as equilibration, but the 10,000-step cap and absent full-resolution restart restrict accessible timescales.

Declare discarded time, convergence diagnostics, production duration, frame interval, and whether statistical analysis uses one trajectory or independent replicas. Do not choose the discard point after inspecting the desired result without acknowledging that exploratory choice.

Temperature is not equilibration

Temperature control is one of the easiest MD outputs to overinterpret. The thermostat acts directly or indirectly on kinetic energy, so temperature often reaches its target faster than configurational variables. A flat temperature line demonstrates that the kinetic controller is operating. It does not show that positions sample a stationary distribution, that populations are equilibrated, or that the model represents the physical material.

Consider a torsional molecule initialized in one conformer. Bussi or Langevin dynamics can establish a canonical-looking kinetic-energy distribution while the torsion never crosses its barrier. The reported temperature is satisfactory, yet conformer populations remain entirely determined by the initial state. Similarly, a crystal under NPT coupling can reach a stable mean volume while a reconstructive phase transition is inaccessible on the simulated timescale.

Monitor observables with different physical meanings: potential energy, kinetic temperature, volume or density, key distances and coordination numbers, order parameters, and state labels. Divide the candidate equilibration region into windows and compare means and distributions. A drifting potential-energy baseline at constant temperature indicates configurational relaxation. A stationary average with changing variance can also indicate nonstationarity.

Thermostat damping provides a useful time scale but not an equilibration guarantee. Ten thermostat time constants may be ample for kinetic energy and irrelevant for a microsecond conformational process. Barostat damping likewise controls cell response, not defect diffusion.

The exact Tako initializer makes the first segment particularly noncanonical: it imposes a fixed $3N$-based kinetic energy after COM removal and ignores constraint-adjusted normalization. Discarding some early time is therefore expected, but the discarded duration must be supported by observables. It should not be chosen solely to remove inconvenient values.

When the run is too short to establish equilibration, label it as a relaxation or trajectory demonstration. Short-time mechanistic observations can still be valuable if phrased conditionally: “under this initialization and model, the structure relaxed toward…” is more defensible than “the equilibrium structure is…”.

Convergence dimensions

MD convergence is multidimensional. A calculation can be numerically stable but statistically unconverged, statistically precise on the wrong potential, or apparently equilibrated because a thermostat suppresses fluctuations.

Timestep convergence

Reduce $\Delta t$ while holding initial state and physical protocol as comparable as possible. Examine NVE drift, temperature distribution, structural averages, transition counts, and any final observable. If Andersen is used, preserve collision rate rather than per-step probability. If Langevin is used, friction is already per unit time and can remain fixed.

Trajectory-length convergence

Plot cumulative averages and block estimates versus production length. A stable cumulative mean can be misleading when the trajectory has not crossed relevant barriers. Monitor state populations and transitions, not only the average of a smooth variable.

System-size convergence

Periodic cells introduce image interactions and quantized collective modes. Diffusion, phase transitions, defects, phonons, dielectric response, and interfacial behavior can depend strongly on cell size. Repeat representative calculations with larger cells where feasible. For molecular clusters, box size and boundary choice can change evaporation or confinement.

Model convergence

Compare calculators or model versions only through like-for-like observables; never mix their absolute energies in one thermodynamic expression. Benchmark representative configurations, including high-energy and distorted states visited during MD, against an appropriate independent reference. A model accurate near a relaxed crystal can fail during melting, diffusion, fracture, chemical reaction, or unusual coordination.

Ensemble and coupling convergence

Vary thermostat/barostat damping within a physically reasonable range. Production averages should not be artifacts of coupling. Compare canonical and NVE segments when time correlations matter. Verify the effective Tako route rather than assuming the displayed labels.

Replica convergence

Independent seeds reveal sampling variability and rare-event sensitivity. Agreement across replicas is stronger evidence than a long smooth trace from one seed. Disagreement is information: it may indicate insufficient duration, metastability, or a broad distribution.

Time correlation and effective sample size

MD frames are correlated. If $A_s$ is an observable sampled every $\Delta t_s$, define the normalized autocorrelation

The integrated autocorrelation time is approximately

with a carefully chosen truncation or window. The statistical inefficiency is $g\approx2\tau_{\mathrm{int}}/\Delta t_s$, and the effective number of independent samples is

The raw frame count overstates the effective sample size. Saving frames more frequently reduces neither $\tau_{\mathrm{int}}$ nor the underlying correlation. It can improve resolution of fast observables but increases storage.

Estimate uncertainty using block averaging, autocorrelation-aware standard errors, or replicate-level statistics. Blocks must be longer than the dominant correlation time. For multimodal trajectories, state transitions may be much slower than local fluctuations, and simple autocorrelation estimates within one trapped basin understate uncertainty.

Thermostats influence correlation. Andersen destroys momentum correlations; Langevin damping changes relaxation; Berendsen suppresses fluctuations; NHC can introduce extended-system oscillations. Thus a time-correlation function is a property of both physical model and thermostat dynamics unless collected under an appropriate production protocol.

Structural observables

Positions and cells support many structural analyses. Examples include bond-length and angle distributions, radial distribution functions, coordination numbers, torsions, orientational correlations, mean-squared displacement, lattice parameters, density, and order parameters. Each requires attention to periodic boundaries, atom identity, cell evolution, and equilibration.

For a radial distribution function,

finite bin width, minimum-image conventions, and normalization must be declared. In a changing NPT cell, the density and accessible shell volume vary. For small or anisotropic cells, shells beyond half the shortest box dimension are not independently sampled under a simple minimum-image treatment.

Mean-squared displacement requires unwrapped coordinates:

Tako frames contain wrapped or calculator-propagated positions but no explicit image counters. A diffusion analysis must reconstruct unwrapped motion carefully, especially with changing cells. The Einstein diffusion coefficient in $d$ dimensions,

requires a verified linear diffusive regime, finite-size assessment, uncertainty, and thermostat sensitivity. A 50 fs default trajectory is not adequate.

Conformational populations require repeated transitions or an enhanced-sampling method. Tako currently provides no umbrella sampling, metadynamics, replica exchange, or free-energy estimator in this MD surface. A trajectory trapped in one basin supports local dynamics, not equilibrium population ratios.

Energetic and thermodynamic observables

Tako exposes potential, kinetic, and total particle energy. For NVE, total-energy drift is an integration diagnostic. For thermostatted modes, energy exchange with the bath is expected; the particle total energy is not conserved. Bussi internally accumulates transferred thermostat energy, and NHC/MTK possess extended variables, but Tako does not return these quantities. Consequently, extended-energy conservation cannot be audited through public artifacts.

Heat capacity from energy fluctuations, compressibility from volume fluctuations, and thermal expansion from temperature/volume response demand correct ensemble sampling. Berendsen weak coupling suppresses the very fluctuations these formulas require. Even with an appropriate sampler, finite-size effects, correlation, quantum nuclear corrections, and model accuracy matter.

The instantaneous pressure is absent from the result. NPT algorithms compute stress internally, but Tako does not persist pressure or stress traces. One cannot validate pressure equilibration from md_trace alone. Volume can be derived from frame cells in Å, but remember the trace-cell unit defect. Enthalpy $H=E+PV$ is likewise not returned and cannot be reconstructed exactly without pressure conventions and complete energy terms.

Absolute energy comparisons across calculators are invalid. Even within one model, compare systems with compatible composition or use a justified thermodynamic cycle. MD averages add sampling uncertainty to electronic/model uncertainty. Report both where possible.

Dynamical observables and thermostat perturbation

Velocity autocorrelation, vibrational density of states, diffusion, viscosity, thermal conductivity, reaction rates, and spectral line shapes depend on real-time dynamics. A thermostat can alter them even if configurational averages remain correct.

Tako does not persist velocities. A velocity autocorrelation cannot be computed exactly from public output. Finite differences of positions amplify noise, require consecutive frames, fail across wrapped boundaries without unwrapping, and represent interval-averaged rather than integrator velocities. Sparse frame intervals can erase high-frequency content and alias spectra.

Andersen collisions directly destroy velocity correlations. Langevin friction imposes a damping timescale. Strong Berendsen coupling scales velocities coherently. Bussi global scaling is gentler for some modes but still perturbs kinetics. Nose–Hoover chains add deterministic frequencies that can couple to system modes. For dynamical observables, an established strategy is to equilibrate canonically, then collect an NVE production trajectory with a validated timestep. Tako’s absence of velocity-preserving restart prevents an exact two-call implementation of that strategy.

Reaction events observed in direct MD are conditional on timescale and initial conditions. A single event does not establish a rate, and the absence of events does not prove a barrier is insurmountable. Rare-event kinetics generally require enhanced sampling, transition-state theory, or many long independent trajectories.

Potential-model validity along the trajectory

MD explores geometries that may be absent from a static benchmark. A stable minimum can evolve toward stretched bonds, close contacts, unusual coordination, changed charge localization, or a different phase. The relevant validation set must cover the visited distribution, including tails that control rare events.

For GFN2, monitor SCF convergence, charge/spin appropriateness, elemental coverage, and whether the electronic state remains qualitatively meaningful. A failed or unconverged electronic step invalidates subsequent dynamics. Implicit solvation, if used, is not explicit solvent friction or molecular collision dynamics.

For MLIPs, training-domain coverage is central. Smooth, plausible trajectories can be silently wrong outside the learned manifold. File size does not track accuracy. Compare representative energies, forces, stresses, structures, and possibly uncertainty indicators against matched higher-level calculations. NPT requires particularly reliable stress. Direct-force models that are not exact energy gradients may show inconsistent energy/force behavior and poor NVE conservation; a thermostat can conceal that inconsistency.

Model validation should be iterative. Select frames spanning energy, temperature, coordination, and cell extremes; recompute reference values; quantify errors; and decide whether conclusions survive. Do not validate only the starting structure.

Finite-size and boundary effects

Periodic MD represents an infinite tiling of the simulation cell. A defect interacts with images. Long-wavelength modes longer than the box are absent. Collective transformations can be artificially coherent. Liquids show finite-size corrections to diffusion. Charged periodic systems require electrostatic conventions. Surfaces need sufficient vacuum and appropriate partial periodicity, but Tako’s MD payload promotes any cell to three-dimensional periodicity.

NPT cell fluctuations are especially size dependent: small cells fluctuate strongly and may undergo shape instabilities that a macroscopic material would suppress. Full MTK exposes nine cell degrees in the implementation, including modes that may be inappropriate for a small or low-symmetry cell. Isotropic coupling is more restrictive but may hide anisotropic stress.

Finite-size validation varies with the observable. Compare multiple supercells for diffusion, lattice response, phase stability, and defect energetics. Check that a radial cutoff fits the box. Confirm that correlation lengths are smaller than the cell. For molecular systems intended to be isolated, a nonperiodic structure avoids image interactions, but then NPT is unavailable.

Ensemble-specific validation signatures

Each effective mode has characteristic diagnostics. Applying one generic “temperature near target” test to every mode misses important failures.

For NVE, inspect particle total energy, not potential energy alone. Potential and kinetic terms should exchange substantially while their sum remains bounded around a stable shadow-Hamiltonian value. Fit drift over more than one window, compare halves of the trajectory, and repeat at smaller timestep. A low-amplitude energy oscillation can be benign; monotonic drift is not. Temperature is free to fluctuate and need not equal the initialization target.

For Berendsen NVT, verify a smooth approach toward target without persistent clamp activation or structural distortion. The instantaneous distribution will be too narrow, so do not use agreement of variance with canonical theory as an acceptance target. Compare a later Bussi or NHC segment if canonical fluctuations matter. A coupling time many timesteps long is generally more interpretable than nearly instantaneous rescaling, but no universal value replaces a sensitivity test.

For Bussi NVT, examine the kinetic-energy distribution as well as the mean. In an ideal unconstrained classical canonical system with $f$ active quadratic momentum degrees,

Tako’s COM/DOF mismatch complicates the expected $f$, particularly for small systems. Treat a distribution discrepancy as a prompt to inspect degree counting, equilibration, constraints, and correlation before blaming the thermostat.

For Langevin, verify configurational stationarity across friction values. The target-temperature mean can be correct even when dynamics are overdamped. Diffusion and correlation times should be reported at the actual friction. A larger friction is not simply “more stable”; it defines different stochastic dynamics.

For Andersen, translate the component-wise per-step probability to an approximate rate. Check that results are insensitive to reasonable changes when only configurational equilibrium matters. Do not use collision-interrupted velocity correlations as physical spectra or transport without an explicit stochastic-dynamics interpretation.

For NHC, look for extended-system oscillations, resonances, poor thermalization of decoupled modes, and sensitivity to chain length and damping. Tako does not expose extended energy, so the strongest standard conservation diagnostic is missing. Constrained small systems are especially problematic because the thermostat target uses $3N$.

For Berendsen NPT, inspect volume and structural response but do not validate fluctuation amplitudes against NPT theory. Ensure the cell scale never approaches pathological values and compare final mean pressure using an independent stress evaluation if possible.

For MTK, inspect cell determinant, lengths, angles, density, energy, temperature, and stress behavior. Sudden shear, cell rotation, or collapse can reveal stress/model instability. The absence of a pressure trace and extended energy in Tako means important diagnostics require external recomputation from saved frames.

Pressure and cell diagnostics in the absence of a pressure trace

Tako’s NPT algorithms consume stress internally but do not return a pressure series. This omission complicates verification. A stable volume is not equivalent to target pressure: an overly stiff barostat, suppressed fluctuations, or a model bias can produce a smooth cell with residual stress.

At minimum, compute cell volume for every coordinate frame,

using frame cells, which are in Å. Do not use md_trace.cell without converting from Bohr. Track lengths, angles, determinant sign/magnitude, and density. Full MTK warrants inspection of all cell components because volume alone can hide shear instability. Masked MTK should change only its selected basis-direction lengths according to the implemented meaning of the mask.

Re-evaluate stress for representative saved frames with the same model and a clearly documented sign convention. Then compute scalar pressure from the trace and, where possible, include the kinetic contribution used by the runtime. Because velocities are absent, exact instantaneous kinetic pressure cannot be reconstructed from saved output; only configurational stress or an approximate average based on temperature may be available. State that limitation.

Check unit magnitude before running. One gigapascal is approximately 0.00624151 eV/ų, while one atmosphere is approximately (6.324×10^{-7}) eV/ų. A target entered as “1” means roughly 160 GPa, not 1 bar or 1 GPa. Pressure-sign errors can drive immediate expansion or collapse.

For Berendsen coupling, the supplied compressibility directly affects response. Compare it with the material’s expected compressibility when known. A default numerical value should not be interpreted as a measured property. For MTK, damping masses depend on target temperature and tau_P rather than the Berendsen compressibility, so that field is recorded but unused.

Finite cells produce large pressure and volume fluctuations. Assess whether these are physical, finite-size enhanced, or numerical. A barostat should not be used to force a mechanically unstable phase to remain near an expected lattice parameter; instability is a result only after model and numerical validation.

Stability failures and diagnostic patterns

An MD failure is often a symptom rather than a single error. Useful categories are:

A failed SCF or missing stress is a hard scientific failure, not a missing optional field. A completed status proves execution, not equilibrium or validity. Preserve logs, requested settings, effective route, model identity, and the first failing geometry.

Statistical uncertainty beyond a standard error

An MD uncertainty statement should separate at least four components: finite correlated sampling, replica-to-replica variation, model error, and protocol sensitivity. A tiny standard error from one long correlated series can coexist with a large model bias. Conversely, several short replicas may reveal metastability that a single nominally precise trajectory hides.

For a stationary unimodal observable, estimate an autocorrelation-aware standard error,

Use block lengths exceeding the integrated correlation time and demonstrate a plateau of the estimated uncertainty with increasing block size. If no plateau appears, the run is too short for that estimator. Report the number of blocks and do not treat block means from overlapping windows as independent.

For multiple replicas, use replica means as the independent units when each run has a different seed. The between-replica standard deviation captures slow-state sensitivity better than pooling every frame. When replicas have unequal effective lengths, weighting requires care; a long trapped trajectory should not automatically dominate shorter trajectories sampling different basins.

For multimodal systems, report state populations and transition counts. A bootstrap over individual frames is invalid because it destroys time correlation. A block bootstrap can help when blocks encompass local correlation but still fails if no interbasin transitions occur. In that case the honest conclusion is that equilibrium populations are unresolved.

Model error requires independent evidence. Select decorrelated frames across the sampled range and compare reference energies, forces, or stresses. Summaries such as force RMSE should be stratified by state, element, coordination, and energy because a global average can conceal failure in the rare configurations that control kinetics. When two models disagree beyond sampling uncertainty, do not average them as if they were replicate measurements; investigate the physical origin.

Protocol sensitivity includes timestep, thermostat strength, barostat choice, cell size, constraints, initialization, and discarded equilibration. Present a small sensitivity table with the target observable, not only runtime diagnostics. A result is robust when the scientific conclusion persists across defensible protocol variants.

For qualitative stability claims, uncertainty still matters. “No decomposition was observed” should be accompanied by simulated time, number of independent replicas, temperature, model, and accessible mechanism. It means only that no event occurred within that finite exposure, not that the lifetime is long experimentally.

Worked reasoning patterns without a software recipe

The following patterns illustrate how to reason from a scientific objective to an interpretable protocol. They are not universal parameter prescriptions and deliberately avoid interface or scripting instructions.

Local molecular flexibility

Suppose the question concerns the distribution of one torsion in a neutral organic molecule near room temperature. GFN2 with an explicitly stated dispersion treatment may provide an exploratory surface if the elements and charge state are appropriate. A fixed-cell molecular structure avoids irrelevant pressure coupling. Bussi or weak Langevin dynamics may provide configurational sampling, but the initial $3N$ rescaling and COM convention are relatively important for a small molecule. The analysis should discard kinetic equilibration, verify repeated torsional transitions, compare independent seeds, and report a population only if both directions are visited sufficiently. If the torsion never crosses, the result describes fluctuations within the starting conformer, not an equilibrium conformer ratio. Classical treatment of high-frequency vibrations and possible proton tunneling belongs in the limitations.

Crystal short-time stability

Suppose the question is whether a proposed periodic crystal remains mechanically coherent under one MLIP at a target temperature. First validate energy, force, and stress on representative distorted cells. A geometry relaxation can reduce excessive initial stress. A short NVT equilibration followed by conservative dynamics would normally separate thermal preparation from energy-drift assessment, but Tako cannot carry velocities exactly between calls. One continuous Bussi run can test finite-time structural integrity, while a separate NVE call with the same structure and declared reinitialization can test timestep stability. Neither proves thermodynamic phase stability. Monitor coordination, cell-fixed density, displacement distributions, and model-domain indicators. Replicate seeds distinguish deterministic fragility from one unusual initial velocity draw.

Pressure-driven cell response

Suppose the question concerns qualitative compression of a bulk cell. The target pressure must be converted to eV/ų, and calculator stress must be benchmarked. Berendsen NPT can screen smooth isotropic response but cannot support fluctuation-derived compressibility. Isotropic MTK samples only volume-like changes; full MTK admits shape and shear; masked MTK changes selected basis-direction lengths. Atomic constraints are incompatible with all current MTK routes. Because control combinations can reroute, the requested thermostat must be chosen so the intended barostat branch actually wins. Analyze cell determinant, lengths, angles, density, and independently recomputed stress. Treat a cell collapse as either physical instability or a numerical/model failure until timestep, damping, stress convention, and finite size are excluded.

Diffusion screening

Suppose the question asks whether an ion is mobile in a host. A trajectory needs unwrapped positions, repeated displacement events, and physical time far beyond vibrational motion. The 10,000-step cap gives only 5 ps at 0.5 fs, often inadequate. Public outputs lack image counters and velocities, and staged calls are discontinuous. A nonzero mean-squared displacement over a few frames may be vibration or a single hop, not diffusion. Any estimated slope needs multiple time origins, a verified diffusive regime, correlation-aware uncertainty, cell-size correction, and replicas. Langevin and Andersen alter transport; strong coupling makes a quantitative diffusion coefficient thermostat dependent. Under current constraints, Tako is better suited to exploratory mobility signals than a publishable transport coefficient unless the accessible timescale is demonstrably sufficient.

Model comparison

Suppose two MLIPs are compared. Use the same initial structures, declared seeds, timestep, and effective mode, but recognize that identical random numbers do not create identical phase-space paths once forces diverge. Compare distributions and observables, not frame-by-frame coordinates after the initial divergence. Benchmark both models against the same independent reference frames. Do not combine their absolute energies or treat the difference between two model means as uncertainty without reference validation. A larger model or smoother trajectory is not necessarily more accurate.

These patterns share a discipline: define what the finite trajectory can test, identify what the current runtime omits, and scale the conclusion to the evidence. The most valuable result can be a negative one—for example, discovering that the accessible time is too short, stress is unreliable, or the effective route differs from the requested ensemble.

Reconstructing provenance from current artifacts

The calculation folder separates settings, streamed frames, scalar trace, result summary, and logs. Reproducibility depends on joining them correctly. The result settings describe normalized requested values. They do not reveal the effective dispatcher branch. The log announces requested ensemble, thermostat, and barostat, also not the effective mode. Therefore the effective route must presently be reconstructed using the routing table and the software revision.

Frames contain the authoritative saved coordinates and Å cell. The md_trace contains the authoritative final scalar point even when no final frame exists, but its cell is in Bohr. energy_trace.json is generated from coordinate frames rather than from the complete scalar trace, so it can omit the final step. A robust analysis should identify the last step in each artifact before combining them.

The final atoms array contains positions but no velocity and no top-level cell. For fixed-cell runs, the input cell remains applicable. For cell dynamics, use the last coordinate frame’s cell if it corresponds to the final step. If the final frame is absent, the trace has an internal-unit cell but no final coordinates; conversion recovers the cell, but exact final phase-space continuation remains impossible.

Model provenance should include more than level_of_theory. MLIP runtime and model name are returned, but a durable record should also identify weight hashes or versioned assets. GFN2 records should retain charge, unpaired/spin settings, dispersion, solvation, and SCF policy from the settings artifact. A result generated after an asset update may not be bitwise reproducible from a generic model label.

Logs are necessary for detecting SCF difficulty, warnings, and failures, but successful logging does not establish physical validity. Screenshots and plots should link back to machine-readable artifacts and exact steps. When citing a frame, include calculation ID or workspace path, frame step, physical time, and whether its cell came from frame or converted trace data.

Reproducibility record

The shared record skeleton and reporting standard are defined in Reproducibility and Reporting. A reproducible Tako MD report should additionally include at least:

Archive the canonical workspace and machine-readable traces. Record whether the last coordinate frame equals the final scalar step. A screenshot is orientation evidence, not a substitute for settings or numerical artifacts.

A minimum scientific acceptance checklist

Before treating a trajectory as evidence, require affirmative answers to the following questions:

1. Is the calculator appropriate for every important configuration visited, not just the start?
2. Are charge, spin, dispersion, solvation, model weights, and stress capabilities declared?
3. Is the effective runtime route known from the requested combination?
4. Is the timestep converged for stability and the target observable?
5. Are initial overlaps, residual forces, and cell strain controlled?
6. Are constraints compatible with the effective mode and reflected in interpretation?
7. Is the COM/DOF temperature convention acceptable for the system size?
8. Has equilibration been demonstrated with more than temperature alone?
9. Is production long enough relative to autocorrelation and transition times?
10. Are uncertainty and effective sample size reported rather than raw frame count?
11. Are finite-size and boundary effects assessed?
12. Are thermostat/barostat perturbations compatible with the observable?
13. Are frame and trace units handled correctly?
14. Is a missing final coordinate frame recognized when the interval does not divide the step count?
15. Is the absence of public restart respected rather than hiding discontinuous concatenation?

Failure of one item does not always make exploratory dynamics useless, but it changes the permissible claim. A 50 fs default run may demonstrate that a model executes and a structure remains intact. It does not establish equilibrium, diffusion, stability over experimental times, or a thermodynamic property.

Reporting language proportional to evidence

MD prose should distinguish observation, inference, and extrapolation. “The O–H distance remained between specified bounds during 5 ps of Bussi dynamics under model X” is an observation. “The local bonded structure is metastable on the sampled model surface over that interval” is a bounded inference. “The material is stable at room temperature” is an extrapolation requiring far more evidence.

For equilibrium language, prefer “the production segment produced a stationary estimate under the effective Bussi-NVT route” over “the system equilibrated” unless stationarity, transitions, replicas, and uncertainty support the latter. For kinetics, state thermostat perturbation and avoid translating one observed event directly into a rate. For absence of events, report exposure time and replica count.

When requested and effective settings differ, do not hide the discrepancy. A complete statement might say that the archived request named NPT/Nose–Hoover/no barostat, while the audited dispatcher selected isotropic MTK NPT. This is scientifically more useful than reproducing a misleading settings table.

When artifacts disagree in endpoint coverage, specify which one supports each number. Final energy and temperature can come from the final scalar trace; the last available geometry may be an earlier frame. If cell values are taken from md_trace, state that they were converted from Bohr. Never silently align arrays by index when their step labels differ.

Report failures and exclusions. If trajectories were discarded for SCF failure, cell collapse, leaving the model domain, or nonstationarity, give counts and prespecified criteria. Selective removal of unfavorable replicas biases estimates. A completed browser status is an execution record, not a scientific acceptance flag.

Good reporting makes uncertainty visible without making the result unreadable. Put the question, effective method, physical duration, evidence level, headline observable with uncertainty, and main limitation near the beginning. Place routing detail, convergence tables, and provenance in method/reference sections. This organization lets a human reader understand the claim before confronting every control while preserving the audit trail needed for expert review.