Integration Schemes

This section describes the numerical integration methods used in Cavity HOOMD for propagating both cavity and molecular degrees of freedom, including adaptive timestepping for handling disparate timescales and sudden protocol changes.

Overview

Cavity-coupled molecular dynamics requires careful numerical integration due to:

  1. Multiple timescales: Femtosecond cavity oscillations to picosecond molecular motion

  2. Stochastic thermostats: Both deterministic and random forces

  3. Time-varying coupling: Sudden or gradual changes in coupling strength

  4. Energy conservation: Maintaining accuracy while handling stiff equations

Cavity HOOMD employs split-operator schemes that separate deterministic propagation from stochastic thermostatting, combined with adaptive timestepping to handle protocol changes.

Velocity Verlet for Cavity Modes

Basic Algorithm

The cavity mode follows Langevin dynamics with canonical coordinates \((q_\alpha, p_\alpha)\) and unit mass. The equations of motion are:

\[\dot{q}_\alpha = p_\alpha, \quad \dot{p}_\alpha = f_\alpha(t) - \gamma_\mathrm{c} p_\alpha + \xi_\alpha(t)\]

where:

  • \(f_\alpha(t) = -\omega_\alpha^2(q_\alpha + \lambda\mu_\alpha/\omega_\mathrm{c})\): Effective force (harmonic + coupling)

  • \(\gamma_\mathrm{c} = 1/(2\tau_\mathrm{c})\): Friction coefficient

  • \(\xi_\alpha(t)\): Gaussian white noise with \(\langle\xi_\alpha^2\rangle = 6\gamma_\mathrm{c}k_BT/\Delta t\)

Split-Operator Scheme

First Half-Step (Deterministic):

Advance momentum by half timestep:

\[p_\alpha(t + \Delta t/2) = p_\alpha(t) + \frac{1}{2}f_\alpha(t)\Delta t\]

Position Update:

Update position using half-step momentum:

\[q_\alpha(t + \Delta t) = q_\alpha(t) + p_\alpha(t + \Delta t/2)\Delta t\]

Second Half-Step (Stochastic):

Modify force with friction and noise:

\[f_\alpha(t + \Delta t) \to f_\alpha(t + \Delta t) - \gamma_\mathrm{c} p_\alpha(t + \Delta t/2) + \xi_\alpha(t)\]

Complete momentum update:

\[p_\alpha(t + \Delta t) = p_\alpha(t + \Delta t/2) + \frac{1}{2}f_\alpha(t + \Delta t)\Delta t\]

Noise Variance

The random force variance satisfies the fluctuation-dissipation theorem:

\[\langle\xi_\alpha^2\rangle = \frac{6\gamma_\mathrm{c}k_BT}{\Delta t}\]

Implementation:

Sample \(\xi_\alpha \sim \mathcal{N}(0, \sqrt{6\gamma_\mathrm{c}k_BT/\Delta t})\) from a normal distribution at each timestep.

Physical Interpretation

Properties:

  • Maintains cavity at temperature \(T\)

  • Thermalizes kinetic energy at rate \(\gamma_\mathrm{c}\)

  • Represents photon loss (cavity lifetime \(\tau_\mathrm{c} = 1/(2\gamma_\mathrm{c})\))

  • Typical: \(\gamma_\mathrm{c} = 0.5\) ps⁻¹ for ultrafast energy exchange

Time-Reversible Core:

Without thermostat (\(\gamma_\mathrm{c} = 0\), \(\xi = 0\)), the scheme is time-reversible and symplectic.

Velocity Verlet for Molecules

Two-Step Propagator

Molecular coordinates \((\mathbf{R}_i, \mathbf{P}_i)\) with mass \(M_i\) evolve with Bussi-Parrinello thermostat. Each timestep \(\Delta t\) consists of:

First Half-Step:

Advance momenta by half timestep:

\[\mathbf{P}_i(t + \Delta t/2) = \mathbf{P}_i(t) + \frac{1}{2}\mathbf{F}_i(t)\Delta t\]

where \(\mathbf{F}_i(t)\) is the total force (molecular + cavity coupling).

First Velocity Rescaling:

