npt_nose_hoover_isotropic_step

torch_sim.integrators.npt.npt_nose_hoover_isotropic_step(state, model, *, dt, kT, external_pressure)[source]

Perform a complete NPT integration step with Nose-Hoover chain thermostats.

Implements the MTK (Martyna-Tobias-Klein) NPT scheme from Tuckerman et al. (2006) [10] with Nose-Hoover chains from Martyna et al. (1996) [3].

Equations of motion (Tuckerman et al. 2006, Eqs. 1-6):

\[\begin{split}\dot{\mathbf{r}}_i &= \frac{\mathbf{p}_i}{m_i} + \frac{p_\epsilon}{W}\,\mathbf{r}_i \\ \dot{\mathbf{p}}_i &= \mathbf{F}_i - \alpha\,\frac{p_\epsilon}{W}\,\mathbf{p}_i \\ \dot{\epsilon} &= \frac{p_\epsilon}{W} \\ \dot{p}_\epsilon &= G_\epsilon = \alpha\,(2K) + \text{Tr}(\boldsymbol{\sigma}_{\text{int}})\,V - P_{\text{ext}}\,V\,d\end{split}\]

where \(\epsilon = (1/d)\ln(V/V_0)\) is the logarithmic cell coordinate, \(\alpha = 1 + d/N_f\), \(d=3\) is spatial dimension, and \(N_f = 3N - 3\) the degrees of freedom.

Symmetric propagator (Trotter factorization):

\[e^{i\mathcal{L}\Delta t} = e^{i\mathcal{L}_{\text{NHC-baro}}\frac{\Delta t}{2}} \;e^{i\mathcal{L}_{\text{NHC-part}}\frac{\Delta t}{2}} \;e^{i\mathcal{L}_2\frac{\Delta t}{2}} \;e^{i\mathcal{L}_1\Delta t} \;e^{i\mathcal{L}_2\frac{\Delta t}{2}} \;e^{i\mathcal{L}_{\text{NHC-part}}\frac{\Delta t}{2}} \;e^{i\mathcal{L}_{\text{NHC-baro}}\frac{\Delta t}{2}}\]

Position update \(e^{i\mathcal{L}_1\Delta t}\):

\[\mathbf{r}_i \leftarrow \mathbf{r}_i + \bigl(e^{v_\epsilon\Delta t} - 1\bigr)\,\mathbf{r}_i + \Delta t\,\mathbf{v}_i\,e^{v_\epsilon\Delta t/2} \,\frac{\sinh(v_\epsilon\Delta t/2)}{v_\epsilon\Delta t/2}\]

Momentum update \(e^{i\mathcal{L}_2\Delta t/2}\):

\[\mathbf{p}_i \leftarrow \mathbf{p}_i\,e^{-\alpha v_\epsilon\Delta t/2} + \frac{\Delta t}{2}\,\mathbf{F}_i\, e^{-\alpha v_\epsilon\Delta t/4} \,\frac{\sinh(\alpha v_\epsilon\Delta t/4)} {\alpha v_\epsilon\Delta t/4}\]

where \(v_\epsilon = p_\epsilon / W\) is the cell velocity.

Variable mapping (equation -> code):

Equation symbol

Code variable

\(\mathbf{r}_i\) (positions)

state.positions

\(\mathbf{p}_i\) (momenta)

state.momenta

\(m_i\) (masses)

state.masses

\(\mathbf{F}_i\) (forces)

state.forces

\(\epsilon\) (log-cell coordinate)

state.cell_position

\(p_\epsilon\) (cell momentum)

state.cell_momentum

\(W\) (cell mass)

state.cell_mass

\(\alpha\) (1 + d/Nf)

alpha (local)

\(v_\epsilon\) (cell velocity)

cell_velocities (local)

\(V_0\) (reference volume)

det(state.reference_cell)

\(G_\epsilon\) (cell force)

cell_force_val

\(P_{\text{ext}}\) (target pressure)

external_pressure

\(k_BT\) (thermal energy)

kT

\(\Delta t\) (timestep)

dt

If the center of mass motion is removed initially, it remains removed throughout the simulation, so the degrees of freedom decreases by 3.

Parameters:
Returns:

Updated state after complete integration step

Return type:

NPTNoseHooverIsotropicState

References