Molecular Dynamics Integrators#

Molecular dynamics (MD) simulations propagate atomic positions and velocities forward in time using numerical integrators. ALCHEMI Toolkit-Ops provides GPU-accelerated implementations of standard MD integrators and thermostats via NVIDIA Warp, with full PyTorch autograd support for machine learning applications.

Tip

For most applications, start with velocity_verlet_position_update() for NVE simulations, or langevin_baoab_half_step() for NVT simulations. These integrators are time-reversible, symplectic, and provide excellent energy conservation or temperature control.

Overview of Available Integrators#

ALCHEMI Toolkit-Ops provides integrators for different statistical ensembles:

Integrator	Ensemble	Conservation	Best For
Velocity Verlet	NVE	Energy	Testing, production NVE runs
Langevin (BAOAB)	NVT	Temperature	Canonical sampling, equilibration
Nosé-Hoover Chain	NVT	Temperature	Deterministic thermostat, long runs
NPT (MTK)	NPT	Temperature, Pressure	Constant pressure simulations
NPH (MTK)	NPH	Enthalpy, Pressure	Adiabatic constant pressure
Velocity Rescaling	-	-	Quick equilibration (non-canonical)

All integrators support:

Single-system and batched calculations (via batch_idx; most integrators also support atom_ptr – see individual function docs for details)
Automatic differentiation (positions, velocities, forces)
Both mutating (in-place) and non-mutating (out) variants
Float32 and float64 precision

Quick Start#

Velocity Verlet (NVE)

import warp as wp
from nvalchemiops.dynamics.integrators import (
    velocity_verlet_position_update,
    velocity_verlet_velocity_finalize
)

# Setup
positions = wp.array(pos_np, dtype=wp.vec3d, device="cuda:0")
velocities = wp.array(vel_np, dtype=wp.vec3d, device="cuda:0")
forces = wp.array(force_np, dtype=wp.vec3d, device="cuda:0")
masses = wp.array(mass_np, dtype=wp.float64, device="cuda:0")
dt = wp.array([0.001], dtype=wp.float64, device="cuda:0")

# MD loop
for step in range(num_steps):
    # Step 1: Update positions and half-step velocities
    velocity_verlet_position_update(positions, velocities, forces, masses, dt)

    # Step 2: Recalculate forces at new positions
    forces = compute_forces(positions)  # User-defined

    # Step 3: Finalize velocity update
    velocity_verlet_velocity_finalize(velocities, forces, masses, dt)

Langevin (NVT)

import warp as wp
from nvalchemiops.dynamics.integrators import (
    langevin_baoab_half_step,
    langevin_baoab_finalize
)

# Setup (NVT parameters)
temperature = wp.array([1.0], dtype=wp.float64, device="cuda:0")  # kT in energy units
friction = wp.array([1.0], dtype=wp.float64, device="cuda:0")  # friction coefficient

# MD loop
for step in range(num_steps):
    # Step 1: BAOAB half-step (B-A-O-A)
    langevin_baoab_half_step(
        positions, velocities, forces, masses, dt,
        temperature, friction, random_seed=step
    )

    # Step 2: Recalculate forces
    forces = compute_forces(positions)

    # Step 3: Final B step
    langevin_baoab_finalize(velocities, forces, masses, dt)

Velocity Initialization

import warp as wp
from nvalchemiops.dynamics.utils import (
    initialize_velocities,
    compute_kinetic_energy,
    compute_temperature
)

# Target temperature (k_B*T in energy units)
temperature = wp.array([1.0], dtype=wp.float64, device="cuda:0")

# Scratch arrays for COM removal (required when remove_com=True)
total_momentum = wp.zeros(1, dtype=wp.vec3d, device="cuda:0")
total_mass = wp.zeros(1, dtype=wp.float64, device="cuda:0")
com_velocities = wp.zeros(1, dtype=wp.vec3d, device="cuda:0")

# Initialize velocities from Maxwell-Boltzmann distribution
initialize_velocities(
    velocities, masses, temperature,
    total_momentum, total_mass, com_velocities,
    random_seed=42,
    remove_com=True  # Remove center-of-mass motion
)

# Verify temperature
ke = compute_kinetic_energy(velocities, masses)
T_out = wp.zeros(1, dtype=wp.float64, device="cuda:0")
num_atoms_per_system = wp.array([100], dtype=wp.int32, device="cuda:0")
compute_temperature(ke, T_out, num_atoms_per_system)
print(f"Target: {temperature.numpy()[0]}, Actual: {T_out.numpy()[0]}")

FIRE2 Optimization

import warp as wp
from nvalchemiops.dynamics.optimizers import fire2_step

# One system of num_atoms atoms (FIRE2 requires batch_idx)
positions = wp.array(pos_np, dtype=wp.vec3d, device="cuda:0")
velocities = wp.zeros(num_atoms, dtype=wp.vec3d, device="cuda:0")
forces = wp.zeros(num_atoms, dtype=wp.vec3d, device="cuda:0")
batch_idx = wp.zeros(num_atoms, dtype=wp.int32, device="cuda:0")

# Per-system adaptive state and scratch (all single-entry arrays)
alpha = wp.array([0.09], dtype=wp.float64, device="cuda:0")
dt = wp.array([0.005], dtype=wp.float64, device="cuda:0")
nsteps_inc = wp.array([0], dtype=wp.int32, device="cuda:0")
vf = wp.array([0.0], dtype=wp.float64, device="cuda:0")
v_sumsq = wp.array([0.0], dtype=wp.float64, device="cuda:0")
f_sumsq = wp.array([0.0], dtype=wp.float64, device="cuda:0")
max_norm = wp.array([0.0], dtype=wp.float64, device="cuda:0")
max_force = wp.zeros(1, dtype=wp.float64, device="cuda:0")

for _ in range(max_steps):
    forces = compute_forces(positions)  # User-defined

    max_force.zero_()
    wp.launch(
        max_force_component_kernel,
        dim=forces.shape[0],
        inputs=[forces, max_force],
        device=forces.device,
    )
    if wp.utils.array_sum(max_force) < force_tol:
        break

    fire2_step(
        positions, velocities, forces,
        batch_idx=batch_idx,
        alpha=alpha, dt=dt, nsteps_inc=nsteps_inc,
        vf=vf, v_sumsq=v_sumsq, f_sumsq=f_sumsq, max_norm=max_norm,
    )

You own the loop and the convergence test; fire2_step adapts the timestep and mixing and applies one bounded step per call.

The optimization examples read the largest per-atom force magnitude without converting the force array to NumPy. Define this reduction once; the max stays on the device and only the final scalar returns to Python:

@wp.kernel
def max_force_component_kernel(
    forces: wp.array(dtype=wp.vec3d),
    max_force: wp.array(dtype=wp.float64),
):
    i = wp.tid()
    wp.atomic_max(max_force, 0, wp.length(forces[i]))

Allocate max_force once. Before each convergence check, zero it, launch the kernel, and read the result with wp.utils.array_sum(max_force).

Batch Mode: Simulating Multiple Systems#

All integrators support three execution modes for efficient multi-system simulations:

Single System Mode (Default)#

Standard mode for simulating one system:

dt = wp.array([0.001], dtype=wp.float64, device="cuda:0")
velocity_verlet_position_update(positions, velocities, forces, masses, dt)

Batch Mode with `batch_idx` (Atomic Operations)#

For systems with varying atom counts, where each atom is tagged with its system ID. Launches with dim=num_atoms_total (one thread per atom):

# 3 systems: 30, 40, and 30 atoms
batch_idx = wp.array([0]*30 + [1]*40 + [2]*30, dtype=wp.int32, device="cuda:0")

# Per-system timesteps
dt = wp.array([0.001, 0.002, 0.0015], dtype=wp.float64, device="cuda:0")

velocity_verlet_position_update(
    positions, velocities, forces, masses, dt, batch_idx=batch_idx
)

Use batch_idx when:

Systems have similar sizes
You want maximum parallelism (one thread per atom)
Memory access patterns are coalesced

Batch Mode with `atom_ptr` (Sequential Per-System)#

CSR-style pointers defining atom ranges, where each thread processes one complete system. Launches with dim=num_systems:

# Same 3 systems as above: [0:30], [30:70], [70:100]
atom_ptr = wp.array([0, 30, 70, 100], dtype=wp.int32, device="cuda:0")

# Per-system timesteps
dt = wp.array([0.001, 0.002, 0.0015], dtype=wp.float64, device="cuda:0")

velocity_verlet_position_update(
    positions, velocities, forces, masses, dt, atom_ptr=atom_ptr
)

Use atom_ptr when:

Systems have very different sizes
You need per-system operations (reductions, thermostat chains)
Each system needs independent sequential processing

Integrator Details#

Velocity Verlet (NVE)#

The velocity Verlet algorithm is a second-order symplectic integrator that exactly conserves energy in the absence of numerical error:

