Integrator Algorithms#

This page documents exactly what each integrator and optimizer in API Reference computes: the update scheme, the governing equations, and the units of each control parameter. Every scheme is a thin nvalchemi wrapper (a torch.library.custom_op() for torch.compile / autograd compatibility) around a GPU kernel in nvalchemiops; the kernels that realise each method are collected in a collapsible panel at the end of every section, so the implemented equations can be inspected directly.

Note

Units. nvalchemi works in atomic-scale units throughout: length in angstrom (\(\mathrm{\AA}\)), energy in electronvolt (\(\mathrm{eV}\)), force in \(\mathrm{eV}/\mathrm{\AA}\), mass in daltons (amu), and time in femtoseconds (\(\mathrm{fs}\)). Temperature is supplied in kelvin and enters the kernels as \(k_B T\) in \(\mathrm{eV}\). Pressure and stress are in \(\mathrm{eV}/\mathrm{\AA}^3\) (positive for compression; no bar/GPa conversion is applied — see Conventions). Thermostat/barostat coupling times \(\tau\) are in \(\mathrm{fs}\) and the Langevin friction \(\gamma\) in \(\mathrm{fs}^{-1}\).

Velocity Verlet (NVE)#

The microcanonical integrator uses velocity Verlet, which is symplectic and time-reversible and therefore conserves the total energy \(H = \mathrm{KE} + \mathrm{PE}\) to within integration error over long trajectories. The step is split around the single force evaluation:

\[\begin{split}\mathbf{r}(t + \Delta t) &= \mathbf{r}(t) + \mathbf{v}(t)\,\Delta t + \tfrac{1}{2}\,\frac{\mathbf{F}(t)}{m}\,\Delta t^{2} \\ \mathbf{v}(t + \tfrac{\Delta t}{2}) &= \mathbf{v}(t) + \tfrac{1}{2}\,\frac{\mathbf{F}(t)}{m}\,\Delta t \\ \mathbf{v}(t + \Delta t) &= \mathbf{v}(t + \tfrac{\Delta t}{2}) + \tfrac{1}{2}\,\frac{\mathbf{F}(t + \Delta t)}{m}\,\Delta t\end{split}\]

The only control parameter is the timestep \(\Delta t\) (fs); a typical value is 0.5–1.0 fs. There is no thermostat, so any energy drift directly measures model and integration error.

Underlying nvalchemiops kernels

nvalchemiops.dynamics.integrators.velocity_verlet_position_update(positions, velocities, forces, masses, dt, batch_idx=None, atom_ptr=None, device=None)[source]#

Perform velocity Verlet position update step (in-place).

Updates positions to \(\mathbf{r}(t + \Delta t)\) and velocities to the half step \(\mathbf{v}(t + \Delta t/2)\). After calling this function, recalculate forces at the new positions, then call velocity_verlet_velocity_finalize().

The kernel evaluates, with \(\mathbf{a}(t) = \mathbf{F}(t)/m\):

\[ \begin{align}\begin{aligned}\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t)\,\Delta t + \tfrac{1}{2}\mathbf{a}(t)\,\Delta t^2\\\mathbf{v}(t + \tfrac{\Delta t}{2}) = \mathbf{v}(t) + \tfrac{1}{2}\mathbf{a}(t)\,\Delta t\end{aligned}\end{align} \]

Parameters:

positions (wp.array(dtype=wp.vec3f or wp.vec3d)) – Atomic positions. Shape (N,). MODIFIED in-place.
velocities (wp.array(dtype=wp.vec3f or wp.vec3d)) – Atomic velocities. Shape (N,). MODIFIED in-place.
forces (wp.array(dtype=wp.vec3f or wp.vec3d)) – Forces on atoms. Shape (N,).
masses (wp.array(dtype=wp.float32 or wp.float64)) – Atomic masses. Shape (N,).
dt (wp.array(dtype=wp.float32 or wp.float64)) – Timestep(s). Shape (1,) for single system, (B,) for batched.
batch_idx (wp.array(dtype=wp.int32), optional) – System index for each atom. For batched mode (atomic operations).
atom_ptr (wp.array(dtype=wp.int32), optional) – CSR-style pointers. Shape (num_systems + 1,). For batched mode (sequential per-system).
device (str, optional) – Warp device. If None, inferred from positions.

Return type:

None

See also

velocity_verlet_velocity_finalize: Complete the velocity update

nvalchemiops.dynamics.integrators.velocity_verlet_velocity_finalize(velocities, forces_new, masses, dt, batch_idx=None, atom_ptr=None, device=None)[source]#

Finalize velocity Verlet velocity update (in-place).

Completes the velocity update using forces evaluated at the new positions, with \(\mathbf{a}(t + \Delta t) = \mathbf{F}_\text{new}/m\):

\[\mathbf{v}(t + \Delta t) = \mathbf{v}(t + \tfrac{\Delta t}{2}) + \tfrac{1}{2}\mathbf{a}(t + \Delta t)\,\Delta t\]

Parameters:

velocities (wp.array(dtype=wp.vec3f or wp.vec3d)) – Half-step velocities. Shape (N,). MODIFIED in-place to full-step.
forces_new (wp.array(dtype=wp.vec3f or wp.vec3d)) – Forces evaluated at new positions \(\mathbf{r}(t + \Delta t)\). Shape (N,).
masses (wp.array(dtype=wp.float32 or wp.float64)) – Atomic masses. Shape (N,).
dt (wp.array(dtype=wp.float32 or wp.float64)) – Timestep(s). Shape (1,) for single, (B,) for batched.
batch_idx (wp.array(dtype=wp.int32), optional) – System index for each atom. For batched mode (atomic operations).
atom_ptr (wp.array(dtype=wp.int32), optional) – CSR-style pointers. Shape (num_systems + 1,). For batched mode (sequential per-system).
device (str, optional) – Warp device. If None, inferred from velocities.

Return type:

None

See also

velocity_verlet_position_update: First step of velocity Verlet

Langevin dynamics (NVT, BAOAB)#

The canonical-ensemble integrator uses the BAOAB splitting of Leimkuhler & Matthews (2012), which gives high configurational-sampling accuracy. Each step applies the sequence B-A-O-A before the force evaluation and a final B after it, where B is a half force kick, A is a half drift, and O is an Ornstein-Uhlenbeck velocity update that acts as an exact thermostat:

\[\begin{split}\text{B:}\quad & \mathbf{v} \leftarrow \mathbf{v} + \tfrac{\Delta t}{2m}\mathbf{F} \\ \text{A:}\quad & \mathbf{r} \leftarrow \mathbf{r} + \tfrac{\Delta t}{2}\mathbf{v} \\ \text{O:}\quad & \mathbf{v} \leftarrow c_1\,\mathbf{v} + \sqrt{(1 - c_1^{2})\,\frac{k_B T}{m}}\;\boldsymbol{\xi}, \qquad c_1 = e^{-\gamma\,\Delta t}\end{split}\]

