Molecular Dynamics

Molecular dynamics (MD) replaces a single optimized geometry with a time-ordered trajectory. Given atomic positions, masses, velocities, and a model for the potential energy, an integrator advances the system through phase space. The resulting sequence can reveal motion, relaxation, diffusion, conformational exchange, temperature response, and fluctuations that no stationary structure can show. It can also produce misleading results when the potential is inappropriate, the timestep is unstable, the run is too short, the ensemble is not what the settings imply, or correlated samples are treated as independent observations.

Tako currently couples its calculation models to the dynamics implementation in ase.rs. The browser-facing settings are more permissive than the effective runtime routing. Consequently, a requested ensemble, thermostat, and barostat can name one combination while the dispatcher executes another. That distinction is central to interpretation. The exact routing table appears later in this chapter, and it should be consulted before assigning an ensemble label to a trajectory.

Statistical design, equilibration, convergence, observables, and reporting standards are developed separately in Molecular Dynamics: Statistics and Validation.

The scientific object: a trajectory on a potential-energy surface

For $N$ nuclei with Cartesian positions $\mathbf R=(\mathbf R_1,\ldots,\mathbf R_N)$ and masses $m_i$, classical MD integrates Newton’s equations

The function $U(\mathbf R)$ is the potential-energy surface supplied by the selected calculator. In a Born–Oppenheimer interpretation, the electronic problem is solved or approximated at each nuclear geometry, the nuclei move on one electronic surface, and nonadiabatic transitions are absent. GFN2-xTB supplies an approximate self-consistent electronic energy and forces. A machine-learned interatomic potential supplies an energy/force model inferred from training data. These surfaces are not interchangeable and do not share an absolute energy reference. A trajectory is therefore always conditional on its calculator, model version, numerical settings, charge and spin treatment where applicable, and the region of configuration space actually visited.

Classical nuclei are another approximation. Zero-point motion, tunneling, nuclear quantum statistics, and electronic nonadiabaticity are not recovered merely by choosing a low timestep or a sophisticated thermostat. Classical equipartition can overpopulate high-frequency molecular vibrations at temperatures where quantum occupation would be small. Proton transfer, isotope effects, cryogenic motion, and electronically excited dynamics require particular caution. Tako does not currently provide path-integral dynamics, surface hopping, multiple-time-step quantum dynamics, or a quantum thermostat.

The primary numerical object is not a single energy but a discretized path

where $\mathbf P$ denotes momenta and $\mathbf h$ the periodic cell when present. Tako exposes positions and cells for stored frames, plus scalar kinetic, potential, and total energies and an instantaneous kinetic temperature. It does not expose velocities, momenta, extended-Hamiltonian terms, or a public restart object. Thus the persisted trajectory is sufficient for many structural analyses but not for exact continuation or every time-correlation observable.

Phase space, conserved quantities, and ensembles

The phase-space point for fixed-cell classical dynamics is $(\mathbf R,\mathbf P)$. The Hamiltonian is

An ideal microcanonical, or NVE, trajectory conserves particle number $N$, volume $V$, and total energy $E$. It is generated by Hamilton’s equations without a thermostat or barostat. Numerically, a finite-timestep symplectic method does not conserve the exact continuous Hamiltonian at every step; it approximately conserves a nearby shadow Hamiltonian. Bounded oscillation of total energy is expected. Persistent drift, abrupt jumps, or systematic heating indicate integration, force, constraint, or model problems.

The canonical NVT ensemble fixes $N$, $V$, and target temperature $T$. Its equilibrium density is proportional to $\exp(-\beta H)$ with $\beta=(k_B T)^{-1}$. A thermostat is a dynamical sampling device, not merely a temperature correction. Different thermostats produce different transient behavior, dynamical perturbations, ergodicity properties, and fluctuation distributions. Berendsen weak coupling suppresses fluctuations and does not generate the exact canonical ensemble. Bussi stochastic velocity rescaling is designed to recover canonical kinetic-energy statistics. Langevin and Andersen dynamics are stochastic and alter real-time correlations. Nose–Hoover chains add deterministic extended variables and can sample the canonical distribution under suitable ergodicity conditions.

The isothermal–isobaric NPT ensemble permits volume or cell changes at target pressure $P$ and temperature $T$. Its implementation requires stress from the calculator and additional cell degrees of freedom. Berendsen pressure coupling is useful for controlled relaxation but does not produce correct NPT fluctuations. Martyna–Tobias–Klein (MTK) dynamics introduces extended thermostat and barostat variables intended to sample an isothermal–isobaric distribution. Isotropic, full-cell, and axis-masked variants represent materially different mechanical assumptions.

An ensemble name is scientifically meaningful only when the effective equations, conserved or stationary distribution, cell freedoms, and sampling stage are stated. A short thermostatted trajectory may be an equilibration procedure rather than a sample of equilibrium. An NVE production segment following NVT equilibration may preserve time correlations better than a strongly stochastic production run. Conversely, an NVE segment begun from a poorly equilibrated configuration merely conserves the wrong energy shell.

From continuous equations to a discrete integrator

Write velocity as $\mathbf v_i=d\mathbf R_i/dt$ and acceleration as $\mathbf a_i=\mathbf F_i/m_i$. Taylor expansion gives

and a corresponding expansion around $t+\Delta t$. Combining forward and backward information yields the velocity-Verlet update:

The first half-kick uses forces at the old geometry. The drift constructs a trial geometry. Forces are recomputed there, and the second half-kick completes the velocity. Tako’s ase.rs implementation follows this structure and caches forces between ordinary Verlet-like steps. When constraints alter the trial positions, velocities are reconstructed from the realized displacement before the second half-kick.

Velocity Verlet is time reversible in exact arithmetic, symplectic for an unconstrained conservative Hamiltonian, and second-order accurate globally. Those properties explain its widespread use. They do not make an arbitrary timestep safe. The shortest important period in the system sets the demanding scale. Light-atom bond stretches, especially X–H modes, require smaller steps than heavy-atom lattice motion. Stiff repulsive encounters created by a poor initial geometry can demand a much smaller step than the equilibrated structure. No thermostat repairs a step that fails to resolve the forces.

The finite-step trajectory follows a shadow Hamiltonian

under ideal conditions. Energy may oscillate around a stable mean even when the exact $H$ fluctuates. Reducing $\Delta t$ should reduce the amplitude of this integration error approximately quadratically in a regular regime. A timestep-convergence study is therefore a physical validation: compare energy drift, temperature distribution, structural observables, and event rates at successively smaller steps. Agreement of a single final energy is insufficient.

Timestep selection and energy drift

Tako accepts a positive timestep in femtoseconds and defaults to 0.5 fs. The runtime converts femtoseconds to atomic time before propagation. There is no automatic frequency analysis, hydrogen-mass repartitioning, SHAKE/RATTLE bond constraint, multiple-time-step integrator, or adaptive timestep. The user must connect the timestep to the fastest active motion and the smoothness of the chosen force model.

A defensible test starts from the same coordinates and velocity seed, runs short NVE segments at $\Delta t$, $\Delta t/2$, and if necessary $\Delta t/4$, and compares a normalized drift such as

or a relative measure scaled by the magnitude of kinetic energy. A linear fit should exclude the initial transient only when that exclusion is declared. Report both the oscillation envelope and fitted drift because an unstable integrator may show bursts hidden by a single endpoint difference. For thermostatted dynamics, total particle energy is not conserved; inspect the expected temperature distribution, coupling behavior, and, when available, the appropriate extended energy instead. Tako does not expose extended energy, so an NHC or MTK conservation test cannot presently be completed from public results alone.

Warning signs include monotonically rising temperature in NVE, alternating explosive bond lengths, atoms traversing large fractions of the cell in one step, discontinuous force changes, SCF failures concentrated at later steps, or sensitivity of qualitative outcomes to a modest timestep reduction. A stochastic thermostat can mask deterministic energy drift by removing injected heat. Therefore thermostat stability does not establish integrator stability. Validate the underlying model and timestep with a conservative segment whenever the scientific question permits.

The default run length is only 100 steps, or 50 fs at 0.5 fs. That may demonstrate execution or resolve several high-frequency vibrations, but it is usually too short for thermodynamic averages, diffusion, conformer populations, solvent reorganization, phase stability, or reliable autocorrelation estimates. Increasing the number of printed frames without increasing physical time does not improve statistical sampling.

Initial velocities in statistical mechanics

For a classical canonical distribution, each unconstrained Cartesian momentum component has a Gaussian distribution

Independent Gaussian draws produce a Maxwell–Boltzmann speed distribution. A finite sample does not have exactly the target kinetic energy. Many programs optionally rescale the draw, remove center-of-mass translation, remove overall rotation for an isolated molecule, and project constraints. Each operation changes the distribution and the appropriate degree count. Exact rescaling replaces canonical kinetic-energy fluctuations with a fixed initial value; this can be acceptable initialization, but it should not be confused with an equilibrium sample.

For constraints removing $n_c$ independent Cartesian components, a basic translational degree count is

If center-of-mass translation is constrained, three more degrees are normally removed for a nondegenerate multi-atom system. Removing overall rotation can remove up to three additional degrees for a nonlinear isolated molecule or two for a linear molecule. Periodic solids are generally not treated with molecular rotational subtraction. The correct count depends on the actual constraints and statistical ensemble.

The instantaneous kinetic temperature is conventionally

It is an estimator, not a statement that each frame is thermally equilibrated. For small systems it fluctuates strongly. Temperature is an ensemble property inferred over time or replicas; a single instantaneous value can be far from the target without indicating failure.