\[\begin{split} \begin{aligned} \mathbf{r}(t + \Delta t) &= \mathbf{r}(t) + \mathbf{v}(t) \Delta t + \frac{1}{2} \mathbf{a}(t) \Delta t^2 \\ \mathbf{v}(t + \Delta t) &= \mathbf{v}(t) + \frac{1}{2}[\mathbf{a}(t) + \mathbf{a}(t + \Delta t)] \Delta t \end{aligned} \end{split}\]

Key Properties:

Time-reversible and symplectic
Excellent long-term energy conservation
Requires two force evaluations per step (before and after position update)

References:

Swope et al. (1982). J. Chem. Phys. 76, 637

Langevin Dynamics (NVT)#

Langevin dynamics adds friction and random forces to maintain constant temperature. We implement the BAOAB splitting scheme for optimal configurational sampling:

\[\begin{split} B: \mathbf{v} \leftarrow \mathbf{v} + \frac{\Delta t}{2m}\mathbf{F} \\ A: \mathbf{r} \leftarrow \mathbf{r} + \frac{\Delta t}{2}\mathbf{v} \\ O: \mathbf{v} \leftarrow e^{-\gamma \Delta t} \mathbf{v} + \sqrt{\frac{k_B T (1 - e^{-2\gamma \Delta t})}{m}} \boldsymbol{\xi} \\ A: \mathbf{r} \leftarrow \mathbf{r} + \frac{\Delta t}{2}\mathbf{v} \\ B: \mathbf{v} \leftarrow \mathbf{v} + \frac{\Delta t}{2m}\mathbf{F} \end{split}\]

where \(\gamma\) is the friction coefficient and \(\boldsymbol{\xi} \sim \mathcal{N}(0, 1)\).

Key Properties:

Maintains canonical (NVT) ensemble
Friction coefficient \(\gamma\) controls thermalization rate
Stochastic (requires random seed)
BAOAB splitting provides optimal sampling

References:

Leimkuhler & Matthews (2013). J. Chem. Phys. 138, 174102

Nosé-Hoover Chain (NVT)#

Deterministic thermostat using extended phase space with chain of thermostats:

\[\begin{split} \begin{aligned} \dot{\mathbf{r}}_i &= \mathbf{v}_i \\ \dot{\mathbf{v}}_i &= \frac{\mathbf{F}_i}{m_i} - \dot{\eta}_1 \mathbf{v}_i \\ \dot{\eta}_1 &= \frac{2 \cdot KE - N_{\text{DOF}} k_B T}{Q_1} \\ \dot{\eta}_k &= \frac{Q_{k-1} \dot{\eta}_{k-1}^2 - k_B T}{Q_k} \quad (k > 1) \end{aligned} \end{split}\]

Key Properties:

Deterministic (no random forces)
Rigorously canonical ensemble
Chain length typically 3-5 for good ergodicity
Requires thermostat masses \(Q_k\) (computed via nhc_compute_masses)

References:

Martyna, Tobias, Klein (1994). J Chem Phys, 101, 4177

Velocity Rescaling#

Simple rescaling of velocities to match target temperature:

\[ \mathbf{v}_i \leftarrow \mathbf{v}_i \cdot \sqrt{\frac{T_{\text{target}}}{T_{\text{current}}}} \]

Key Properties:

Very fast equilibration
Does NOT produce canonical ensemble
Useful for initial equilibration before switching to proper thermostat
Can cause artifacts if used for production runs

NPT (Isothermal-Isobaric)#

Constant temperature and pressure simulations using Martyna-Tobias-Klein (MTK) equations. The convenience run_npt_step applies particle thermostat friction inside its NPT velocity half-steps; it is not a drop-in primitive for integrator splits that apply particle Nosé-Hoover velocity scaling separately.

from nvalchemiops.dynamics.integrators import run_npt_step

# Run a complete NPT step
run_npt_step(
    positions, velocities, forces, masses, dt,
    cell, cell_velocities,
    target_temperature, target_pressure,
    nhc_positions, nhc_velocities, nhc_masses,
    barostat_mass, dof
)

Key Properties:

Maintains constant temperature and pressure
Supports isotropic (scalar), orthorhombic (3 components), and fully anisotropic (9 components) pressure control
Uses an NPT velocity primitive that includes particle thermostat drag
Use nph_velocity_half_step or another no-thermostat velocity primitive when composing a separate particle NHC Trotter scaling step

NPH (Isenthalpic-Isobaric)#

Constant enthalpy and pressure simulations without thermostat:

