PaperlayerView exact revision on arXiv

Source-linked AI summary

Langevin Thermostat for Rigid Body Dynamics

Ruslan L. Davidchack, Richard Handel, M. V. Tretyakov

arXiv:0901.4725v1physics.comp-phmath.NAmath.PR

TL;DR

The paper addresses canonical sampling for rigid molecules by developing quaternion-based Langevin thermostats that directly thermalize rotational and translational degrees of freedom. It constructs second-order quasi-symplectic integrators and tests them on TIP4P water, finding faster relaxation when both translation and rotation are thermostated, while gradient-Langevin relaxation is faster but computationally less efficient.

  • Problem

    Stochastic thermostats must thermalize all types of degrees of freedom to provide direct equipartition control in rigid-body simulations.

  • Method

    The paper derives quaternion Langevin and gradient-Langevin thermostats and proposes second-order quasi-symplectic numerical integrators for them.

  • Results

    With γ = 2 −8 ps−1 and Γ = 3 −40 ps−1, Langevin relaxation is almost twice as fast when both translational and rotational degrees of freedom are thermostated.

  • Takeaways & Limitations

    The Langevin thermostat with quasi-symplectic integrators was computationally significantly more efficient than the tested gradient-Langevin approach.

Abstract

from arXiv · show

We present a new method for isothermal rigid body simulations using the quaternion representation and Langevin dynamics. It can be combined with the traditional Langevin or gradient (Brownian) dynamics for the translational degrees of freedom to correctly sample the NVT distribution in a simulation of rigid molecules. We propose simple, quasi-symplectic second-order numerical integrators and test their performance on the TIP4P model of water. We also investigate the optimal choice of thermostat parameters.

I. INTRODUCTION

The paper motivates stochastic thermostats for canonical sampling, emphasizing their direct thermalization of individual degrees of freedom. It introduces quaternion-based Langevin rigid-body equations and second-order integrators that couple rotational and translational dynamics.

  • Canonical molecular dynamics specifies temperature T instead of conserving total energy E, requiring modified equations with a thermostat.
  • Deterministic thermostats can struggle to equilibrate systems whose fast and slow parts exchange energy slowly.Nosé-Hoover-type methods directly control global quantities and rely on efficient internal energy transfer.
  • Stochastic thermostats independently thermalize degrees of freedom, providing direct equipartition control without relying on efficient internal energy transfer.
  • The paper introduces quaternion Langevin equations for rigid-body dynamics and effective second-order quasi-symplectic integrators.The equations can be coupled with Langevin or Brownian translational dynamics.
  • The proposed methods are tested on TIP4P water while examining thermostat parameters and discretization error.

II. EQUATIONS OF MOTION

The rigid-body model uses center-of-mass coordinates and unit quaternions for translation and rotation. Its Hamiltonian dynamics preserve quaternion normalization, which the stochastic thermostats are designed to retain while sampling the Gibbsian measure.

  • Rigid molecules are represented by center-of-mass coordinates r and rotational quaternion coordinates q with |q_j| = 1.
  • The Hamiltonian includes center-of-mass and angular momenta together with potential interaction energy and rotational kinetic energy.
  • The rotational dynamics use constant 4-by-4 matrices S_l to construct quaternion-based motion.
  • Hamiltonian evolution preserves unit quaternion length when the initial conditions satisfy |q_j(0)| = 1.
  • The stochastic thermostats are formulated as ergodic SDEs with Gibbsian invariant measure while preserving the quaternion constraint.The inverse temperature is β = 1/(k_BT).

A. Langevin-type equations