Exact Tako velocity initialization

Tako does not preserve incoming velocities for a public MD call. Before selecting either the GFN2 or MLIP calculator, the WASM bridge creates a new velocity vector from a deterministic Gaussian generator seeded by mdSeed. Each component is first drawn with a width proportional to (m_i^{-1/2}). The bridge then removes the mass-weighted center-of-mass velocity unconditionally and rescales the remaining velocities to

This exact behavior is implemented in crates/tako_calculation_wasm/src/lib.rs. It has four consequences that must be made explicit in scientific reporting.

First, mdFixCom=false does not retain a randomly drawn center-of-mass velocity. Initial COM motion is removed for every mode. The setting controls only some later algorithms.

Second, the rescaling uses $3N$ even though three COM components were just removed. The remaining internal velocity subspace therefore receives the kinetic energy associated with all $3N$ components. If temperature is later computed with denominator $3N$, the printed initial temperature can equal the target, but the active $3N-3$ components carry excess energy relative to a COM-constrained canonical distribution.

Third, the initializer does not use the constraint-adjusted degree count. Selective atomic constraints have already been attached as FixAtoms or FixCartesian, but the target remains $3N$, and the freshly generated velocities are not explicitly projected through the general fixed-velocity routine before the rescaling. Constraints are applied during propagation, so the first integration step can remove components after an overinclusive kinetic normalization.

Fourth, exact rescaling destroys the canonical fluctuation in initial kinetic energy. Stochastic or extended-system thermostats can subsequently establish their stationary distribution, but an equilibration interval is necessary. A short trajectory beginning immediately after this initialization should not be analyzed as if its first frame were canonical.

The reported temperature uses ase.rs degrees_of_freedom(), which subtracts explicit constraint removals but not the removed COM components. This mismatch between initialization, active subspace, and temperature denominator is especially significant for small molecules and heavily constrained structures.

The user-facing registry entry for this limitation is in Known Limitations.

NVE and velocity-Verlet dynamics

The cleanest interpretation of Tako’s NVE route is fixed-cell velocity Verlet on the selected potential-energy surface. It has no thermostat and no cell evolution. It is appropriate for testing force continuity and timestep stability, observing approximately energy-conserving motion, or collecting dynamics where thermostat perturbations would be undesirable. It is not necessarily an equilibrium ensemble: the starting configuration and velocities must already represent the intended energy shell.

In Tako routing, any request with ensemble=nve becomes velocity Verlet, regardless of the thermostat or barostat fields. In addition, any request with thermostat=none becomes velocity Verlet before NPT routing is considered. This means a nominal npt + none + mtkFull request performs fixed-cell Verlet, although its returned settings retain NPT and MTK labels. The scientific method section of a report must therefore state the effective route, not merely repeat the submitted names.

The mdTemperatureK value still matters for NVE because Tako uses it to generate and rescale initial velocities. It is an initialization temperature, not a maintained bath temperature. Subsequent kinetic temperature can change as potential and kinetic energy exchange. The mdFixCom setting does not trigger explicit COM removal within the ordinary Verlet step. Initial COM has already been removed. For a closed unconstrained system with translationally invariant forces, total momentum should remain nearly constant, but constraints, numerical effects, or a non-invariant model can complicate that expectation.

For validation, inspect total energy versus time, momentum or COM drift if velocities are available externally, structural integrity, and timestep convergence. Tako does not persist velocities, so an exact momentum trace cannot be reconstructed from public artifacts. Finite-difference velocities inferred from stored positions are approximate and depend on frame interval, periodic unwrapping, and cell motion; they should not be silently substituted for integrator velocities.

Berendsen temperature coupling

Berendsen weak coupling rescales velocities toward a target temperature. Tako computes the instantaneous kinetic temperature from current kinetic energy and constrained Cartesian DOF, then uses

with the scale clamped to the range 0.9–1.1. The rescaling occurs before each velocity-Verlet step. If the current temperature is nonpositive, the scaling routine returns without intervention. When mdFixCom is true, COM velocity is removed after rescaling, and fixed velocity components are zeroed.

The coupling time tau_T controls how rapidly temperature approaches the target. A very short time produces aggressive, nonphysical damping and can conceal poor initialization or unstable integration. A long time perturbs motion less but equilibrates slowly. The clamp adds another nonlinearity when the requested correction would be extreme.

Berendsen coupling does not generate the correct canonical kinetic-energy distribution. It suppresses natural fluctuations and can produce the flying-ice-cube pathology, where energy is redistributed among modes in unphysical ways. It is best treated as an equilibration or stabilization method. Production averages that depend on canonical fluctuations, heat capacity, or correct kinetic-energy statistics should use a sampler designed for the canonical ensemble after the system has reached the target regime.