from nvalchemiops.dynamics.integrators import run_nph_step

# Run a complete NPH step
run_nph_step(
    positions, velocities, forces, masses, dt,
    cell, cell_velocities,
    target_pressure,
    barostat_mass, dof
)

Key Properties:

Maintains constant pressure without temperature control
Useful for adiabatic simulations at fixed pressure
Supports isotropic and anisotropic pressure modes

Geometry Optimization#

FIRE (Fast Inertial Relaxation Engine)#

Accelerated gradient descent for finding energy minima:

import warp as wp
from nvalchemiops.dynamics.optimizers import fire_step

# FIRE control parameters (per-system arrays)
alpha = wp.array([0.1], dtype=wp.float64, device="cuda:0")
dt = wp.array([0.1], dtype=wp.float64, device="cuda:0")
alpha_start = wp.array([0.1], dtype=wp.float64, device="cuda:0")
f_alpha = wp.array([0.99], dtype=wp.float64, device="cuda:0")
dt_min = wp.array([1e-3], dtype=wp.float64, device="cuda:0")
dt_max = wp.array([1.0], dtype=wp.float64, device="cuda:0")
maxstep = wp.array([0.1], dtype=wp.float64, device="cuda:0")
n_steps_positive = wp.array([0], dtype=wp.int32, device="cuda:0")
n_min = wp.array([5], dtype=wp.int32, device="cuda:0")
f_dec = wp.array([0.5], dtype=wp.float64, device="cuda:0")
f_inc = wp.array([1.1], dtype=wp.float64, device="cuda:0")

# Scratch arrays
uphill_flag = wp.array([0], dtype=wp.int32, device="cuda:0")
vf = wp.array([0.0], dtype=wp.float64, device="cuda:0")
vv = wp.array([0.0], dtype=wp.float64, device="cuda:0")
ff = wp.array([0.0], dtype=wp.float64, device="cuda:0")
max_force = wp.zeros(1, dtype=wp.float64, device="cuda:0")

for step in range(max_steps):
    # Compute forces
    forces = compute_forces(positions)

    # FIRE step (all parameters are arrays)
    fire_step(
        positions, velocities, forces, masses,
        alpha, dt, alpha_start, f_alpha, dt_min, dt_max,
        maxstep, n_steps_positive, n_min, f_dec, f_inc,
        uphill_flag, vf, vv, ff
    )

    # Check convergence: largest per-atom force magnitude below the tolerance
    max_force.zero_()
    wp.launch(
        max_force_component_kernel,
        dim=forces.shape[0],
        inputs=[forces, max_force],
        device=forces.device,
    )
    fmax = wp.utils.array_sum(max_force)
    if fmax < force_tolerance:
        break

Key Properties:

Adaptive timestep and mixing parameter
Much faster than steepest descent
Suitable for local minimization (not global search)

FIRE2 (Improved FIRE)#

Improved FIRE optimizer with adaptive damping and velocity mixing:

from nvalchemiops.dynamics.optimizers import fire2_step, fire2_update

# FIRE2 with Warp arrays
fire2_step(
    positions, velocities, forces,
    batch_idx=batch_idx,  # Required for FIRE2
    alpha=alpha, dt=dt, nsteps_inc=nsteps_inc,
    vf=vf, v_sumsq=v_sumsq, f_sumsq=f_sumsq, max_norm=max_norm
)

# Mix/update only for custom final apply phases
fire2_update(
    velocities, forces, batch_idx,
    alpha, dt, nsteps_inc,
    vf, v_sumsq, f_sumsq, max_norm,
)

PyTorch Interface:

For PyTorch users, FIRE2 has dedicated high-level adapters:

from nvalchemiops.torch import fire2_step_coord, fire2_step_coord_cell

# Coordinate-only optimization
fire2_step_coord(
    positions, velocities, forces, batch_idx,
    alpha, dt, nsteps_inc
)

# Variable-cell optimization (coordinates + cell DOFs)
fire2_step_coord_cell(
    positions, velocities, forces,
    cell, cell_velocities, cell_force,
    batch_idx,
    alpha, dt, nsteps_inc,
    cell_force_scale=1.0,
)

Key Properties:

Uses batch_idx for batched operations (required)
fire2_update performs FIRE2 reduction, adaptive parameter update, and velocity mixing without applying positions
fire2_update(compute_max_norm=False) skips extended-DOF max-norm reduction for custom final apply phases that recompute their own physical displacement norm
Applies coupled affine cell/coordinate updates for variable-cell optimization through fire2_step_coord_cell
Normalizes stress-derived cell forces by the number of atoms in each system; cell_force_scale=1.0 is the default extra multiplier
Improved convergence compared to original FIRE
PyTorch adapters handle tensor conversion automatically