with \(\boldsymbol{\xi}\) a standard normal draw. Control parameters are the timestep \(\Delta t\) (fs), the target temperature (K, entering as \(k_B T\)), and the friction \(\gamma\) (\(\mathrm{fs}^{-1}\)), which sets the thermostat coupling strength. The stochastic O step is seeded by random_seed for reproducibility.

Underlying nvalchemiops kernels

nvalchemiops.dynamics.integrators.langevin_baoab_half_step(positions, velocities, forces, masses, dt, temperature, friction, random_seed, batch_idx=None, atom_ptr=None, device=None)[source]#

Perform BAOAB Langevin half-step (B-A-O-A sequence) in-place.

This function performs the first four operations of the BAOAB splitting: B (velocity kick), A (drift), O (Ornstein-Uhlenbeck thermostat), A (drift). The kernel evaluates the following updates in order:

\[ \begin{align}\begin{aligned}B:\quad \mathbf{v} \leftarrow \mathbf{v} + \tfrac{\Delta t}{2m}\mathbf{F}\\A:\quad \mathbf{r} \leftarrow \mathbf{r} + \tfrac{\Delta t}{2}\mathbf{v}\\O:\quad \mathbf{v} \leftarrow c_1\,\mathbf{v} + c_2\,\boldsymbol{\xi}\\A:\quad \mathbf{r} \leftarrow \mathbf{r} + \tfrac{\Delta t}{2}\mathbf{v}\end{aligned}\end{align} \]

with damping \(c_1 = e^{-\gamma \Delta t}\), noise amplitude \(c_2 = \sqrt{k_B T\,(1 - c_1^2)/m}\), and \(\boldsymbol{\xi} \sim \mathcal{N}(0, 1)\).

After calling this function, recalculate forces at the new positions, then call langevin_baoab_finalize() to complete the step.

Parameters:

positions (wp.array(dtype=wp.vec3f or wp.vec3d)) – Atomic positions. Shape (N,). MODIFIED in-place.
velocities (wp.array(dtype=wp.vec3f or wp.vec3d)) – Atomic velocities. Shape (N,). MODIFIED in-place.
forces (wp.array(dtype=wp.vec3f or wp.vec3d)) – Forces on atoms. Shape (N,).
masses (wp.array(dtype=wp.float32 or wp.float64)) – Atomic masses. Shape (N,).
dt (wp.array(dtype=wp.float32 or wp.float64)) – Timestep(s). Shape (1,) for single, (B,) for batched.
temperature (wp.array(dtype=wp.float32 or wp.float64)) – Temperature \(k_B T\). Shape (1,) for single, (B,) for batched.
friction (wp.array(dtype=wp.float32 or wp.float64)) – Friction coefficient \(\gamma\). Shape (1,) for single, (B,) for batched.
random_seed (int) – Random seed for stochastic forces.
batch_idx (wp.array(dtype=wp.int32), optional) – System index for each atom. For batched mode (atomic operations).
atom_ptr (wp.array(dtype=wp.int32), optional) – CSR-style pointers. Shape (num_systems + 1,). For batched mode (sequential per-system).
device (str, optional) – Warp device. If None, inferred from positions.

Return type:

None

Example

Single system NVT simulation:

import warp as wp
import numpy as np

# Setup
positions = wp.array(np.random.randn(100, 3), dtype=wp.vec3d, device="cuda:0")
velocities = wp.array(np.random.randn(100, 3), dtype=wp.vec3d, device="cuda:0")
forces = wp.array(np.random.randn(100, 3), dtype=wp.vec3d, device="cuda:0")
masses = wp.array(np.ones(100), dtype=wp.float64, device="cuda:0")

dt = wp.array([0.001], dtype=wp.float64, device="cuda:0")
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")

# BAOAB half-step
langevin_baoab_half_step(
    positions, velocities, forces, masses, dt, temperature, friction,
    random_seed=42
)

Complete BAOAB step:

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)

Batched mode:

# With batch_idx (3 systems)
batch_idx = wp.array([0]*30 + [1]*40 + [2]*30, dtype=wp.int32, device="cuda:0")
dt = wp.array([0.001, 0.002, 0.0015], dtype=wp.float64, device="cuda:0")
temperature = wp.array([1.0, 1.5, 1.2], dtype=wp.float64, device="cuda:0")
friction = wp.array([1.0, 1.0, 1.0], dtype=wp.float64, device="cuda:0")

langevin_baoab_half_step(
    positions, velocities, forces, masses, dt, temperature, friction,
    random_seed=42, batch_idx=batch_idx
)

See also

langevin_baoab_finalize: Complete the BAOAB step

nvalchemiops.dynamics.integrators.langevin_baoab_finalize(velocities, forces_new, masses, dt, batch_idx=None, atom_ptr=None, device=None)[source]#

Finalize BAOAB Langevin step (final B step) in-place.

Completes the BAOAB sequence with the final velocity half-kick using forces \(\mathbf{F}_\text{new}\) calculated at the new positions:

\[B:\quad \mathbf{v} \leftarrow \mathbf{v} + \tfrac{\Delta t}{2m}\mathbf{F}_\text{new}\]

Parameters:

velocities (wp.array(dtype=wp.vec3f or wp.vec3d)) – Atomic velocities. Shape (N,). MODIFIED in-place.
forces_new (wp.array(dtype=wp.vec3f or wp.vec3d)) – Forces at new positions. Shape (N,).
masses (wp.array(dtype=wp.float32 or wp.float64)) – Atomic masses. Shape (N,).
dt (wp.array(dtype=wp.float32 or wp.float64)) – Timestep(s). Shape (1,) for single, (B,) for batched.
batch_idx (wp.array(dtype=wp.int32), optional) – System index for each atom. For batched mode (atomic operations).
atom_ptr (wp.array(dtype=wp.int32), optional) – CSR-style pointers. Shape (num_systems + 1,). For batched mode (sequential per-system).
device (str, optional) – Warp device. If None, inferred from velocities.

Return type:

None

Nose-Hoover chain thermostat (NVT)#

The deterministic, time-reversible alternative for the canonical ensemble is a Nosé-Hoover chain (NHC) using the Martyna-Tobias-Klein equations with Yoshida-Suzuki factorisation. A chain of chain_length fictitious thermostat variables (default 3) is coupled to the physical DOFs; the chain masses are set from the coupling time \(\tau_T\):

\[Q_0 = N_\text{dof}\,k_B T\,\tau_T^{2}, \qquad Q_k = k_B T\,\tau_T^{2}\ \ (k > 0)\]

Each step propagates the chain by a half step (rescaling the particle velocities), performs the half velocity kick and full position drift, then propagates the chain again after the force evaluation. Control parameters are \(\Delta t\) (fs), temperature (K), and \(\tau_T\) (fs) — typically 10–100 \(\times \Delta t\). The scheme is ergodic for coupled systems but can fail to thermalise stiff, near-harmonic modes, where Langevin is preferable.

Underlying nvalchemiops kernels

nvalchemiops.dynamics.integrators.nhc_compute_masses(ndof, target_temp, tau, chain_length, masses, num_systems=1, device=None, dtype=wp.float64)[source]#

Compute Nosé-Hoover chain masses using GPU kernel.

Computes \(Q_k\) values for the Nosé-Hoover chain:

\[\begin{split}Q_0 &= N_\text{dof} \cdot k_B T \cdot \tau^2 \\ Q_k &= k_B T \cdot \tau^2 \quad \text{for } k > 0\end{split}\]

Parameters:

ndof (wp.array(dtype=wp.int32)) – Number of degrees of freedom per system. Shape (1,) for single system, (num_systems,) for batched.
target_temp (wp.array) – Target temperature \(k_B T\) per system. Shape (1,) for single system, (num_systems,) for batched.
tau (wp.array) – Time constant \(\tau\) per system. Shape (1,) for single system, (num_systems,) for batched.
chain_length (int) – Number of thermostats in the chain.
masses (wp.array) – Chain masses output. Caller must pre-allocate. Shape (chain_length,) for single system, (num_systems, chain_length) for batched.
num_systems (int, optional) – Number of systems for batched mode. Default: 1.
device (str, optional) – Warp device. If None, inferred from masses.
dtype (dtype, optional) – Data type for the masses. Default: wp.float64.

Returns:

Chain masses. Shape (chain_length,) for single system, (num_systems, chain_length) for batched.

Return type:

wp.array

nvalchemiops.dynamics.integrators.nhc_thermostat_chain_update(velocities, masses, eta, eta_dot, eta_mass, target_temp, dt, ndof, ke2, total_scale, step_scale, dt_chain, nloops=1, batch_idx=None, num_systems=1, device=None, compute_ke=True)[source]#

Propagate Nosé-Hoover chain and scale velocities (in-place).

Propagates the chain over the Yoshida-Suzuki sub-steps, accumulating the per-step velocity scale factor \(\exp(-\tfrac{1}{2}\,\Delta t_c\,\dot{\eta}_0)\) into total_scale, then rescales every particle velocity by it:

\[\mathbf{v}_i \leftarrow s\,\mathbf{v}_i, \qquad s = \prod_j \exp\!\left(-\tfrac{1}{2}\,\Delta t_{c,j}\,\dot{\eta}_0\right)\]

Uses Yoshida-Suzuki factorization for time-reversible integration. All computations are performed on GPU using Warp kernels.

Parameters:

velocities (wp.array(dtype=wp.vec3f or wp.vec3d)) – Atomic velocities. Shape (N,). MODIFIED in-place.
masses (wp.array(dtype=wp.float32 or wp.float64)) – Atomic masses. Shape (N,).
eta (wp.array(dtype=wp.float64)) – Thermostat chain positions. Non-batched: Shape (chain_length,). Batched: Shape (num_systems, chain_length) as wp.array2d. MODIFIED in-place.
eta_dot (wp.array(dtype=wp.float64)) – Thermostat chain velocities. Same shape as eta. MODIFIED in-place.
eta_mass (wp.array(dtype=wp.float64)) – Thermostat chain masses. Same shape as eta.
target_temp (wp.array(dtype=wp.float64)) – Target temperature \(k_B T\). Shape (1,) or (num_systems,).
dt (wp.array(dtype=wp.float32 or wp.float64)) – Time step. Shape (1,) or (num_systems,).
ndof (wp.array(dtype=wp.float64)) – Degrees of freedom. Shape (1,) or (num_systems,).
ke2 (wp.array) – \(2\,\mathrm{KE} = \sum_i m_i \lVert \mathbf{v}_i \rVert^2\) per system. When compute_ke is True (default) this is a scratch array zeroed and filled internally from velocities/masses before use. When compute_ke is False it is an INPUT supplying a caller-precomputed \(2\,\mathrm{KE}\); it is not zeroed. Either way the chain forward sweep rescales it in place. Shape (1,) for single system, (num_systems,) for batched.
total_scale (wp.array) – Scratch array for accumulated velocity scale factor. Must be initialized to ones by caller (wp.ones). Shape (1,) for single system, (num_systems,) for batched.
step_scale (wp.array) – Scratch array for per-step velocity scale factor. Shape (1,) for single system, (num_systems,) for batched.
dt_chain (wp.array) – Scratch array for weighted time steps. Shape (1,) for single system, (num_systems,) for batched.
nloops (int, optional) – Number of Yoshida-Suzuki integration sub-steps. Default: 1. Use nloops=3 or 5 for higher accuracy.
batch_idx (wp.array(dtype=wp.int32), optional) – System index for each atom. Required for batched mode.
num_systems (int, optional) – Number of systems for batched mode. Default: 1.
device (str, optional) – Warp device. If None, inferred from velocities.
compute_ke (bool, optional) – If True (default), recompute \(2\,\mathrm{KE}\) from velocities/masses internally (byte-identical to legacy behavior). If False, use the caller-supplied value already in ke2 instead of recomputing it. See nhc_compute_2ke for the matching reduction helper.

Return type:

None

nvalchemiops.dynamics.integrators.nhc_velocity_half_step(velocities, forces, masses, dt, batch_idx=None, atom_ptr=None, device=None)[source]#

Half-step velocity update (in-place).

\[\mathbf{v}_i \leftarrow \mathbf{v}_i + \tfrac{1}{2}\,\frac{\mathbf{F}_i}{m_i}\,\Delta t\]

Parameters:

velocities (wp.array(dtype=wp.vec3f or wp.vec3d)) – Atomic velocities. Shape (N,). MODIFIED in-place.
forces (wp.array(dtype=wp.vec3f or wp.vec3d)) – Forces on atoms. Shape (N,).
masses (wp.array(dtype=wp.float32 or wp.float64)) – Atomic masses. Shape (N,).
dt (wp.array(dtype=wp.float32 or wp.float64)) – Time step. Shape (1,) or (B,).
batch_idx (wp.array(dtype=wp.int32), optional) – System index for each atom. For batched mode (atomic operations).
atom_ptr (wp.array(dtype=wp.int32), optional) – CSR-style pointers. Shape (num_systems + 1,). For batched mode (sequential per-system).
device (str, optional) – Warp device. If None, inferred from velocities.

Return type:

None

nvalchemiops.dynamics.integrators.nhc_position_update(positions, velocities, dt, batch_idx=None, atom_ptr=None, device=None)[source]#

Full-step position update (in-place).

\[\mathbf{r}_i \leftarrow \mathbf{r}_i + \mathbf{v}_i\,\Delta t\]

Parameters:

positions (wp.array(dtype=wp.vec3f or wp.vec3d)) – Atomic positions. Shape (N,). MODIFIED in-place.
velocities (wp.array(dtype=wp.vec3f or wp.vec3d)) – Atomic velocities. Shape (N,).
dt (wp.array(dtype=wp.float32 or wp.float64)) – Time step. Shape (1,) or (B,).
batch_idx (wp.array(dtype=wp.int32), optional) – System index for each atom. For batched mode (atomic operations).
atom_ptr (wp.array(dtype=wp.int32), optional) – CSR-style pointers. Shape (num_systems + 1,). For batched mode (sequential per-system).
device (str, optional) – Warp device. If None, inferred from positions.

Return type:

None

nvalchemiops.dynamics.integrators.nhc_compute_chain_energy(eta, eta_dot, eta_mass, target_temp, ndof, ke_chain, pe_chain, batch_idx=None, num_systems=1, device=None)[source]#

Compute Nosé-Hoover chain kinetic and potential energy.

For conservation checks, the extended system Hamiltonian is:

\[H_\text{ext} = \text{KE}_\text{particles} + \text{PE} + \text{KE}_\text{chain} + \text{PE}_\text{chain}\]

where:

\[ \begin{align}\begin{aligned}\text{KE}_\text{chain} = \sum_k \tfrac{1}{2} Q_k \dot{\eta}_k^2\\\text{PE}_\text{chain} = N_\text{dof} \cdot kT \cdot \eta_0 + kT \sum_{k>0} \eta_k\end{aligned}\end{align} \]

Parameters:

eta (wp.array(dtype=wp.float64)) – Thermostat chain positions. Shape (M,) or (B, M).
eta_dot (wp.array(dtype=wp.float64)) – Thermostat chain velocities. Shape (M,) or (B, M).
eta_mass (wp.array(dtype=wp.float64)) – Thermostat chain masses. Shape (M,) or (B, M).
target_temp (wp.array(dtype=wp.float64)) – Target temperature \(k_B T\). Shape (1,) or (B,).
ndof (wp.array(dtype=wp.float64)) – Degrees of freedom. Shape (1,) or (B,).
ke_chain (wp.array) – Output kinetic energy of the chain. Shape (1,) for single system, (num_systems,) for batched.
pe_chain (wp.array) – Output potential energy of the chain. Shape (1,) for single system, (num_systems,) for batched.
batch_idx (wp.array(dtype=wp.int32), optional) – Not used directly, but included for API consistency.
num_systems (int, optional) – Number of systems for batched mode. Default: 1.
device (str, optional) – Warp device. If None, inferred from eta.

Returns:

(ke_chain, pe_chain) each with shape (1,) or (B,).

Return type:

tuple[wp.array, wp.array]

nvalchemiops.dynamics.integrators.nhc_compute_2ke(velocities, masses, ke2, batch_idx=None, device=None)[source]#

Fill ke2[s] with \(\sum_{i \in s} m_i \lVert \mathbf{v}_i \rVert^2\) (i.e. \(2\,\mathrm{KE}\) per system) — the same reduction that nhc_thermostat_chain_update performs internally, exposed standalone.

Lets a caller compute \(2\,\mathrm{KE}\) separately, optionally post-process it, and feed it back via nhc_thermostat_chain_update(..., compute_ke=False). Using this kernel for the reduction gives machine-precision parity with the value the chain update would otherwise compute itself.

Parameters:

velocities (wp.array(dtype=wp.vec3f or wp.vec3d)) – Atomic velocities. Shape (N,).
masses (wp.array(dtype=wp.float32 or wp.float64)) – Atomic masses. Shape (N,).
ke2 (wp.array) – Output \(2\,\mathrm{KE}\). Shape (1,) for single system, (num_systems,) for batched (its length sets the number of systems). Zeroed internally before accumulation. MODIFIED in-place.
batch_idx (wp.array(dtype=wp.int32), optional) – System index for each atom. Required for batched mode.
device (str, optional) – Warp device. If None, inferred from velocities.

Returns:

The ke2 array, filled with \(2\,\mathrm{KE}\) per system.

Return type:

wp.array

A simpler stochastic-velocity-rescaling thermostat is also available:

nvalchemiops.dynamics.integrators.velocity_rescale(velocities, scale_factor, batch_idx=None, atom_ptr=None, device=None)[source]#

Rescale velocities to achieve target temperature (in-place).

\[\mathbf{v}_i \leftarrow s\,\mathbf{v}_i\]

where \(s\) is the scale factor, typically \(s = \sqrt{T_{\text{target}} / T_{\text{current}}}\).

Parameters:

velocities (wp.array(dtype=wp.vec3f or wp.vec3d)) – Atomic velocities. Shape (N,). MODIFIED in-place.
scale_factor (wp.array) – Scaling factor(s). Shape (1,) for single system, (B,) for batched. Typically \(\sqrt{T_{\text{target}} / T_{\text{current}}}\).
batch_idx (wp.array(dtype=wp.int32), optional) – System index for each atom. For batched mode (atomic operations).
atom_ptr (wp.array(dtype=wp.int32), optional) – CSR-style pointers. Shape (num_systems + 1,). For batched mode (sequential per-system).
device (str, optional) – Warp device. If None, inferred from velocities.

Return type:

None

Example

>>> # Compute current temperature
>>> ke = wp.empty(1, dtype=wp.float64, device=device)
>>> compute_kinetic_energy(velocities, masses, ke)
>>> T_current = wp.empty(1, dtype=wp.float64, device=device)
>>> compute_temperature(ke, T_current, num_atoms=100)
>>>
>>> # Compute rescaling factor
>>> factor = _compute_rescale_factor(T_current, T_target)
>>> scale = wp.array([factor], dtype=wp.float32, device=device)
>>>
>>> # Apply rescaling
>>> velocity_rescale(velocities, scale)

Isothermal-isobaric (NPT)#

NPT couples the cell to a Martyna-Tobias-Klein barostat and the particles and barostat to two independent Nosé-Hoover chains, sampling constant temperature and pressure. The model must provide a stress output. The instantaneous pressure tensor combines the kinetic and virial contributions,

\[\mathbf{P} = \frac{\mathbf{K} + \mathbf{W}}{V},\]

and the barostat inertia is fixed by the coupling time \(\tau_P\),

\[W = (N_f + d)\,k_B T\,\tau_P^{2}.\]

The cell strain rate \(\dot{\varepsilon}\) is advanced by half steps from the pressure imbalance,

\[\dot{\varepsilon}\ \mathrel{+}=\ \frac{\Delta t}{2}\,\frac{V}{W}\, (P_\text{inst} - P_\text{ext}),\]

interleaved with the two NHC updates, the particle velocity/position updates, and the cell update. Control parameters: \(\Delta t\) (fs), temperature (K), target pressure (\(\mathrm{eV}/\mathrm{\AA}^3\)), barostat time \(\tau_P\) (fs), and thermostat time \(\tau_T\) (fs). pressure_coupling selects isotropic, anisotropic (orthorhombic), or triclinic cell fluctuations.

Underlying nvalchemiops kernels

nvalchemiops.dynamics.integrators.compute_pressure_tensor(velocities, masses, virial_tensors, cells, kinetic_tensors, pressure_tensors, volumes, batch_idx=None, device=None, compute_kinetic=True)[source]#

Compute full pressure tensor \(\mathbf{P} = (\mathbf{K} + \mathbf{W}) / V\), with kinetic term \(\mathbf{K}\), virial term \(\mathbf{W}\), and cell volume \(V\).

Parameters:

velocities (wp.array(dtype=wp.vec3f or wp.vec3d)) – Particle velocities. Shape (N,).
masses (wp.array) – Particle masses. Shape (N,).
virial_tensors (wp.array(dtype=vec9f or vec9d)) – Virial tensor from forces (physics sign). Shape (B,).
cells (wp.array(dtype=wp.mat33f or wp.mat33d)) – Cell matrices. Shape (B,).
kinetic_tensors (wp.array(dtype=scalar, ndim=2)) – Kinetic tensor \(\mathbf{K}_s = \sum_{i \in s} m_i (\mathbf{v}_i \otimes \mathbf{v}_i)\), vec9 layout. Shape (B, 9). When compute_kinetic is True (default) this is a scratch array zeroed and filled internally from velocities/masses. When compute_kinetic is False it is an INPUT supplying a caller-precomputed kinetic tensor; it is not zeroed or recomputed.
pressure_tensors (wp.array(dtype=vec9f or vec9d)) – Output pressure tensor. Shape (B,).
volumes (wp.array(dtype=scalar)) – Pre-computed cell volumes. Shape (B,). Caller must pre-compute via compute_cell_volume.
batch_idx (wp.array(dtype=wp.int32), optional) – System index for each atom. If None, assumes single system.
device (str, optional) – Warp device.
compute_kinetic (bool, optional) – If True (default), recompute the kinetic tensor from velocities/masses internally (byte-identical to legacy behavior). If False, use the caller-supplied value already in kinetic_tensors instead of recomputing it. See compute_kinetic_tensor for the matching reduction helper.

Returns:

Pressure tensor. Shape (B,).

Return type:

wp.array(dtype=vec9f or vec9d)

nvalchemiops.dynamics.integrators.compute_scalar_pressure(pressure_tensors, scalar_pressures, device=None)[source]#

Compute scalar pressure from pressure tensor.

\[P = (P_{xx} + P_{yy} + P_{zz}) / 3\]

Parameters:

pressure_tensors (wp.array(dtype=vec9f or vec9d)) – Pressure tensor. Shape (B,).
scalar_pressures (wp.array(dtype=scalar)) – Output scalar pressure. Shape (B,).
device (str, optional) – Warp device.

Returns:

Scalar pressure. Shape (B,).

Return type:

wp.array

nvalchemiops.dynamics.integrators.compute_kinetic_tensor(velocities, masses, kinetic_tensors, batch_idx=None, device=None)[source]#

Fill kinetic_tensors[s] with \(\sum_{i \in s} m_i (\mathbf{v}_i \otimes \mathbf{v}_i)\) (vec9 layout) — the same accumulation compute_pressure_tensor performs internally, exposed standalone.

Lets a caller compute the kinetic tensor separately, optionally post-process it, and feed it back via compute_pressure_tensor(..., compute_kinetic=False). Using this kernel for the reduction gives machine-precision parity with the value compute_pressure_tensor would otherwise compute itself.

Parameters:

velocities (wp.array(dtype=wp.vec3f or wp.vec3d)) – Particle velocities. Shape (N,).
masses (wp.array) – Particle masses. Shape (N,).
kinetic_tensors (wp.array(dtype=scalar, ndim=2)) – Output kinetic tensor. Shape (B, 9). Zeroed internally before accumulation. MODIFIED in-place.
batch_idx (wp.array(dtype=wp.int32), optional) – System index for each atom. If None, assumes single system.
device (str, optional) – Warp device.

Returns:

The kinetic_tensors array, filled with \(\sum_i m_i (\mathbf{v}_i \otimes \mathbf{v}_i)\) per system.

Return type:

wp.array

nvalchemiops.dynamics.integrators.compute_cell_kinetic_energy(cell_velocities, cell_masses, kinetic_energy, cells_inv=None, volumes=None, device=None)[source]#

Compute kinetic energy of cell degrees of freedom.

\(\mathrm{KE}_{\text{cell}} = \frac{1}{2} W \|\dot{\varepsilon}\|_F^2\)

Parameters:

cell_velocities (wp.array(dtype=wp.mat33f or wp.mat33d)) – Strain-rate matrices \(\dot{\varepsilon} = p_g/W\). Shape (B,).
cell_masses (wp.array) – Barostat masses. Shape (B,). Note that this is the per-DOF mass \(W = (d + 1) k_B T \tau_p^2\) consumed by the strain-rate update, not the full barostat mass \(d \cdot W\).
kinetic_energy (wp.array(dtype=scalar)) – Output cell kinetic energy. Shape (B,).
cells_inv (wp.array, optional) –

Deprecated since version 0.3.1: Ignored; cell_velocities is the strain rate eps_dot = p_g/W.
volumes (wp.array, optional) –

Deprecated since version 0.3.1: Ignored; cell_velocities is the strain rate eps_dot = p_g/W.
device (str, optional) – Warp device.

Returns:

Cell kinetic energy. Shape (B,).

Return type:

wp.array

nvalchemiops.dynamics.integrators.compute_barostat_mass(target_temperature, tau_p, num_atoms, masses_out, device=None)[source]#

Compute barostat mass(es) for desired pressure fluctuation timescale.

The barostat mass determines the inertia of cell volume/shape fluctuations in NPT/NPH simulations. It is computed from the Martyna-Tobias-Klein equations [MTK1994] to give a characteristic pressure relaxation time \(\tau_p\):

\[W = (N_f + d) \cdot k_B T \cdot \tau_p^2\]

where: - N_f = 3N is the number of degrees of freedom (N atoms in 3D) - d = 3 is the dimensionality (for 3D simulations) - k_B T is the thermal energy (in reduced units, k_B = 1) - \(\tau_p\) is the pressure relaxation time

Larger W gives slower pressure equilibration (more stable but slower to reach target pressure). Smaller W gives faster equilibration but may cause oscillations.

Parameters:

target_temperature (wp.array) – Per-system target temperatures. Shape (num_systems,). Caller must pre-broadcast if a single value applies to all systems.
tau_p (wp.array) – Per-system pressure relaxation times. Shape (num_systems,). Typical values: 0.5-2.0 ps. Caller must pre-broadcast if a single value applies to all systems.
num_atoms (wp.array(dtype=wp.int32)) – Per-system atom counts. Shape (num_systems,). Caller must pre-broadcast if a single value applies to all systems.
masses_out (wp.array) – Output barostat masses W. Shape (num_systems,).
device (str, optional) – Warp device.

Returns:

Barostat mass(es) W. Shape (num_systems,).

Return type:

wp.array

Examples

Batched systems (using wp.arrays):

>>> temps = wp.array([1.0, 2.0], dtype=wp.float64, device="cuda:0")
>>> tau = wp.array([1.0, 1.0], dtype=wp.float64, device="cuda:0")
>>> n_atoms = wp.array([100, 200], dtype=wp.int32, device="cuda:0")
>>> W = wp.empty(2, dtype=wp.float64, device="cuda:0")
>>> compute_barostat_mass(temps, tau, n_atoms, W)
>>> print(W.numpy())  # [303.0, 1206.0]

Notes

The formula assumes k_B = 1 (reduced units). Scale \(\tau_p\) accordingly for real units.
For isotropic barostat, a single mass controls all cell dimensions.
For anisotropic/triclinic barostat [SSM2004], the same mass is typically used for all cell velocity components.
All input arrays must have the same length (num_systems). The caller is responsible for broadcasting scalar values to arrays before calling.

References

[MTK1994]

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

[SSM2004]

Shinoda, Shiga, Mikami, Phys. Rev. B 69, 134103 (2004)

nvalchemiops.dynamics.integrators.compute_barostat_potential_energy(target_pressures, volumes, potential_energy, device=None)[source]#

Compute barostat potential energy.

\[U = P_{\text{ext}} V\]

where \(P_{\text{ext}}\) is the target pressure and \(V\) the cell volume.

Parameters:

target_pressures (wp.array) – External/target pressures. Shape (B,).
volumes (wp.array) – Cell volumes. Shape (B,).
potential_energy (wp.array) – Output barostat potential energy. Shape (B,).
device (str, optional) – Warp device.

Returns:

Barostat potential energy. Shape (B,).

Return type:

wp.array

nvalchemiops.dynamics.integrators.npt_thermostat_half_step(eta, eta_dot, kinetic_energy, target_temperature, thermostat_masses, num_atoms_per_system, chain_length, dt, device=None)[source]#

Perform Nose-Hoover chain thermostat half-step for NPT.

Propagates the thermostat chain with the leading driving force

\[G_1 = \frac{2\,\mathrm{KE} - N_f\,k_B T}{Q_1}, \qquad G_k = \frac{Q_{k-1}\,\dot{\eta}_{k-1}^2 - k_B T}{Q_k}\]

and advances the chain positions by \(\eta_k \mathrel{+}= \tfrac{\Delta t}{2}\,\dot{\eta}_k\), where \(N_f = 3N\), \(Q_k\) are the thermostat_masses and \(\dot{\eta}_k\) the chain velocities.

Parameters:

eta (wp.array2d) – Thermostat positions. Shape (B, chain_length). MODIFIED in-place.
eta_dot (wp.array2d) – Thermostat velocities. Shape (B, chain_length). MODIFIED in-place.
kinetic_energy (wp.array) – Kinetic energy per system. Shape (B,).
target_temperature (wp.array) – Target temperatures. Shape (B,).
thermostat_masses (wp.array2d) – Thermostat masses. Shape (B, chain_length).
num_atoms_per_system (wp.array(dtype=wp.int32)) – Number of atoms per system. Shape (B,).
chain_length (int) – Number of thermostats in chain.
dt (wp.array) – Full time step dt per system. Shape (B,). The half-step and quarter-step factors are applied internally.
device (str, optional) – Warp device.

Return type:

None

nvalchemiops.dynamics.integrators.npt_barostat_half_step(cell_velocities, pressure_tensors, target_pressures, volumes, cell_masses, kinetic_energy, num_atoms_per_system, dt, device=None)[source]#

NPT barostat half-step: \(\dot{\varepsilon} += \frac{\Delta t}{2} \frac{V}{W} (P_{\text{inst}} - P_{\text{ext}})\).

Mode is dispatched by target_pressures dtype (scalar / vec3 / vec9). Barostat-NHC coupling is applied separately by the caller.

Parameters:

cell_velocities (wp.array(dtype=wp.mat33f or wp.mat33d)) – Strain-rate matrices \(\dot{\varepsilon}\). Shape (B,). MODIFIED in-place.
pressure_tensors (wp.array(dtype=vec9f or vec9d)) – Current pressure tensors from virial. Shape (B,).
target_pressures (wp.array) – External pressure(s); dtype selects mode (scalar/vec3/vec9).
volumes (wp.array(dtype=scalar)) – Cell volumes V, barostat masses W, system KE. Shape (B,).
cell_masses (wp.array(dtype=scalar)) – Cell volumes V, barostat masses W, system KE. Shape (B,).
kinetic_energy (wp.array(dtype=scalar)) – Cell volumes V, barostat masses W, system KE. Shape (B,).
num_atoms_per_system (wp.array(dtype=wp.int32)) – Atoms per system. Shape (B,).
dt (wp.array(dtype=scalar)) – Full time step; half-step factor applied internally.
device (str, optional) – Warp device. Default: inferred from cell_velocities.

Return type:

None

See also

nph_barostat_half_step: NPH (no thermostat).
npt_velocity_half_step: Particle velocity update with NHC drag.

References

See compute_barostat_mass for the MTK1994 and SSM2004 citations.

nvalchemiops.dynamics.integrators.npt_velocity_half_step(velocities, masses, forces, cell_velocities, volumes=None, eta_dots=None, num_atoms=None, dt=None, batch_idx=None, num_atoms_per_system=None, cells_inv=None, mode='isotropic', device=None)[source]#

Perform half-step velocity update for NPT ensemble.

Updates particle velocities accounting for: 1. Forces from the potential energy surface 2. Coupling to barostat (cell velocity / strain rate) 3. Coupling to thermostat (Nose-Hoover chain)

Mathematical Formulation#

Isotropic mode (default): \(\mathbf{v} \leftarrow \mathbf{v} \exp\!\left(-\Delta t \left((1 + 1/N_{\text{atoms}}) \frac{\operatorname{Tr}(\dot{\varepsilon})}{3} + \dot{\eta}_1\right)\right)\)
Anisotropic / triclinic mode: \(\mathbf{v} \leftarrow \mathbf{v} \exp\!\left(-\Delta t \left(\dot{\varepsilon} + \frac{\operatorname{Tr}(\dot{\varepsilon})}{3 N_{\text{atoms}}} I + \dot{\eta}_1 I\right)\right)\)

where eps_dot = cell_velocities = p_g/W and \(\dot{\eta}_1\) is the first thermostat-chain velocity. This primitive already applies particle thermostat drag through eta_dots[:, 0]; callers that apply particle NHC velocity scaling as a separate Trotter operator should use nph_velocity_half_step() or another no-thermostat velocity primitive instead.

param velocities:

Particle velocities. Shape (N,). MODIFIED in-place.

type velocities:

wp.array(dtype=wp.vec3f or wp.vec3d)

param masses:

Particle masses. Shape (N,).

type masses:

wp.array(dtype=scalar)

param forces:

Forces on particles. Shape (N,).

type forces:

wp.array(dtype=wp.vec3f or wp.vec3d)

param cell_velocities:

Strain-rate matrices \(\dot{\varepsilon} = p_g/W\). Shape (B,).

type cell_velocities:

wp.array(dtype=wp.mat33f or wp.mat33d)

param volumes:

Deprecated since version 0.3.1: Ignored; cell_velocities is the strain rate eps_dot = p_g/W.

type volumes:

wp.array, optional

param eta_dots:

Thermostat chain velocities. Shape (B, chain_length).

type eta_dots:

wp.array2d(dtype=scalar)

param num_atoms:

Atom count for single-system mode. Shape (1,).

type num_atoms:

wp.array(dtype=wp.int32)

param dt:

Full time step per system. Shape (B,). The half-step factor is applied internally.

type dt:

wp.array(dtype=scalar)

param batch_idx:

System index for each atom. Required for batched simulations.

type batch_idx:

wp.array(dtype=wp.int32), optional

param num_atoms_per_system:

Number of atoms per system. Required for batched simulations.

type num_atoms_per_system:

wp.array(dtype=wp.int32), optional

param cells_inv:

Deprecated since version 0.3.1: Ignored; cell_velocities is the strain rate eps_dot = p_g/W.

type cells_inv:

wp.array(dtype=wp.mat33f or wp.mat33d), optional

param mode:

Pressure control mode. One of:

"isotropic" (default): scalar coupling on \(\operatorname{Tr}(\dot{\varepsilon})/3\)
"anisotropic": diagonal coupling on \(\dot{\varepsilon}\)
"triclinic": full-tensor coupling on \(\dot{\varepsilon}\)

type mode:

str, optional

param device:

Warp device.

type device:

str, optional

Examples

Single system (isotropic):

>>> npt_velocity_half_step(
...     velocities, masses, forces, cell_velocities,
...     eta_dots=eta_dots, num_atoms=num_atoms, dt=dt,
... )

See also

npt_barostat_half_step: Cell velocity update step.
npt_position_update: Position update step.

Parameters:

velocities (array)
masses (array)
forces (array)
cell_velocities (array)
volumes (array)
eta_dots (array)
num_atoms (array)
dt (array)
batch_idx (array)
num_atoms_per_system (array)
cells_inv (array)
mode (str)
device (str)

Return type:

None

nvalchemiops.dynamics.integrators.npt_position_update(positions, velocities, cells, cell_velocities, dt, cells_inv=None, batch_idx=None, device=None)[source]#

Update particle positions for NPT integration in-place.

Applies \(\mathbf{r} \leftarrow \mathbf{r} + \Delta t\,(\mathbf{v} + \dot{\varepsilon}\,\mathbf{r})\) where cell_velocities is the strain-rate tensor \(\dot{\varepsilon} = p_g / W\).

Parameters:

positions (wp.array(dtype=wp.vec3f or wp.vec3d)) – Particle positions. Shape (N,). Modified in-place.
velocities (wp.array(dtype=wp.vec3f or wp.vec3d)) – Particle velocities. Shape (N,).
cells (wp.array(dtype=wp.mat33f or wp.mat33d)) – Cell matrices. Shape (B,).
cell_velocities (wp.array(dtype=wp.mat33f or wp.mat33d)) – Strain-rate matrices \(\dot{\varepsilon}\). Shape (B,).
dt (wp.array(dtype=scalar)) – Full time step per system. Shape (B,).
cells_inv (wp.array, optional) –

Deprecated since version 0.3.1: Ignored; cell_velocities is the strain rate eps_dot = p_g/W.
batch_idx (wp.array(dtype=wp.int32), optional) – System index for each atom. Required for batched simulations.
device (str, optional) – Warp device. Default: inferred from positions.

Return type:

None

See also

nvalchemiops.dynamics.integrators.npt.npt_position_update_out(): Non-mutating variant.
nvalchemiops.dynamics.integrators.npt.npt_cell_update(): Cell matrix update step.

nvalchemiops.dynamics.integrators.npt_cell_update(cells, cell_velocities, dt, device=None)[source]#

Update cell matrices: \(\mathbf{h}_{\text{new}} = \mathbf{h} + \Delta t\,\dot{\varepsilon}\,\mathbf{h}\).

Parameters:

cells (wp.array) – Cell matrices. MODIFIED in-place.
cell_velocities (wp.array) – Cell velocity matrices.
dt (wp.array(dtype=scalar)) – Full time step per system. Shape (B,).
device (str, optional) – Warp device.

Return type:

None

Isenthalpic-isobaric (NPH)#

NPH uses the same MTK barostat as NPT but omits the thermostat, so temperature fluctuates and the enthalpy \(H = E + PV\) is the conserved quantity. It is useful for measuring the adiabatic response of a system to an applied pressure. Control parameters are \(\Delta t\) (fs), target pressure (\(\mathrm{eV}/\mathrm{\AA}^3\)), and barostat time \(\tau_P\) (fs).

Underlying nvalchemiops kernels

nvalchemiops.dynamics.integrators.nph_barostat_half_step(cell_velocities, pressure_tensors, target_pressures, volumes, cell_masses, kinetic_energy, num_atoms_per_system, dt, device=None)[source]#

Perform barostat half-step for NPH ensemble (no thermostat coupling).

NPH (isenthalpic-isobaric) simulations maintain constant pressure and enthalpy. Unlike NPT, there is no thermostat - the temperature evolves naturally on the constant-enthalpy surface.

This function updates the cell velocity matrix \(\dot{h}\) based on the pressure difference. The pressure control mode is automatically detected from the target_pressures array dtype:

Isotropic (scalar dtype): Uniform scaling in all directions
Anisotropic/Orthorhombic (vec3 dtype): Independent x, y, z control
Triclinic (vec9 dtype): Full stress tensor control

Mathematical Formulation#

The cell velocity follows the MTK equations without thermostat:

Isotropic mode:

\[\ddot{h} = \frac{V}{W}(P - P_{\text{ext}})\]

Anisotropic mode:

\[\ddot{h}_{ii} = \frac{V}{W}(P_{ii} - P_{\text{ext},ii})\]

Triclinic mode:

\[\ddot{h}_{ij} = \frac{V}{W}(P_{ij} - P_{\text{ext},ij})\]

param cell_velocities:

Strain-rate matrices \(\dot{\varepsilon}\). Shape (B,). MODIFIED in-place.

type cell_velocities:

wp.array(dtype=wp.mat33f or wp.mat33d)

param pressure_tensors:

Current pressure tensors from virial. Shape (B,).

type pressure_tensors:

wp.array(dtype=vec9f or vec9d)

param target_pressures:

External/target pressure(s). The dtype determines the mode:

wp.float32 or wp.float64: Isotropic. Shape (B,).
wp.vec3f or wp.vec3d: Anisotropic [Pxx, Pyy, Pzz]. Shape (B,).
vec9f or vec9d: Full stress tensor. Shape (B,).

type target_pressures:

wp.array

param volumes:

Cell volumes V. Shape (B,).

type volumes:

wp.array(dtype=scalar)

param cell_masses:

Barostat masses W. Shape (B,).

type cell_masses:

wp.array(dtype=scalar)

param kinetic_energy:

System kinetic energies. Shape (B,).

type kinetic_energy:

wp.array(dtype=scalar)

param num_atoms_per_system:

Number of atoms per system. Shape (B,).

type num_atoms_per_system:

wp.array(dtype=wp.int32)

param dt:

Full time step per system. Shape (B,). The half-step factor is applied internally.

type dt:

wp.array(dtype=scalar)

param device:

Warp device.

type device:

str, optional

Examples

Isotropic NPH:

>>> target_P = wp.array([1.0], dtype=wp.float32, device="cuda:0")
>>> nph_barostat_half_step(
...     cell_velocities, pressure_tensors, target_P,
...     volumes, cell_masses, kinetic_energy,
...     num_atoms_per_system, dt=0.001
... )

Anisotropic NPH:

>>> target_P = wp.array([[1.0, 2.0, 1.5]], dtype=wp.vec3f, device="cuda:0")
>>> nph_barostat_half_step(
...     cell_velocities, pressure_tensors, target_P, ...
... )

See also

npt_barostat_half_step: Barostat with thermostat coupling.
run_nph_step: Complete NPH integration step.

Parameters:

cell_velocities (array)
pressure_tensors (array)
target_pressures (array)
volumes (array)
cell_masses (array)
kinetic_energy (array)
num_atoms_per_system (array)
dt (array)
device (str)

Return type:

None

nvalchemiops.dynamics.integrators.nph_velocity_half_step(velocities, masses, forces, cell_velocities, volumes=None, num_atoms=None, dt=None, batch_idx=None, num_atoms_per_system=None, cells_inv=None, mode='isotropic', device=None)[source]#

Perform half-step velocity update for NPH ensemble (no thermostat).

Updates particle velocities accounting for: 1. Forces from the potential energy surface 2. Coupling to barostat (cell velocity / strain rate)

Unlike NPT, there is no thermostat coupling - the temperature evolves naturally on the constant-enthalpy surface.

Mathematical Formulation#

Isotropic mode: \(\mathbf{v} \leftarrow \mathbf{v} \exp\!\left(-\Delta t (1 + 1/N_{\text{atoms}}) \frac{\operatorname{Tr}(\dot{\varepsilon})}{3}\right)\)
Anisotropic / triclinic mode: \(\mathbf{v} \leftarrow \mathbf{v} \exp\!\left(-\Delta t \left(\dot{\varepsilon} + \frac{\operatorname{Tr}(\dot{\varepsilon})}{3 N_{\text{atoms}}} I\right)\right)\)

param velocities:

Particle velocities. MODIFIED in-place.

type velocities:

wp.array(dtype=wp.vec3f or wp.vec3d)

param masses:

Particle masses.

type masses:

wp.array(dtype=scalar)

param forces:

Forces on particles.

type forces:

wp.array(dtype=wp.vec3f or wp.vec3d)

param cell_velocities:

Strain-rate matrices \(\dot{\varepsilon} = p_g/W\).

type cell_velocities:

wp.array(dtype=wp.mat33f or wp.mat33d)

param volumes:

Deprecated since version 0.3.1: Ignored; cell_velocities is the strain rate eps_dot = p_g/W.

type volumes:

wp.array, optional

param num_atoms:

Atom count for single-system mode. Shape (1,).

type num_atoms:

wp.array(dtype=wp.int32)

param dt:

Full time step per system. Shape (B,). The half-step factor is applied internally.

type dt:

wp.array(dtype=scalar)

param batch_idx:

System index for each atom.

type batch_idx:

wp.array, optional

param num_atoms_per_system:

Number of atoms per system.

type num_atoms_per_system:

wp.array, optional

param cells_inv:

Deprecated since version 0.3.1: Ignored; cell_velocities is the strain rate eps_dot = p_g/W.

type cells_inv:

wp.array(dtype=wp.mat33f or wp.mat33d), optional

param mode:

Pressure control mode:

"isotropic": scalar coupling on \(\operatorname{Tr}(\dot{\varepsilon})/3\)
"anisotropic": diagonal coupling on \(\dot{\varepsilon}\)
"triclinic": full-tensor coupling on \(\dot{\varepsilon}\)

type mode:

str, optional

param device:

Warp device.

type device:

str, optional

See also

nph_barostat_half_step: Cell velocity update.
npt_velocity_half_step: Velocity update with thermostat.

Parameters:

velocities (array)
masses (array)
forces (array)
cell_velocities (array)
volumes (array)
num_atoms (array)
dt (array)
batch_idx (array)
num_atoms_per_system (array)
cells_inv (array)
mode (str)
device (str)

Return type:

None

FIRE relaxation#

FIRE (Fast Inertial Relaxation Engine; Bitzek et al., 2006) is a geometry optimiser, not a thermostatted integrator: it drives coordinates to a local energy minimum using a damped-MD trajectory with an adaptive timestep. After each force evaluation it computes the power \(P = \sum_i \mathbf{F}_i \cdot \mathbf{v}_i\) and mixes each velocity toward the force direction,

\[\mathbf{v} \leftarrow (1 - \alpha)\,\mathbf{v} + \alpha\,\sqrt{\frac{\mathbf{v}\cdot\mathbf{v}} {\mathbf{F}\cdot\mathbf{F}}}\;\mathbf{F}.\]

While \(P > 0\) (moving downhill) the timestep grows and \(\alpha\) shrinks; when \(P \le 0\) the velocity is zeroed and the timestep is cut. Displacements are capped at maxstep (\(\mathrm{\AA}\)) and the timestep is clamped to [dt_min, dt_max] (fs). FIREVariableCell extends the same mixing to the cell degrees of freedom using NPH-style cell propagation at zero target pressure.

FIRE2 (Shuang et al., 2020) improves the restart conditions and the mixing rule; it uses a distinct set of hyperparameters (delaystep, dtgrow, dtshrink, alpha0, …) and places the whole step before the force evaluation.