The original weak-coupling method is described by Berendsen et al., J. Chem. Phys. 81, 3684–3690 (1984), DOI: 10.1063/1.448118. Tako’s clamp and its degree-count/COM choices are implementation details beyond the idealized method and should be recorded when reproducing a run.

Bussi stochastic velocity rescaling

Bussi, Donadio, and Parrinello introduced stochastic velocity rescaling to retain the convenience of global scaling while sampling the canonical kinetic-energy distribution. Instead of deterministically forcing $K$ toward its mean, the scaling factor contains Gaussian and chi-square random variables chosen so the stationary kinetic-energy distribution is correct.

Tako’s ase.rs state computes

where $f=\texttt{atoms.degrees\_of\_freedom()}$. At each step it draws one normal variate and a chi-square variate with $f-1$ degrees, constructs a global scaling factor, rescales all velocities, optionally removes COM motion, zeros fixed components, and accumulates the thermostat’s transferred energy internally. The step then continues with velocity Verlet.

Bussi requires nonzero initial kinetic energy. Tako’s initializer normally guarantees that for a nonempty, not-fully-constrained system. A system with no active degree or velocities annihilated by constraints can still fail.

Tako defaults mdFixCom to true and overwrites the ase.rs constructor’s mode-specific default. However, the Bussi target $K_0$ uses the constraint-adjusted $f$ without subtracting COM degrees. If COM is removed after scaling, the thermostat and active subspace are not counted consistently. The mismatch becomes appreciable for small $N$.

Stochastic velocity rescaling is often a good production thermostat for equilibrium structural averages because it perturbs all velocities coherently and preserves the canonical kinetic distribution. It still changes real-time dynamics. Transport coefficients and frequency-dependent correlations require sensitivity tests with weaker coupling or an NVE production segment. The foundational reference is Bussi, Donadio, and Parrinello, J. Chem. Phys. 126, 014101 (2007), DOI: 10.1063/1.2408420.

Langevin dynamics

Langevin dynamics augments Newton’s equation with friction and random force:

with noise covariance satisfying fluctuation–dissipation balance,

Tako accepts friction in fs⁻¹, converts it to atomic units, and uses a finite-step scheme with correlated stochastic increments for positions and velocities. The algorithm draws six Gaussian variates per atom per step, constructs the random displacement and velocity terms, performs force/friction half updates, drifts positions, recomputes forces, and completes the update.

When mdFixCom is true and $N>1$, the implementation removes the random center-of-mass components with a finite-size correction and removes COM velocity again after the final update. This is one of the modes where the setting has an ongoing effect. Constraints adjust forces and realized positions, but the stochastic generator is initially constructed for every Cartesian component; heavily constrained dynamics deserves explicit verification.

Large friction rapidly thermalizes but overdamps motion, changes diffusion, broadens or suppresses spectral features, and destroys inertial dynamics. Small friction perturbs motion less but may equilibrate slowly. A thermostat value should therefore be justified by the observable. Equilibrium configurational averages can be robust over a range of friction after adequate sampling, while dynamical observables are directly friction dependent.

The seed deterministically controls the Langevin noise stream in addition to controlling the separately generated initial velocities. Independent replicas require distinct declared seeds. Repeating a single seed demonstrates numerical reproducibility, not statistical precision. Langevin’s stochastic formulation traces to Langevin’s 1908 treatment and modern finite-step algorithms; a commonly relevant integration analysis is Vanden-Eijnden and Ciccotti, Chem. Phys. Lett. 429, 310–316 (2006), DOI: 10.1016/j.cplett.2006.07.086.

Andersen thermostat

The Andersen thermostat represents collisions with an ideal heat bath by occasionally replacing velocities with fresh Maxwell–Boltzmann draws. In the usual conceptual description, particles experience stochastic collisions at a rate related to a collision frequency. Between collisions they follow conservative dynamics.

Tako implements a discrete probability per Cartesian velocity component per MD step. For every atom and each of (x,y,z), it draws an independent collision decision. When accepted, that one component is replaced by a Gaussian velocity of width $\sqrt{k_B T/m_i}$. The implementation therefore does not make one collision decision per atom and replace its complete three-dimensional velocity.

The user-facing mdAndersenProbability is clamped to [0,1]. It is not a rate in fs⁻¹. Changing the timestep while holding the per-step probability fixed changes the physical collision frequency approximately as

for a component. Thus timestep convergence should either keep the intended rate fixed by changing $p$, or explicitly acknowledge that thermostat strength changed.

With mdFixCom=true, Tako removes COM velocity after the first collision stage and again after the final force kick. Andersen collisions efficiently destroy correlations and can improve canonical configurational sampling in systems with ergodicity problems. They also strongly disrupt momentum conservation, hydrodynamic behavior, diffusion, and vibrational time correlations. Andersen dynamics is generally unsuitable when realistic kinetics or transport is the target.