L-BFGS (Limited-memory Quasi-Newton)#

L-BFGS approximates the inverse Hessian from the last few position and force differences to pick a search direction, then steps along it as far as a maxstep trust region allows. It reads no energies; the only per-step input is one force evaluation.

Choosing between FIRE2 and L-BFGS. Both methods take one force evaluation per step. FIRE2 has little optimizer state; L-BFGS stores 2 * history_size history vectors. Their evaluation counts and step costs depend on the system, force model, precision, and optimizer settings. Use the measurements in the benchmark section as a starting point, then compare both methods on your workload.

Prepare the state once. lbfgs_prepare_state allocates, initializes and validates every array and returns an LBFGSState; calling it again is how you restart. A step allocates nothing, so it stays capturable in a CUDA graph. The arrays stay yours — LBFGSState is a plain dataclass with a validate() method, so you can build one from buffers you already have. See the nvalchemiops.dynamics.optimizers.lbfgs module documentation for shapes and the required initial contents.

Every array follows the coordinate dtype: an fp32 state is fp32 throughout, an fp64 state fp64 throughout, and a mixed one is rejected rather than silently half-converted. curvature_eps follows from that, defaulting to 1e-6 for float32 and 1e-10 for float64 — a pair resolvable in fp64 is noise in fp32, so one value cannot serve both.

You own the loop and the stopping rule. Each lbfgs_step consumes exactly one force evaluation: it updates the history, restarts a direction that stops descending, and takes one bounded step. No tolerance, no terminal status — exactly as with FIRE2.

import warp as wp
from nvalchemiops.dynamics.optimizers import lbfgs_prepare_state, lbfgs_step

state = lbfgs_prepare_state(num_atoms, num_systems, device=device)
max_force = wp.zeros(1, dtype=wp.float64, device=device)

for _ in range(max_steps):
    forces = model(positions)

    # Test *before* stepping: these forces describe the positions you have,
    # and after a step they describe the previous point.
    max_force.zero_()
    wp.launch(
        max_force_component_kernel,
        dim=forces.shape[0],
        inputs=[forces, max_force],
        device=forces.device,
    )
    if wp.utils.array_sum(max_force) < force_tol:
        break

    lbfgs_step(
        positions=positions,
        forces=forces,
        state=state,
        batch_idx=batch_idx,
        maxstep=0.2,      # A, largest displacement in one step
    )

The PyTorch binding mutates the same state in place; the JAX binding returns a new one, since JAX arrays are immutable, and both state classes are registered pytrees there.

Two consequences of testing before stepping. A restart direction is the force normalized, so the trust region caps it at exactly maxstep however small the force — stepping a geometry that has already arrived kicks it. And nothing can stall: a short enough step along a descent direction always makes progress, and an uphill direction is replaced by steepest descent.

No energy is read. The step length comes from the trust region, not from an Armijo test comparing energies — which would compare total energies against a direction built from forces, and reject good steps for any model whose forces are not the gradient of its reported energy. The cost is that the energy may rise on a step; the force still converges.

Batching. Systems are identified by a sorted batch_idx and relax independently, each with its own history and trust-region step. Every system steps on every call, so a batch is done when its slowest member is; compact the batch on the host to retire finished systems earlier.

Variable cell. The packed layout interleaves each system’s atoms with its two cell entries, so the LBFGSState must be sized for num_atoms + 2 * num_systems degrees of freedom, and a second LBFGSCellState carries the chart. Align the cell first with nvalchemiops.dynamics.utils.cell_filter.align_cell — the same requirement, helper and six-component convention fire2_step_coord_cell has; preparation checks it. The one L-BFGS-specific rule is that the reference cell is fixed for the whole relaxation, because stored pairs compare cell coordinates across steps; reset the state if you must re-reference. Apply your thresholds to the Cartesian forces and the stress, never to the packed norms. Systems with no atoms are rejected on this path (there is no scale kappa could take), though empty input on the coordinate-only path remains a no-op.

The full contract — alignment, the six-component convention, the fixed reference, topology, and ragged and empty-system behaviour — is stated once in the nvalchemiops.dynamics.optimizers.lbfgs module documentation. This page and both bindings defer to it.

