Using nonequilibrium fluctuation theorems to understand and correct errors in equilibrium and nonequilibrium discrete Langevin dynamics simulations

Common algorithms for computationally simulating Langevin dynamics must discretize the stochastic differential equations of motion. These resulting finite time step integrators necessarily have several practical issues in common: Microscopic reversibility is violated, the sampled stationary distribution differs from the desired equilibrium distribution, and the work accumulated in nonequilibrium simulations is not directly usable in estimators based on nonequilibrium work theorems. Here, we show that even with a time-independent Hamiltonian, finite time step Langevin integrators can be thought of as a driven, nonequilibrium physical process. Once an appropriate work-like quantity is defined -- here called the shadow work -- recently developed nonequilibrium fluctuation theorems can be used to measure or correct for the errors introduced by the use of finite time steps. In particular, we demonstrate that amending estimators based on nonequilibrium work theorems to include this shadow work removes the time step dependent error from estimates of free energies. We also quantify, for the first time, the magnitude of deviations between the sampled stationary distribution and the desired equilibrium distribution for equilibrium Langevin simulations of solvated systems of varying size. While these deviations can be large, they can be eliminated altogether by Metropolization or greatly diminished by small reductions in the time step. Through this connection with driven processes, further developments in nonequilibrium fluctuation theorems can provide additional analytical tools for dealing with errors in finite time step integrators.Comment: 11 pages, 4 figure

Adaptive multi-stage integrators for optimal energy conservation in molecular simulations

We introduce a new Adaptive Integration Approach (AIA) to be used in a wide range of molecular simulations. Given a simulation problem and a step size, the method automatically chooses the optimal scheme out of an available family of numerical integrators. Although we focus on two-stage splitting integrators, the idea may be used with more general families. In each instance, the system-specific integrating scheme identified by our approach is optimal in the sense that it provides the best conservation of energy for harmonic forces. The AIA method has been implemented in the BCAM-modified GROMACS software package. Numerical tests in molecular dynamics and hybrid Monte Carlo simulations of constrained and unconstrained physical systems show that the method successfully realises the fail-safe strategy. In all experiments, and for each of the criteria employed, the AIA is at least as good as, and often significantly outperforms the standard Verlet scheme, as well as fixed parameter, optimized two-stage integrators. In particular, the sampling efficiency found in simulations using the AIA is up to 5 times better than the one achieved with other tested schemes

Numerical implementation of the exact dynamics of free rigid bodies

In this paper the exact analytical solution of the motion of a rigid body with arbitrary mass distribution is derived in the absence of forces or torques. The resulting expressions are cast into a form where the dependence of the motion on initial conditions is explicit and the equations governing the orientation of the body involve only real numbers. Based on these results, an efficient method to calculate the location and orientation of the rigid body at arbitrary times is presented. This implementation can be used to verify the accuracy of numerical integration schemes for rigid bodies, to serve as a building block for event-driven discontinuous molecular dynamics simulations of general rigid bodies, and for constructing symplectic integrators for rigid body dynamics.Comment: Shortened paper with updated references, 28 pages, 3 figure

Palindromic 3-stage splitting integrators, a roadmap

The implementation of multi-stage splitting integrators is essentially the same as the implementation of the familiar Strang/Verlet method. Therefore multi-stage formulas may be easily incorporated into software that now uses the Strang/Verlet integrator. We study in detail the two-parameter family of palindromic, three-stage splitting formulas and identify choices of parameters that may outperform the Strang/Verlet method. One of these choices leads to a method of effective order four suitable to integrate in time some partial differential equations. Other choices may be seen as perturbations of the Strang method that increase efficiency in molecular dynamics simulations and in Hybrid Monte Carlo sampling.Comment: 20 pages, 8 figures, 2 table

Optimized Verlet-like algorithms for molecular dynamics simulations

New explicit velocity- and position-Verlet-like algorithms of the second order are proposed to integrate the equations of motion in many-body systems. The algorithms are derived on the basis of an extended decomposition scheme at the presence of a free parameter. The nonzero value for this parameter is obtained by reducing the influence of truncated terms to a minimum. As a result, the new algorithms appear to be more efficient than the original Verlet versions which correspond to a particular case when the introduced parameter is equal to zero. Like the original versions, the proposed counterparts are symplectic and time reversible, but lead to an improved accuracy in the generated solutions at the same overall computational costs. The advantages of the new algorithms are demonstrated in molecular dynamics simulations of a Lennard-Jones fluid.Comment: 5 pages, 2 figures; submitted to Phys. Rev.

Efficient algorithms for rigid body integration using optimized splitting methods and exact free rotational motion

Hamiltonian splitting methods are an established technique to derive stable and accurate integration schemes in molecular dynamics, in which additional accuracy can be gained using force gradients. For rigid bodies, a tradition exists in the literature to further split up the kinetic part of the Hamiltonian, which lowers the accuracy. The goal of this note is to comment on the best combination of optimized splitting and gradient methods that avoids splitting the kinetic energy. These schemes are generally applicable, but the optimal scheme depends on the desired level of accuracy. For simulations of liquid water it is found that the velocity Verlet scheme is only optimal for crude simulations with accuracies larger than 1.5%, while surprisingly a modified Verlet scheme (HOA) is optimal up to accuracies of 0.4% and a fourth order gradient scheme (GIER4) is optimal for even higher accuracies.Comment: 2 pages, 1 figure. Added clarifying comments. Accepted for publication in the Journal of Chemical Physic