The Langevin thermostat adds stochastic translational and rotational coupling to the heat bath while preserving quaternion normalization and targeting the Gibbsian measure. The paper specializes the coefficients for simpler additive-noise computation and analyzes their parameter-dependent relaxation.

  • The Langevin-type equations are written in Itô form for translational and rotational rigid-body variables.
  • The friction coefficients γ and Γ control translational and rotational coupling to the heat bath, while g, G, b, and B define noise coupling.
  • Independent noise is assigned to degrees of freedom, with diagonal noise matrices b and B assumed for the subsequent analysis.
  • The stochastic dynamics preserve the quaternion constraint |Q_j(t)| = 1 for every molecule and time.
  • The paper restricts computational treatment to constant b_ii and B_ii coefficients because additive-noise systems yield simpler numerical methods.
  • The Langevin thermostat is defined with γ and Γ as free heat-bath coupling strengths, and rotational angular velocities are related to quaternion derivatives through S_iQ.

B. A mixture of gradient system and Langevin-type equation

The gradient-Langevin thermostat combines a gradient subsystem for translation with Langevin dynamics for rotation while preserving the relevant invariant measure. Its parameter ν sets relative translational evolution speed, but its numerical scheme is more costly than the Langevin schemes.

  • The gradient-Langevin approach uses a gradient system for translational dynamics and a Langevin equation for rotational dynamics.
  • The center-of-mass momenta can be omitted from differential-equation simulation because they are independent Gaussian variables under the Gibbsian measure.
  • Its invariant measure is the Gibbsian measure integrated over translational momenta, and quaternion normalization is preserved.
  • The parameter ν controls the gradient subsystem’s speed relative to rotational Langevin dynamics rather than defining an intrinsic dynamical time.
  • For fixed rotational coordinates, the gradient translational subsystem can be viewed as an overdamped limit of the translational Langevin subsystem.

III. NUMERICAL INTEGRATORS

The paper develops quasi-symplectic numerical methods for Langevin-type rigid-body dynamics and characterizes their structural and weak-convergence properties. The methods preserve or approximate key phase-space structures while supporting weak-order analysis.

  • Method construction: The paper constructs second-order numerical methods for both the Langevin and gradient-Langevin thermostats.The methods target stochastic rigid-body systems while retaining the structural properties needed for simulations near Hamiltonian dynamics.
  • Phase-volume behavior: exp(−n(3γ + Γ) · t) gives the Jacobian determinant of the Langevin flow, quantifying phase-volume contraction under translational and rotational friction.The phase volume is preserved when γ = 0 and Γ = 0; otherwise, with γΓ > 0, the flow is contractive.
  • Quasi-symplectic structure: Quasi-symplectic integrators are designed to remain close to symplectic methods when stochastic Langevin dynamics approaches Hamiltonian dynamics.They degenerate to a symplectic method when the Langevin system becomes Hamiltonian, while reflecting phase-volume contractivity through their Jacobian behavior.
  • Weak convergence: Weak numerical methods are used because they suffice for evaluating averages and are simpler than mean-square methods for direct trajectory simulation.Weak convergence is defined through expectations of sufficiently smooth test functions over a fixed time interval.

A. Numerical schemes for the Langevin thermostat

Two quasi-symplectic schemes are constructed for the Langevin thermostat by combining stochastic Hamiltonian integration with damping and rotational updates. Both preserve the quaternion constraint, achieve weak order two, and exactly reproduce the continuous system’s Jacobian.

  • Scheme construction: Two Langevin integrators are derived using splitting techniques together with a numerical method for the deterministic rigid-body Hamiltonian system.The first splits the dynamics into a Hamiltonian system with additive noise and a deterministic linear damping system; the second is constructed for strong heat-bath coupling.
  • Rotational update: The rotational update uses the mapping Ψ_l(t; q, π), which advances quaternion and conjugate-momentum variables through trigonometric rotations.The map is defined by Q = cos(χ_l t)q + sin(χ_l t)S_lq and Π = cos(χ_l t)π + sin(χ_l t)S_lπ.
  • Jacobian preservation: The one-step Jacobian of each proposed Langevin integrator is exactly equal to the Jacobian of the original stochastic system.This equality is established for both the first and second schemes, rather than merely being an asymptotic approximation.
  • Structural properties: Both Langevin schemes are quasi-symplectic, preserve the quaternion constraint |Q_j| = 1, and have weak order two.These properties are stated for the first scheme in Proposition 1 and for the second scheme in Proposition 2.
  • Parameter regime: The second scheme is usually preferable when γ and/or Γ are large, although matrix-factorization work adds implementation cost.The additional cost is described as negligible for molecular systems where force calculations dominate computational effort.