Memory. The history dominates: 2 * history_size vectors of num_dofs each, roughly 180 bytes per degree of freedom at history_size = 6 with float32 coordinates. A history_size between 3 and 7 is typical.

Temperature Control Utilities#

Computing Temperature#

import warp as wp
from nvalchemiops.dynamics.utils import (
    compute_kinetic_energy,
    compute_temperature
)

# Compute kinetic energy
ke = compute_kinetic_energy(velocities, masses)

# Convert to temperature (assumes k_B = 1)
T_out = wp.zeros(1, dtype=wp.float64, device=velocities.device)
num_atoms_per_system = wp.array([100], dtype=wp.int32, device=velocities.device)
compute_temperature(ke, T_out, num_atoms_per_system)

Cell Utilities#

For periodic boundary conditions and variable-cell simulations:

from nvalchemiops.dynamics.utils import (
    compute_cell_volume,
    compute_cell_inverse,
    scale_positions_with_cell,
    wrap_positions_to_cell,
    cartesian_to_fractional,
    fractional_to_cartesian,
)

# Compute cell volume
volume = compute_cell_volume(cell)

# Wrap positions into primary cell
wrap_positions_to_cell(positions, cell)

# Scale positions when cell changes (preserving fractional coordinates)
scale_positions_with_cell(positions, cell_old, cell_new)

Constraint Utilities (SHAKE/RATTLE)#

Holonomic constraints for fixing bond lengths:

import warp as wp
from nvalchemiops.dynamics.utils import shake_constraints, rattle_constraints

max_error = wp.zeros(1, dtype=wp.float64, device=positions.device)

# After position update: correct positions to satisfy constraints
shake_constraints(
    positions, positions_old, masses,
    bond_atom_i, bond_atom_j, bond_lengths_sq,
    max_error, num_iter=10,
)

# After velocity update: correct velocities to satisfy constraints
rattle_constraints(
    positions, velocities, masses,
    bond_atom_i, bond_atom_j,
    max_error, num_iter=10,
)

Common Pitfalls#

Mutating vs Non-Mutating Functions#

# WRONG: Using mutating function incorrectly
new_positions = velocity_verlet_position_update(...)  # Returns None!

# CORRECT: Use non-mutating variant (requires pre-allocated output arrays)
positions_out = wp.zeros_like(positions)
velocities_out = wp.zeros_like(velocities)
positions_out, velocities_out = velocity_verlet_position_update_out(
    positions, velocities, forces, masses, dt,
    positions_out, velocities_out,
)

# CORRECT: Use mutating variant in-place
velocity_verlet_position_update(positions, velocities, forces, masses, dt)
# positions and velocities are now modified

Timestep as Array#

All integrators require timestep as a Warp array, not a scalar:

# WRONG
velocity_verlet_position_update(positions, velocities, forces, masses, 0.001)

# CORRECT
dt = wp.array([0.001], dtype=wp.float64, device="cuda:0")
velocity_verlet_position_update(positions, velocities, forces, masses, dt)

Batch Mode Mutual Exclusivity#

Cannot use both batch_idx and atom_ptr simultaneously:

# WRONG
velocity_verlet_position_update(
    positions, velocities, forces, masses, dt,
    batch_idx=batch_idx,
    atom_ptr=atom_ptr  # Error: provide one or the other
)

# CORRECT
velocity_verlet_position_update(
    positions, velocities, forces, masses, dt,
    batch_idx=batch_idx
)

Molecular Dynamics Integrators#

Overview of Available Integrators#

Quick Start#

Batch Mode: Simulating Multiple Systems#

Single System Mode (Default)#

Batch Mode with batch_idx (Atomic Operations)#

Batch Mode with atom_ptr (Sequential Per-System)#

Integrator Details#

Velocity Verlet (NVE)#

Langevin Dynamics (NVT)#

Nosé-Hoover Chain (NVT)#

Velocity Rescaling#

NPT (Isothermal-Isobaric)#

NPH (Isenthalpic-Isobaric)#

Geometry Optimization#

FIRE (Fast Inertial Relaxation Engine)#

FIRE2 (Improved FIRE)#

L-BFGS (Limited-memory Quasi-Newton)#

Temperature Control Utilities#

Computing Temperature#

Cell Utilities#

Constraint Utilities (SHAKE/RATTLE)#

Common Pitfalls#

Mutating vs Non-Mutating Functions#

Timestep as Array#

Batch Mode Mutual Exclusivity#

Further Reading#

Batch Mode with `batch_idx` (Atomic Operations)#

Batch Mode with `atom_ptr` (Sequential Per-System)#