Rescale momenta by stochastic factor \(\alpha_\upsilon\) computed from current kinetic energy:

\[K_{t+\Delta t/2} = \sum_i \frac{|\mathbf{P}_i(t + \Delta t/2)|^2}{2M_i}\]

The scaling factor is:

\[\alpha_\upsilon = \left[e^{-\Delta t/\tau_\mathrm{b}} + \frac{T}{2K_{t+\Delta t/2}}(1 - e^{-\Delta t/\tau_\mathrm{b}})(R_\Gamma + R_N^2) + 2R_N\sqrt{\frac{T}{2K_{t+\Delta t/2}}e^{-\Delta t/\tau_\mathrm{b}}(1 - e^{-\Delta t/\tau_\mathrm{b}})}\right]^{1/2}\]

where \(R_\Gamma\) and \(R_N\) are independent standard normal random variables.

Apply rescaling:

\[\mathbf{P}_i(t + \Delta t/2) \to \alpha_\upsilon \mathbf{P}_i(t + \Delta t/2)\]

Position Update:

Update positions:

\[\mathbf{R}_i(t + \Delta t) = \mathbf{R}_i(t) + \frac{\mathbf{P}_i(t + \Delta t/2)}{M_i}\Delta t\]

Force Recalculation:

Compute forces at new positions \(\mathbf{F}_i(t + \Delta t)\).

Second Half-Step:

Advance momenta:

\[\mathbf{P}_i(t + \Delta t) = \mathbf{P}_i(t + \Delta t/2) + \frac{1}{2}\mathbf{F}_i(t + \Delta t)\Delta t\]

Second Velocity Rescaling:

Apply a freshly computed rescaling factor \(\alpha'_\upsilon\):

\[\mathbf{P}_i(t + \Delta t) \to \alpha'_\upsilon \mathbf{P}_i(t + \Delta t)\]

Bussi-Parrinello Implementation

Key Parameters:

  • \(\tau_\mathrm{b} = 1.0\) ps: Thermostat time constant

  • Two rescalings per timestep ensure proper canonical sampling

  • Rescaling factors depend on instantaneous kinetic energy

Advantages:

  • Exactly samples canonical ensemble

  • Preserves short-time velocity autocorrelation

  • Minimal perturbation to microscopic dynamics

  • Ideal for supercooled liquids where dynamics are critical

Force Composition

The total force on particle \(i\) is:

\[\mathbf{F}_i = \mathbf{F}_i^\mathrm{mol} + \mathbf{F}_i^\mathrm{cavity}\]

Molecular Force:

\[\mathbf{F}_i^\mathrm{mol} = -\nabla_{\mathbf{R}_i} V\]

Includes bonds, angles, dihedrals, Lennard-Jones, Coulomb, etc.

Cavity Force:

\[\mathbf{F}_i^\mathrm{cavity} = -\sum_\alpha \left(\epsilon_{0,\alpha} q_{0,\alpha} + \frac{\epsilon_{0,\alpha}^2}{K_\alpha}D_\alpha\right) \frac{\partial d_{n\alpha}}{\partial \mathbf{R}_i}\]

where \(d_{n\alpha}\) is the component of molecular dipole moment in direction \(\alpha\).

Adaptive Timestepping

Motivation

Simulations with disparate timescales and sudden coupling changes require dynamic timestep adjustment:

  1. Femtosecond Rabi oscillations vs nanosecond structural dynamics

  2. Instantaneous step changes in coupling constant (sudden quench experiments)

  3. Numerical stability during protocol switches

Adaptive timestepping dynamically adjusts \(\Delta t\) to maintain accuracy while maximizing efficiency.

Error Estimation

Local Truncation Error:

Estimate error using explicit Euler step as reference. For particle \(i\):

Velocity-Verlet predictor:

\[\mathbf{R}_{i,\mathrm{VV}}(t + \Delta t) = \mathbf{R}_i(t) + \mathbf{v}_i(t)\Delta t + \frac{1}{2}\mathbf{a}_i(t)(\Delta t)^2\]

where \(\mathbf{v}_i = \mathbf{P}_i/M_i\) and \(\mathbf{a}_i = \mathbf{F}_i/M_i\).

Explicit Euler:

\[\mathbf{R}_{i,\mathrm{E}}(t + \Delta t) = \mathbf{R}_i(t) + \mathbf{v}_i(t)\Delta t\]

Per-Particle Error:

\[\delta_i(t) = |\mathbf{R}_{i,\mathrm{VV}} - \mathbf{R}_{i,\mathrm{E}}| = \frac{1}{2}|\mathbf{a}_i(t)|(\Delta t)^2\]

System-Wide Error:

Root-mean-square over all \(N\) particles:

\[\delta(t) = \sqrt{\frac{1}{N}\sum_{i=1}^N \delta_i(t)^2}\]

Timestep Adjustment

Update Rule:

\[\Delta t' = \Delta t \sqrt{\frac{\varepsilon^*}{\delta(t)}}\]

where \(\varepsilon^*\) is the target error tolerance.

Behavior:

  • \(\delta > \varepsilon^*\): Decrease timestep (high forces)

  • \(\delta < \varepsilon^*\): Increase timestep (weak forces)

  • \(\delta \approx \varepsilon^*\): Maintain timestep (optimal)

Practical Limits:

Impose bounds to prevent extreme changes:

\[\Delta t_\min \leq \Delta t' \leq \Delta t_\max\]

Typical: \(\Delta t_\min = 10^{-5}\) fs, \(\Delta t_\max = 2.0\) fs.

Protocol-Aware Error Control

Dynamic Error Tolerance

Sudden coupling changes cause numerical instabilities even with adaptive timestepping. Implement error tolerance ramping:

Reset and Ramp:

When coupling constant changes by \(\Delta\lambda \geq \varepsilon_\lambda\), reset error tolerance:

\[\begin{split}\varepsilon(t) = \begin{cases} \varepsilon^*, & t < t_0 \\ \varepsilon^*[1 - (1 - f_0)e^{-(t-t_0)/\tau^*}], & t \geq t_0 \end{cases}\end{split}\]

where:

  • \(t_0\): Time of coupling change

  • \(f_0\): Initial fraction (typically \(10^{-3}\))

  • \(\tau^*\): Ramping time constant (typically 50 ps)

  • \(\varepsilon^*\): Steady-state tolerance (typically 5.0)

Physical Interpretation:

  • Start with very conservative timestep (\(\Delta t \sim 10^{-4}\) fs) immediately after switch

  • Gradually increase timestep as system equilibrates

  • Reach normal timestep (\(\Delta t \sim 1.5\) fs) after ramping period

  • Balances stability during transient vs efficiency in stable regime

Implementation Details

Threshold Detection:

Monitor coupling constant each step:

\[\text{if } |\lambda(t) - \lambda(t - \Delta t)| \geq \varepsilon_\lambda \text{ then reset } \varepsilon(t)\]

Ramping Function:

Exponential approach to steady state:

\[\varepsilon(t) = \varepsilon^* - (\varepsilon^* - f_0\varepsilon^*)e^{-(t-t_0)/\tau^*}\]

Disabling Ramping:

For well-equilibrated continuation runs, set \(f_0 = 0\) (no ramping) and use fixed \(\varepsilon^*\) throughout.

Practical Parameters

Timestep Selection

Nyquist Criterion:

Resolve fastest oscillation:

\[\Delta t < \frac{\pi}{\omega_\max}\]

For \(\omega_\mathrm{c} = 2000\) cm⁻¹ \(\approx\) 377 fs period:

\[\Delta t < \frac{377\text{ fs}}{\pi} \approx 120\text{ fs}\]

In practice, use \(\Delta t \approx T_\mathrm{min}/10\) for accuracy.

Stiffness Considerations:

Bonded forces are stiffer than non-bonded. Use smaller timesteps when:

  • High bond spring constants

  • Strong coupling to cavity

  • Large dipole moment gradients

Performance Considerations

Computational Cost

Timestep vs Accuracy Tradeoff:

  • Smaller \(\Delta t\): Better accuracy, more steps

  • Larger \(\Delta t\): Faster simulation, risk instability

  • Adaptive: Best of both, overhead ~5%

Adaptive Overhead:

Error estimation adds minimal cost:

  • Compute error: \(\mathcal{O}(N)\) per step

  • Typically 3-5% overhead

  • Worthwhile for stability and efficiency

Optimal Strategy:

Use adaptive timestepping with ramping for:

  • Time-varying coupling protocols

  • Non-equilibrium experiments

  • Systems with unknown force scales

Use fixed timesteps for:

  • Well-characterized equilibrium systems

  • Maximum performance on simple systems

  • Benchmarking and validation

GPU Optimization

Memory Access Patterns:

  • Coalesced reads of positions and forces

  • Reduction for error calculation

  • Minimal branching in adaptive logic

Parallel Error Estimation:

  • Compute \(\delta_i\) for each particle in parallel

  • Reduce to \(\delta\) using shared memory

  • Update timestep on host

Typical Performance:

  • Fixed timestep: ~10⁶ steps/hour (N=500, NVIDIA A100)

  • Adaptive timestep: ~9.5×10⁵ steps/hour (5% overhead)

Validation and Testing

Energy Conservation

NVE Test:

Run without thermostats, monitor total energy:

\[\Delta E_\mathrm{rel} = \frac{|E(t) - E(0)|}{|E(0)|}\]

Acceptance Criteria:

  • \(\Delta E_\mathrm{rel} < 10^{-4}\) over 100 ps: Excellent

  • \(\Delta E_\mathrm{rel} < 10^{-3}\) over 100 ps: Acceptable

  • \(\Delta E_\mathrm{rel} > 10^{-2}\): Investigate (timestep too large, numerical issues)

Canonical Ensemble

Equipartition Test:

Check kinetic energy distribution:

\[\langle E_\mathrm{kin}\rangle = \frac{1}{2}N_f k_B T\]

where \(N_f\) is degrees of freedom.

Statistical Test:

Over long trajectory, \(E_\mathrm{kin}\) should follow expected distribution with correct mean and variance.

Time-Reversibility

Reversibility Test:

  1. Run forward for time \(t\)

  2. Reverse velocities: \(\mathbf{v} \to -\mathbf{v}\)

  3. Run backward for time \(t\)

  4. Check return to initial state

Expected:

  • With thermostats: Not reversible (stochastic)

  • Without thermostats: Should return to \(\mathbf{r}(0)\) within numerical precision

Timestep Convergence

Convergence Study:

Run same simulation with \(\Delta t\), \(\Delta t/2\), \(\Delta t/4\). Check observable \(O(\Delta t)\) converges as \(\mathcal{O}(\Delta t^2)\) (second-order scheme).

Best Practices

Initialization

Starting Configuration:

  • Generate initial positions on lattice or from equilibrated structure

  • Assign velocities from Maxwell-Boltzmann distribution at \(T\)

  • Run short equilibration with strong thermostat

Cavity Initialization:

  • Sample \(q_\alpha\), \(p_\alpha\) from thermal distribution at \(T\)

  • For harmonic oscillator: \(q \sim \mathcal{N}(0, k_BT/\omega^2)\), \(p \sim \mathcal{N}(0, k_BT)\)

Equilibration Protocol

Stage 1: Thermalization (0-10 ps):

  • Strong Bussi thermostat (\(\tau_\mathrm{b} = 0.1\) ps)

  • Fixed small timestep (\(\Delta t = 0.5\) fs)

  • Monitor temperature convergence

Stage 2: Equilibration (10-100 ps):

  • Standard thermostat (\(\tau_\mathrm{b} = 1.0\) ps)

  • Adaptive timestepping ON

  • Monitor energy, pressure, temperature

Stage 3: Production (>100 ps):

  • Begin data collection

  • Continue with production settings

  • Enable ramping if protocol involves coupling changes

Troubleshooting

Simulation Explodes:

  • Decrease \(\Delta t\) or \(\varepsilon^*\)

  • Check for particle overlaps in initial configuration

  • Enable error ramping for coupling changes

Energy Drift:

  • Reduce timestep

  • Check force cutoffs and neighbor lists

  • Verify thermostat parameters

Poor Statistics:

  • Increase number of trajectories

  • Longer equilibration before production

  • Check convergence of time averages

Next Sections

Continue to:

Related Theory: