Theory · 3 methods

Thermostats

What CAPS computes, as it computes it: the equation, its symbols with CAPS's defaults, when to use the method, the source file, the tests that check it, and any departure from the cited method.

Stochastic velocity rescaling (CSVR)

CSVR rescales all velocities by one factor each step. The factor is drawn so that the kinetic energy follows the stochastic equation below, which samples the canonical ensemble exactly while disturbing the dynamics as little as a Berendsen thermostat does. CAPS draws the new kinetic energy in one step from its exact distribution, as in the paper's appendix.

$$\mathrm{d}K = (\bar{K} - K)\,\frac{\mathrm{d}t}{\tau} + 2\sqrt{\frac{K\bar{K}}{N_f}}\,\frac{\mathrm{d}W}{\sqrt{\tau}}$$

Bussi et al. 2007, eq. 7

Symbol	Meaning	In CAPS
K	instantaneous kinetic energy	—
K̄	target kinetic energy, N_f k_B T / 2	T
τ	thermostat relaxation time	τ_T · default 100 fs
N_f	degrees of freedom: 3N − 3, held atoms excluded	computed
dW	Wiener noise: one Gaussian and a χ² sum of N_f − 1 per step	mt19937-64 stream from the run's seed

When to use it

The default for NVT and NPT in CAPS. The energy it exchanges with the bath is tracked, so the conserved quantity reports integration error. Prefer Langevin (BAOAB) for far-from-equilibrium starts that need strong damping.

Source: core/src/dynamics.cpp
Tested by: Dynamics.ThermostatsHoldTheTemperature · Dynamics.BussiConservedQuantity
Departure from the reference: none

References

Bussi, G., Donadio, D., Parrinello, M., "Canonical sampling through velocity rescaling", J. Chem. Phys. 126, 014101 (2007). doi:10.1063/1.2408420

Langevin dynamics (BAOAB)

Each step splits into half kicks (B), half drifts (A) and an exact Ornstein–Uhlenbeck update of the velocities (O) in the order B A O A B. The splitting gives accurate configurational averages at larger friction than other Langevin schemes.

$$\text{O:}\quad v \leftarrow c_1 v + c_2\sqrt{\frac{k_\mathrm{B}T}{m}}\,\xi, \qquad c_1 = e^{-\gamma\Delta t}, \qquad c_2 = \sqrt{1 - c_1^2}$$

Leimkuhler & Matthews 2013

Symbol	Meaning	In CAPS
γ	friction, 1/τ_T	τ_T · default 100 fs
ξ	standard Gaussian per component	mt19937-64 stream from the run's seed
T	target temperature	T (or a linear ramp)

When to use it

Robust thermalisation of stiff or strained starts. The friction slows diffusion and relaxation, so measure dynamics (MSD, relaxation times) with CSVR instead.

Source: core/src/dynamics.cpp
Tested by: Dynamics.ThermostatsHoldTheTemperature
Departure from the reference: none

References

Leimkuhler, B., Matthews, C., "Rational construction of stochastic numerical methods for molecular sampling", Appl. Math. Res. Express 2013, 34–56 (2013). doi:10.1093/amrx/abs010

Nosé–Hoover chains and MTK pressure coupling

A chain of three thermostat variables couples to the particles' kinetic energy, each link thermostatting the one before; the equations are deterministic and time-reversible, and the extended system's energy is conserved. For constant pressure the logarithm of the volume becomes a dynamical variable with its own mass and chain (Martyna–Tobias–Klein), driven by the difference between the instantaneous and target pressure. CAPS integrates both exactly as LAMMPS's fix nvt and fix npt iso do: chain half-steps (Trotter splitting) around velocity Verlet, the cell scaled in two half-steps around the drift.

$$\begin{gathered}\dot{p}_i = F_i - \left(\dot{\eta}_1 + \left(1 + \frac{3}{N_f}\right)\dot{\epsilon}\right)p_i \\ Q_1\ddot{\eta}_1 = \sum_i \frac{p_i^2}{m_i} - N_f k_\mathrm{B}T - Q_1\dot{\eta}_1\dot{\eta}_2 \\ W\ddot{\epsilon} = 3V(P - P_0) + \frac{3}{N_f}\sum_i \frac{p_i^2}{m_i} - W\dot{\epsilon}\,\dot{\eta}_{\mathrm{p}1}\end{gathered}$$

Martyna, Klein & Tuckerman 1992; Martyna, Tobias & Klein 1994

Symbol	Meaning	In CAPS
Q₁, Qₖ	chain masses	N_f k_B T τ_T², k_B T τ_T² (τ_T the damping time)
W	volume mass	(N + 1) k_B T τ_P²
chain length	links per chain	3, as LAMMPS

When to use it

Production NVT and NPT runs that must match LAMMPS (or GROMACS Nose-Hoover/MTTK) exactly; for equilibration from a poor start the Bussi thermostat and stochastic cell rescaling relax faster.

Source: core/src/dynamics.cpp (nhc_temp, nhc_press, omega_step, remap)
Tested by: Dynamics.NoseHooverChainsAndMtkBarostat · scripts/check_lammps.sh (fix nvt, fix npt iso: positions within 1.5e-5 Å after 400 steps)
Departure from the reference: isotropic pressure coupling only; not with bond constraints or r-RESPA

References

Martyna, G. J., Klein, M. L., Tuckerman, M., "Nosé–Hoover chains: the canonical ensemble via continuous dynamics", J. Chem. Phys. 97, 2635–2643 (1992). doi:10.1063/1.463940
Martyna, G. J., Tobias, D. J., Klein, M. L., "Constant pressure molecular dynamics algorithms", J. Chem. Phys. 101, 4177–4189 (1994). doi:10.1063/1.467468