The original method is H. C. Andersen, J. Chem. Phys. 72, 2384–2393 (1980), DOI: 10.1063/1.439486. A Tako report should additionally state that collisions are component-wise and give both $p$ and $\Delta t$.

Nose–Hoover-chain dynamics

The Nose–Hoover approach embeds the physical system in an extended deterministic dynamics whose additional variables exchange energy with particle momenta. A single thermostat variable can be nonergodic for small or nearly harmonic systems. Chains attach further thermostat variables to the first, improving energy exchange across the extended system.

Tako’s NHC propagator uses chain coordinates $\eta_j$, momenta $p_{\eta_j}$, and masses $Q_j$. The first mass is initialized proportional to

while later chain masses are $k_B T\tau_T^2$. The first driving term likewise compares twice the particle kinetic energy with $(3N)k_B T$. The propagator uses a symmetric three-coefficient fourth-order factorization, repeated mdTloop times.

The use of $3N$, rather than atoms.degrees_of_freedom(), is a significant Tako limitation. Explicit Cartesian constraints are permitted in NHC mode, but they do not reduce the thermostat’s target count. mdFixCom is accepted and stored, yet the NHC step does not remove COM velocity. Only Tako’s unconditional initial removal applies. The extended dynamics can therefore regenerate or retain COM drift without the checkbox enforcing ongoing correction.

NHC dynamics can preserve deterministic trajectories and is often chosen for canonical sampling or dynamics less directly randomized than Langevin/Andersen. It still perturbs physical time correlations, especially with short damping time, and ergodicity is not guaranteed. Chain length and subloop count are numerical and sampling parameters. Longer chains are not necessarily more accurate; they add extended modes and cost. Tako clamps chain lengths and loop counts to 1–32.

The relevant foundations include Nosé, J. Chem. Phys. 81, 511–519 (1984), DOI: 10.1063/1.447334; Hoover, Phys. Rev. A 31, 1695–1697 (1985), DOI: 10.1103/PhysRevA.31.1695; and Martyna, Klein, and Tuckerman, J. Chem. Phys. 97, 2635–2643 (1992), DOI: 10.1063/1.463940.

Pressure, stress, and periodic cells

Pressure coupling requires a periodic cell $\mathbf h$ and a calculator capable of returning a $3\times3$ stress tensor. Tako rejects an NPT request when the structure has no cell or has no periodic boundary. The frontend, however, marks all three axes periodic whenever a cell exists, so partial slab or chain periodicity is promoted to bulk periodicity in the MD payload.

The scalar configurational pressure convention used by Berendsen coupling is

with an ideal-gas kinetic contribution added. The sign convention matters: positive target pressure and the calculator’s stress convention must be mutually consistent. MLIP inference requests JAX/ASE-compatible stress and converts it into Hartree per Bohr cubed internally. Missing stress is a hard failure. A model advertised for energy and forces need not be suitable for cell dynamics; its stress must be implemented, smooth, and validated.

Tako accepts target pressure in eV/ų. Useful conversions include

so the default zero is the internal representation of approximately zero external pressure, not one atmosphere. A value that looks numerically small in eV/ų may correspond to a large laboratory pressure. Compressibility is entered in ų/eV, damping times in fs, and friction in fs⁻¹. The runtime converts these to atomic units.

Cell validity requires positive nonzero volume. MTK modes reject degenerate cells and all explicit atomic constraints. Full-cell motion can change lengths, angles, orientation-like components, and shape; it needs a mechanically stable model and enough system size to avoid pathological fluctuations. Isotropic motion preserves shape. Masked MTK evolves selected lattice-vector length directions in a basis derived from the initial cell; it is not a general stress-component mask and does not independently select shear strains.

Berendsen pressure coupling

Tako’s NPT Berendsen route first applies Berendsen velocity scaling, then rescales the cell isotropically, then performs a velocity-Verlet step. The cell scale is

where $\kappa$ is the supplied compressibility and $\tau_P$ the pressure-coupling time. The implementation rejects nonfinite or nonpositive scale factors. Atomic coordinates are scaled with the cell.

This is isotropic weak coupling even when the cell is anisotropic. It cannot independently equilibrate lattice lengths or shear. It drives pressure and volume smoothly but suppresses their equilibrium fluctuations, so compressibility or fluctuation-derived thermodynamic properties cannot be estimated correctly from the resulting trajectory. As with temperature weak coupling, it is most defensible as an equilibration stage.

The compressibility parameter controls response strength in addition to tau_P. A generic default is not a material property. Too large a compressibility or too short a damping time can produce violent volume oscillation or invalid scaling. Too weak a coupling can leave residual pressure for the entire short run. Cell relaxation before MD and a staged pressure equilibration are often safer than asking the barostat to repair a badly strained initial cell.

Tako requires positive compressibility for this route but not for MTK. The returned settings report the supplied value whether or not the effective route uses it. Berendsen et al. (1984), DOI: 10.1063/1.448118, is the primary weak-coupling reference.

MTK isothermal–isobaric dynamics

MTK dynamics augments particle phase space with cell and thermostat/barostat-chain variables and corrects the measure of earlier extended-system NPT formulations. Tako implements isotropic, full, and masked forms. Each contains a Nose–Hoover-chain thermostat for particles and a chain controlling barostat momenta. The particle damping time, barostat damping time, chain lengths, and subloop counts determine extended masses and integration resolution.

The isotropic state represents volume through a logarithmic scalar epsilon:

This preserves cell shape. Full MTK stores a 33 cell momentum and propagates the whole cell matrix. Masked MTK stores scalar momenta along selected normalized initial cell columns. The number of barostat degrees is nine for full-cell mode and the number of enabled axes for masked mode.

MTK propagation requests forces and stress repeatedly, includes kinetic stress contributions, integrates thermostat and barostat chains with symmetric splitting, evolves particle coordinates under cell motion, and updates the cell. It is substantially more demanding than fixed-cell MD. Stress noise, an unsuitable MLIP, loose electronic convergence, or a poor initial cell can destabilize it even when fixed-cell forces look reasonable.

All MTK variants reject atomic constraints in the current ase.rs implementation. mdFixCom has no ongoing effect inside the MTK step. The requested thermostat name is ignored whenever an explicit MTK barostat route wins; MTK always uses its internal NHC machinery. For example, npt + langevin + mtkFull does not combine Langevin particles with an MTK cell. It runs full MTK with NHC particle thermostatting.

The primary formulation is Martyna, Tobias, and Klein, J. Chem. Phys. 101, 4177–4189 (1994), DOI: 10.1063/1.467468. Symmetric extended-system integration is further developed by Tuckerman, Alejandre, López-Rendón, Jochim, and Martyna, J. Phys. A 39, 5629–5651 (2006), DOI: 10.1088/0305-4470/39/19/S18.

Exact requested-to-effective routing in Tako

The runtime lowercases the submitted strings and applies a first-match dispatcher. The order is scientifically consequential. The table below is exhaustive by equivalence class; it covers every possible combination of the three public ensembles, six thermostats, and five barostats.

Three combinations deserve special emphasis:

• npt + none + mtkFull is fixed-cell Verlet because thermostat=none is matched before any barostat.
• npt + noseHoover + none is isotropic MTK, introducing cell motion despite a requested none barostat.
• npt + langevin + none is Berendsen NPT, replacing Langevin and introducing a barostat.

Returned metadata does not record the effective mode. It normalizes and echoes the requested ensemble, thermostat, and barostat. Therefore metadata alone cannot prove which equations were integrated. Until the runtime reports an effective_integrator field or rejects inconsistent combinations, reproducible studies should derive the route from this table and record both requested and effective settings.

The user-facing registry entry for this limitation is in Known Limitations.

Constraints across dynamics modes

Tako converts each atom’s selective-coordinate flags into either a completely fixed-atom constraint or a Cartesian component mask. Forces are adjusted by constraints, and constrained position setters can alter the trial displacement. The nominal degree count subtracts the removed components.

Constraint compatibility is nevertheless mode dependent:

Constraints change statistical mechanics. A fixed atom does not contribute three active kinetic degrees. A partially fixed atom contributes only its free components. Removing forces without changing thermostat DOF overheats the active subspace. A constraint also changes pressure and virial interpretation. The current MTK rejection is safer than silently applying an inconsistent cell/constraint treatment, but it means a constrained slab cannot use Tako’s MTK pressure coupling.

Constraint validation should include checking that fixed coordinates remain fixed to numerical precision, that the reported temperature uses an understood denominator, that no unintended impulse appears at the first step, and that structural observables are normalized over active atoms where appropriate. If a restraint rather than a rigid constraint is scientifically intended, an explicit bias potential would be needed; Tako’s selective flags are hard Cartesian constraints, not harmonic restraints.

Center-of-mass handling by mode

Center-of-mass removal is often used to prevent numerical drift, but it is itself a constraint on phase space. It changes the kinetic degree count and can distort dynamics if applied carelessly to small or externally driven systems.

Tako always removes COM velocity during initialization. Thereafter:

Consequently, the same checkbox does not mean the same equations in every mode. It cannot preserve initial COM motion when false, and it cannot enforce continued removal in NHC or MTK when true.

The reported temperature denominator does not subtract three COM degrees. Thus continuous COM fixing removes kinetic components without matching the printed degree count. For large systems the relative difference is small; for a three-atom molecule it is substantial. A scientifically careful analysis should recompute kinetic temperature with the intended active degrees, but Tako does not expose velocities, so this correction cannot be performed exactly from public trajectory artifacts.

Randomness and seeds

Reproducible stochastic MD requires distinguishing deterministic replay from statistical replication. Tako uses one submitted integer seed in two separate generators. A bridge-level generator creates the initial velocity draw. A different ase.rs generator, initialized independently with the same integer, drives Bussi, Langevin, or Andersen steps. The sequences are deterministic for a given build but are not one continuous random stream.

Public normalization forces mdSeed to a positive integer. The lower-level initializer has a special nonzero replacement for seed zero, but ordinary UI or script settings will normalize zero to the default. For deterministic Verlet, Berendsen, NHC, and MTK, the seed controls initial velocities but not later random forces because those methods are deterministic. For Bussi, Langevin, and Andersen, it controls both initialization and the thermostat’s independent stochastic stream.

One repeated seed is valuable for regression testing, timestep comparison, or isolating a changed model. It does not estimate uncertainty. Statistical work should use multiple independent seeds, identical equilibrated protocol, and a prespecified combination rule. Report between-replica dispersion and within-replica autocorrelation. Trajectories that share an initial equilibrated configuration have correlated early segments even when their thermostat seeds differ.

Tako does not store the RNG state in its public result. ase.rs can put Bussi RNG state into a restart object, but the bridge discards that object. Therefore a nominal continuation using the same seed restarts the random sequence and regenerates velocities; it is not a continuation.

Frames, scalar traces, and time indexing

The ase.rs loop numbers completed integration steps from 1 through $S$. The timestamp is

No step-zero frame is emitted. A frame is due when step % frame_interval == 0. Tako asks the event sink for frames and disables ase.rs owned-frame storage to avoid duplication. Every due frame is streamed with coordinates, cell, potential energy, kinetic energy, total energy, and temperature.

The final scalar observables are always evaluated even if the final step is not frame-aligned. In that case the event contains a scalar trace with no coordinate frame. This produces an important divergence:

• md_trace includes every frame-aligned trace and the final scalar trace.
• frames, trajectory_xyz, md.traj, and the frame-derived energy_trace.json include only frame-aligned coordinate frames.
• If $S$ is not divisible by the frame interval, the persisted coordinate trajectory and energy trace omit the final state while final_energy_*, final_temperature_k, and the last md_trace entry describe it.

Potential energy is not necessarily evaluated on every step. On steps that do not require output, final reporting, or logging, the integrator may request forces only and reuse cached force information where possible. Observables are evaluated on due/final steps. This is efficient, but it means the scalar trace is a sampled record, not a hidden full-resolution energy series.

The frame interval is clamped to the total step count. Thus an interval larger than the run produces a single final frame. The step count is clamped to 1–10,000 in the current browser bridge. A long physical trajectory may therefore require multiple calls, but because restart state and velocities are not exposed, separate calls are statistically discontinuous and cannot be concatenated as one exact trajectory.

The user-facing registry entry for this limitation is in Known Limitations.

Result schema and unit caveats

The molecular-dynamics result contains method provenance, requested settings, step count, final potential and kinetic energies in Hartree and eV, final temperature, streamed frames, scalar MD trace, an XYZ-like trajectory string, and final atomic positions. It does not contain final velocities, a top-level final cell, pressure, stress, volume, density, enthalpy, conserved extended energy, thermostat work, barostat work, active DOF, or effective routing.

Each frame includes step, time_fs, potential energy in Hartree and eV, kinetic and total energy in eV, temperature K, cell vectors, and atom positions. Because the frame shares a generic optimization schema, force-maximum and maximum-step fields are filled with zero; those zeros are placeholders, not measured MD force or displacement diagnostics.

Each md_trace entry includes potential, kinetic, and total energy in Hartree and eV, temperature, time, step, and a cell matrix. There is a unit inconsistency: frame cell vectors are passed through output_cell() and converted from Bohr to Å, but trace cells are emitted directly from ase.rs internal matrices and remain in Bohr without a unit label. Until corrected, do not combine md_trace.cell numerically with frame-cell values under an assumed common unit.

Requested settings are normalized for display but do not include the effective runtime mode. Pressure traces cannot be reconstructed without stress or velocities. Approximate volume can be computed from frame cells, but not safely from the trace without applying the undocumented Bohr conversion. Temperature can be read, but its DOF convention has the constraint/COM limitations described earlier.

The user-facing registry entry for this limitation is in Known Limitations. The exhaustive field list is generated in tako.molecularDynamics Settings; Script handling is in Molecular Dynamics with Tako Script.

Restart and continuation limitations

The underlying ase.rs result constructs a restart containing the final atoms, velocities, cell, step, and, where relevant, Bussi random/transferred-energy state, NHC chain variables, or MTK thermostat/barostat/cell momentum. It can serialize versioned restart files.

The Tako WASM bridge discards this restart object. The public result contains final positions but not velocities or extended variables. There is no restart input in the public MD settings. Cancellation exposes no checkpoint. A later MD call always generates fresh velocities from its temperature and seed.

Therefore:

• two calls cannot be concatenated as one continuous trajectory;
• an equilibration call cannot transfer its exact phase-space state into production;
• stochastic RNG continuation is absent;
• NHC and MTK extended variables reset;
• changing only the step count and reusing the final structure is a new simulation, not a continuation;
• the 10,000-step bridge cap is a hard practical limitation on continuous public trajectories.

Trajectory coordinates may be concatenated for visualization only when the discontinuity is made explicit. They must not be analyzed as a continuous time series for diffusion, correlation, or rates.

The user-facing registry entry for this limitation is in Known Limitations.

What Tako MD can and cannot establish

Tako’s MD surface can establish that a specified calculator and effective integrator produce a particular finite trajectory from Tako’s deterministic initialization. It can support local structural relaxation observations, short-time motion, preliminary canonical structural sampling under suitable thermostats, fixed-cell stability tests, and exploratory cell response when stress is reliable. Its browser execution and explicit artifacts make small studies inspectable.

It cannot by itself establish experimental stability, equilibrium populations across unsampled barriers, macroscopic transport, reaction rates, or thermodynamic properties without the convergence and validation layers developed in Molecular Dynamics: Statistics and Validation. A trajectory is not experimental time: the force model may omit electronic transitions, quantum nuclei, explicit environment, defects, or rare mechanisms.

Current implementation limitations further narrow claims. Requested and effective controls may disagree silently. Initial velocities use inconsistent COM/DOF normalization. NHC counts $3N$ despite constraints. MTK rejects constraints. COM control is mode dependent. Partial periodicity is lost. Pressure and stress traces are absent. Trace-cell units differ from frame-cell units. Final coordinates can be omitted from the frame sequence. Velocities, extended energy, and restart are not public. Continuous runs are capped at 10,000 steps.

These constraints do not make all MD results invalid. They determine the evidence level. A short constrained molecule run with Bussi may still show whether a local geometry remains intact, provided the temperature convention is acknowledged. An NVE timestep series can test force consistency. A small NPT Berendsen run can screen gross cell response. But the language must remain proportional: “screening,” “demonstration,” “within this model and accessible time,” and “requires validation” are often the correct qualifiers.

Future runtime improvements that would materially strengthen science include rejecting contradictory control combinations, returning an explicit effective mode, degree count, pressure/stress/volume and extended-energy traces, using constraint- and COM-consistent initialization, preserving partial periodicity, ensuring the final coordinate frame is stored, exporting velocities and restart state, and supporting staged continuation without phase-space reset.

Primary references and further reading

The following sources support the algorithms and statistical-mechanical framework discussed here:

1. L. Verlet, “Computer ‘Experiments’ on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules,” Physical Review 159, 98–103 (1967), DOI: 10.1103/PhysRev.159.98.
2. W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson, “A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules,” Journal of Chemical Physics 76, 637–649 (1982), DOI: 10.1063/1.442716.
3. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, “Molecular dynamics with coupling to an external bath,” Journal of Chemical Physics 81, 3684–3690 (1984), DOI: 10.1063/1.448118.
4. G. Bussi, D. Donadio, and M. Parrinello, “Canonical sampling through velocity rescaling,” Journal of Chemical Physics 126, 014101 (2007), DOI: 10.1063/1.2408420.
5. H. C. Andersen, “Molecular dynamics simulations at constant pressure and/or temperature,” Journal of Chemical Physics 72, 2384–2393 (1980), DOI: 10.1063/1.439486.
6. S. Nosé, “A unified formulation of the constant temperature molecular dynamics methods,” Journal of Chemical Physics 81, 511–519 (1984), DOI: 10.1063/1.447334.
7. W. G. Hoover, “Canonical dynamics: Equilibrium phase-space distributions,” Physical Review A 31, 1695–1697 (1985), DOI: 10.1103/PhysRevA.31.1695.
8. G. J. Martyna, M. L. Klein, and M. Tuckerman, “Nosé–Hoover chains: The canonical ensemble via continuous dynamics,” Journal of Chemical Physics 97, 2635–2643 (1992), DOI: 10.1063/1.463940.
9. G. J. Martyna, D. J. Tobias, and M. L. Klein, “Constant pressure molecular dynamics algorithms,” Journal of Chemical Physics 101, 4177–4189 (1994), DOI: 10.1063/1.467468.
10. M. E. Tuckerman, J. Alejandre, R. López-Rendón, A. L. Jochim, and G. J. Martyna, “A Liouville-operator derived measure-preserving integrator for molecular dynamics simulations in the isothermal-isobaric ensemble,” Journal of Physics A 39, 5629–5651 (2006), DOI: 10.1088/0305-4470/39/19/S18.
11. E. Vanden-Eijnden and G. Ciccotti, “Second-order integrators for Langevin equations with holonomic constraints,” Chemical Physics Letters 429, 310–316 (2006), DOI: 10.1016/j.cplett.2006.07.086.