B. Numerical scheme for the gradient-Langevin system

A second-order weak scheme is developed for the gradient-Langevin system by combining a second-order additive-noise Runge–Kutta treatment with the rotational Langevin approximation. It preserves the quaternion constraint but requires more force evaluations.

  • Scheme construction: The gradient-Langevin scheme combines a second-order Runge–Kutta method for the additive-noise gradient part with the rotational approximation used in the Langevin scheme.The resulting method is a second-order weak scheme for the coupled gradient-translational and rotational dynamics.
  • Structural properties: The gradient-Langevin numerical scheme preserves the structural property |Q_j| = 1 and has weak order two.These properties are stated in Proposition 3 for scheme (34), (25).
  • Computational cost: Two force calculations per step make the gradient-Langevin scheme approximately twice as expensive as the Langevin schemes.The extra force evaluation is the stated computational cost of the gradient-Langevin approach.

C. Computational errors

The paper estimates ergodic averages by balancing finite-time, numerical-integration, and Monte Carlo errors. The relevant controls are simulation time, time step and method choice, and the number of independent trajectories.

  • Ergodic estimation: Ergodic averages are estimated from long-time trajectories under the assumption that the stochastic systems are ergodic with the Gibbsian invariant measure.The paper further focuses on systems satisfying exponential ergodicity.
  • Error decomposition: The total estimation error has three components: finite-time approximation error, numerical integration error Ch^p, and Monte Carlo error.The decomposition separates errors from approaching the invariant average, discretizing the dynamics, and sampling independent realizations.
  • Error control: Large T controls the finite-time error, h and the numerical method control integration error, and L controls statistical error.T is chosen to approach the invariant distribution and reduce bias, while the number of trajectories regulates Monte Carlo variability.
  • Averaging horizon: The trajectory length required for ergodic averaging must also reduce estimator variance, so it is much larger than the time needed only to approach stationarity.The paper distinguishes the long averaging horizon ˜T from the finite-time horizon T used in distributional approximation.

IV. NUMERICAL INVESTIGATION

The numerical study evaluates thermostat sampling and relaxation for TIP4P water, focusing on parameter choices and relaxation after a temperature switch. Relaxation rates differ across translational temperature, rotational temperature, and potential energy, with potential energy governing equilibration.

  • Study design: The study tests Langevin and gradient-Langevin thermostats, their parameter dependence, and discretization errors using the TIP4P rigid-water model.The investigation varies γ, Γ, and ν and examines integration step size h for three proposed numerical schemes.
  • Evaluation criteria: The thermostat must correctly sample the Gibbs distribution at the target temperature and rapidly relax the system to equilibrium.Sampling accuracy is assessed from long-run averages of kinetic and potential energies, pressure, and related properties.
  • Relaxation experiment: Relaxation is measured after 2000 TIP4P water molecules are switched from T0 = 220 K to T1 = 270 K and evolved to the new equilibrium.The experiment uses an instantaneous thermostat-temperature increase and tracks several measured quantities.
  • Relaxation analysis: The characteristic relaxation time is estimated by fitting an exponential response, and system equilibration is determined by the largest τ A among quantities of interest.The measured quantities include translational and rotational kinetic temperatures and potential energy per molecule.
  • Numerical results: τ Ttr = 0.2 ps, τ Trot = 1.9 ps, and τ U = 3.6 ps for γ = 4.0 ps−1 and Γ = 0, making potential energy the slowest reported quantity.With translational thermostatting only, translational temperature relaxes quickly while rotational temperature and potential energy lag behind.
  • Numerical results: Relaxation is slower at both small and large γ, with the fastest relaxation around γ = 4.0 ps−1.Small γ limits heat flux, whereas large γ can hinder configurational relaxation through the random force.

A. Langevin Thermostats

The study maps how Langevin parameters affect relaxation in TIP4P water and evaluates discretization behavior of two second-order integrators. Relaxation is fastest with translational and rotational thermostating together, while excessive coupling and large time steps introduce reliability or stability problems.

  • Thermostat-parameter dependence: Relaxation speeds for translational and rotational temperatures increase with γ and Γ, respectively, whereas potential-energy relaxation depends nonuniformly on both parameters.The parameter sweep uses five independent runs per logarithmic-grid point and a 0.2 fs time step to reduce discretization effects.
  • Thermostat-parameter dependence: The fastest relaxation occurs when both translational and rotational degrees of freedom are thermostated, with γ = 2–8 ps−1 and Γ = 3–40 ps−1.At γ = 4.0 ps−1 and Γ = 10.0 ps−1, the estimated relaxation times are τ Ttr = 0.28 ps, τ Trot = 0.26 ps, and τ U = 2.0 ps.
  • Integrator performance: Langevin A has larger rotational-temperature discretization error but smaller potential-energy error than Langevin B.Both methods retain linear dependence of measured quantities on h^2 up to about h = 7 fs.
  • Integrator performance: Both Langevin A and B become unstable at about h = 10 fs.The discretization-error coefficients C_A are reported in Table I for the measured quantities under the Fig. 6 conditions.
  • Parameter-dependent numerical behavior: The parameter study also examines ⟨Ttr⟩0, ⟨Trot⟩0, and ⟨U⟩0 together with their discretization-error coefficients for Langevin A and B.The translational-temperature extrapolation converges to the thermostat temperature, while the straight-line fit breaks down at large Γ for Langevin A.

B. Gradient-Langevin thermostat

The gradient-Langevin thermostat efficiently transfers heat between translational gradient dynamics and rotational motion, but its numerical scheme has restrictive stability limits.

  • Numerical stability: At h = 0.2 fs, the gradient-Langevin scheme is unstable when hν exceeds about 200 fs2, limiting ν to approximately 1000 fs in these simulations.The bound applies to the particular TIP4P water system studied.
  • Relaxation behavior: Rotational-temperature relaxation decreases as Γ increases, while potential-energy relaxation shows little dependence on Γ except at very large values.The rotational and potential-energy relaxation times were measured after switching the thermostat temperature from 220 K to 270 K.
  • Relaxation behavior: For ν = 100 fs and Γ = 0, τ Trot = 0.14 ps while τ U = 2.7 ps, indicating efficient heat transfer between translational and rotational dynamics.The rotational relaxation is much faster than in the corresponding Langevin case, whereas potential-energy relaxation is similar.
  • Relaxation behavior: Potential-energy relaxation time decreases with increasing ν without reaching a minimum within the explored range.This dependence differs markedly from that observed for the Langevin system.
  • Discretization error: Measured ⟨Trot⟩h and ⟨U⟩h depend linearly on h2, with estimated slopes CTrot = −0.38 K/fs2 and CU = −0.029 kcal/mol/fs2.For ν used in the discretization-error study, the integrator becomes unstable for h > 1 fs, compared with stability up to about 10 fs for the Langevin A and B integrators.

V. SUMMARY

The paper introduces quaternion-based stochastic thermostats for rigid-body models and evaluates them on TIP4P water. Gradient-Langevin relaxation is faster, but Langevin thermostating with quasi-symplectic integrators is computationally more efficient because of better stability.

  • V. SUMMARY: The proposed Langevin and gradient-Langevin systems thermostat translational and rotational degrees of freedom while preserving unit quaternion length.The thermostats are designed for quaternion-based rigid-body dynamics and couple both types of motion to a heat bath.
  • V. SUMMARY: Gradient-Langevin relaxation times are smaller than Langevin relaxation times.This faster relaxation is accompanied by weaker numerical stability for the proposed gradient-Langevin scheme.
  • V. SUMMARY: Langevin thermostating with the quasi-symplectic integrators was computationally significantly more efficient than gradient-Langevin thermostating with its numerical scheme.The comparison reflects the experimental study and the differing stability properties of the numerical methods.
Loading 0901.4725v1…