Molecular Dynamics with Tako Script

tako.molecularDynamics runs the same WebAssembly molecular-dynamics operation as MD Setup, but it makes the complete normalized input record available in source. Use it when settings must be generated, compared across timestep or seed series, or preserved with downstream analysis. The classical equations, thermostats, barostats, sampling tests, and implementation limitations are documented in Molecular Dynamics. The visible dialog belongs in Run Molecular Dynamics.

Call the operation

Pass one structure and one flat options object. Keep the calculator and every MD control in that object; names use the md prefix.

const structure = await tako.input("input/liquid.extxyz");
const calculator = tako.calculator.mlip({ levelOfTheory: "nequix" });

const result = await tako.molecularDynamics(structure, {
  calculator,
  properties: ["energy", "forces"],
  mdSteps: 2000,
  mdTimestepFs: 0.5,
  mdTemperatureK: 300,
  mdFrameInterval: 10,
  mdEnsemble: "nvt",
  mdThermostat: "bussi",
  mdBarostat: "none",
  mdTautFs: 100,
  mdTaupFs: 1000,
  mdPressureEvPerAngstrom3: 0,
  mdCompressibilityAngstrom3PerEv: 0.01,
  mdFrictionPerFs: 0.01,
  mdAndersenProbability: 0.01,
  mdTchain: 3,
  mdPchain: 3,
  mdTloop: 1,
  mdPloop: 1,
  mdMtkMask: [true, true, true],
  mdFixCom: true,
  mdSeed: 42,
});

await tako.output("output/md-result.json", result);

The normalized maximum is 10,000 steps per call. A second call is a new trajectory: it assigns new velocities from its seed and does not consume the first call’s final velocities, thermostat variables, cell momentum, or random-number state. The public API therefore cannot extend a trajectory continuously.

Select a coherent runtime mode

The three selection fields are normalized independently, but the runtime resolves them by precedence. Returned result.settings contains the normalized requested labels, not the resolved integrator. Use only combinations whose labels communicate the effective mode:

npt with thermostat none becomes fixed-cell velocity Verlet. npt with Nose–Hoover and barostat none becomes isotropic MTK. Other ambiguous NPT combinations generally become Berendsen NPT or an MTK mode. A barostat under NVT is ignored. These are current routing rules, not alternate ensemble definitions.

Validate before dispatch

Validation confirms calculator capability and normalized property names; it does not prove that a thermostat/barostat combination resolves as intended. Add an explicit guard in reusable scripts:

const operation = "molecularDynamics";
const properties = ["energy", "forces"];
const check = calculator.validate(operation, { properties });
if (!check.ok) {
  throw new Error(check.reasons.join("; "));
}

const mode = ["nvt", "bussi", "none"];
const allowed = new Set([
  "nve|none|none",
  "nvt|berendsen|none",
  "nvt|bussi|none",
  "nvt|langevin|none",
  "nvt|noseHoover|none",
  "nvt|andersen|none",
  "npt|berendsen|berendsen",
  "npt|noseHoover|mtkIsotropic",
  "npt|noseHoover|mtkFull",
  "npt|noseHoover|mtkMasked",
]);
if (!allowed.has(mode.join("|"))) {
  throw new Error(`Ambiguous MD routing: ${mode.join(" / ")}`);
}

NPT additionally requires a periodic cell and stress from the calculator. MTK modes reject constraints. Tako currently promotes any structure carrying a cell to full three-dimensional periodicity for calculation input, so do not use public NPT for a partially periodic slab or chain.

Understand initialization

Every call creates mass-weighted Gaussian velocity components from mdSeed, removes total center-of-mass velocity, then rescales kinetic energy using $3N$ components and mdTemperatureK. Input velocities are not continued. The initializer does not replace $3N$ by the constrained degrees of freedom, even though reported temperature uses the constrained count. Constrained systems can therefore show an immediate temperature mismatch.

mdFixCom: false does not preserve initial translation; initialization always removes it. During propagation the flag is applied by Berendsen, Bussi, Langevin, and Andersen modes. Nose–Hoover-chain and MTK propagation do not consult it.

For deterministic modes, the seed affects only the initial state. For Bussi, Langevin, and Andersen it also seeds a second runtime random-number sequence. The same seed is reproducible only with identical atom ordering, masses, constraints, calculator, settings, and runtime build.

Read the result boundary

The promise resolves to one MD object. The public top-level fields are:

There is no step-zero frame. If mdSteps % mdFrameInterval !== 0, the final step appears in md_trace but not in frames or trajectory_xyz. A frame contains step, time, positions, cell, potential/kinetic/total energy, and temperature. Frame cell vectors and positions are in ångström. The inherited fmax and maxstep values are zeros and have no MD meaning.

Each md_trace entry contains step, time in fs, potential/kinetic/total energies in Hartree and eV, temperature in K, and cell. The current trace cell is emitted in the engine’s internal Bohr convention, whereas frame cells use ångström (Known Limitations). Convert trace cells before computing volume:

const BOHR_TO_ANGSTROM = 0.529177210903;
const traceCellAngstrom = result.md_trace.map((entry) => ({
  step: entry.step,
  cell: entry.cell?.map((row) =>
    row.map((value) => value * BOHR_TO_ANGSTROM)
  ),
}));

The result does not provide velocities, momenta, pressure, stress, density, enthalpy, conserved extended energy, thermostat/barostat energy, effective resolved mode, degrees of freedom, or restart data. Derived claims must remain inside that boundary.

Run a timestep series

Independent calls cannot share precisely the same initial velocities because no velocity input is public, but using one structure and seed provides a deterministic comparison within the same runtime. Keep duration short and compare the total-energy trace:

const timesteps = [0.5, 0.25, 0.125];
const probes = [];

for (const dt of timesteps) {
  const probe = await tako.molecularDynamics(structure, {
    calculator,
    properties: ["energy", "forces"],
    mdSteps: Math.round(100 / dt),
    mdTimestepFs: dt,
    mdTemperatureK: 300,
    mdFrameInterval: 10,
    mdEnsemble: "nve",
    mdThermostat: "none",
    mdBarostat: "none",
    mdSeed: 42,
  });
  probes.push({ dt, trace: probe.md_trace });
}

await tako.output("output/timestep-series.json", probes);

Compare total-energy drift per unit time and the scientific observable, not only the final potential energy. The smallest timestep is a reference rather than an exact result; model noise and constraint work can remain.

Run independent replicas

Replica seeds quantify initialization and stochastic-thermostat variability. They do not extend total time along one trajectory:

const seeds = [11, 23, 37, 53];
const replicas = [];

for (const mdSeed of seeds) {
  const replica = await tako.molecularDynamics(structure, {
    calculator,
    properties: ["energy", "forces"],
    mdSteps: 5000,
    mdTimestepFs: 0.5,
    mdTemperatureK: 300,
    mdFrameInterval: 20,
    mdEnsemble: "nvt",
    mdThermostat: "bussi",
    mdBarostat: "none",
    mdTautFs: 100,
    mdSeed,
  });
  replicas.push({ mdSeed, result: replica });
}

await tako.output("output/md-replicas.json", replicas);

Discard equilibration and estimate autocorrelation/block uncertainty outside the operation. A saved-frame count is not an independent-sample count.

Handle failures and partial evidence

The promise rejects on invalid NPT periodicity, missing calculator stress, MTK constraints, an empty MTK mask, calculator failure, or integration failure. Cancellation does not yield a restart checkpoint. Streamed frames may have appeared in the workspace before failure, but the rejected call has no successful returned result; treat those frames as diagnostic partial evidence, not a completed trajectory.

Before publication, preserve the input structure, complete options object, calculator checkpoint identity, seed, runtime version, result JSON, and the analysis code that selects equilibration and computes uncertainty. The exhaustive types, units, availability conditions, and normalization rules are generated in tako.molecularDynamics Settings.