Verlet Method Research Articles

Preservation of adiabatic invariants under symplectic discretization

Symplectic methods, like Verlet's method, are standard tools for long time integration of Hamiltonian systems arising, for example, in molecular dynamics. A reason for their popularity is conservation of energy over very long time up to small fluctuations that scale with the order of the method. We discuss a qualitative feature of Hamiltonian systems with separated time scales that is also preserved under symplectic discretization. Specifically, highly oscillatory degrees of freedom often lead to almost preserved quantities (adiabatic invariants). Using recent results from backward error analysis and normal form theory, we show that a symplectic method preserves those adiabatic invariants. We also discuss step size restrictions necessary to maintain adiabatic invariants in practice.

Asymptotic Error Analysis of the Adaptive Verlet Method

The Adaptive Verlet method and variants are time-reversible schemes for treating Hamiltonian systems subject to a Sundman time transformation. These methods have been observed in computer experiments to exhibit superior numerical stability when implemented in a counterintuitive “reciprocal” formulation. Here we give a theoretical explanation of this behavior by examining the leading terms of the modified equation (backward error analysis) and those of the asymptotic error expansion. With this insight we are able to improve the algorithm by simply correcting the starting stepsize.

Constant temperature simulations using the Langevin equation with velocity Verlet integration

An algorithm, which reduces to velocity Verlet in the limit of zero friction, is obtained for the generalized Langevin equation. The formulation presented is unique in that the velocities are based on a direct second order Taylor expansion of the inertial forces. The resulting finite difference equation is compared with a previous formulation of Verlet-based Langevin dynamics. The equations are implemented to study the physical properties of dense neon and liquid water at constant temperatures as a function of the friction rate γ. The results show canonical ensembles can be obtained at moderate friction rates without recurring to switching functions to control energy losses, or overdamping. The LD methodology is further applied to a protein in solution and compared to MD simulations which use a more conventional scaling function for temperature control. Analyses of the MD trajectories show the protein remains at either lower or higher temperature compared to the solvent, depending upon the treatment of the potential function at the system's boundary. On the contrary, LD simulations show a rapid and even equilibration between protein and solvent in all cases.

Variable step implementation of geometric integrators

We compare experimentally several techniques for combining geometric integrators with variable time steps. In particular, we study modifications of the Verlet method due to Leimkuhler and a technique for symplectic integration based on Poincaré transformations suggested by Hairer and Reich independently. We conclude that it is feasible to develop symplectic variable step size codes that, for Hamiltonian problems, are competitive with standard software. We also analyze the error growth of the new algorithms when integrating periodic orbits.

Symplectic integrator for molecular dynamics of a protein in water

The symplectic integrator is an algorithm for solving equations of motion, preserving the volume in phase space and ensuring a stable simulation. We carried out molecular dynamics simulations of liquid water and a protein in water using several variations of symplectic integrators. It was found that a fourth-order symplectic integrator of Calvo and Sanz-Serna generated a trajectory of much higher accuracy than the conventional Verlet and Gear methods with the same requirements for CPU time.

The Adaptive Verlet Method

We discuss the integration of autonomous Hamiltonian systems via dynamical rescaling of the vector field (reparameterization of time). Appropriate rescalings (e.g., based on normalization of the vector field or on minimum particle separation in an N-body problem) do not alter the time-reversal symmetry of the flow, and it is desirable to maintain this symmetry under discretization. For standard form mechanical systems without rescaling, this can be achieved by using the explicit leapfrog--Verlet method; we show that explicit time-reversible integration of the reparameterized equations is also possible if the parameterization depends on positions or velocities only. For general rescalings, a scalar nonlinear equation must be solved at each step, but only one force evaluation is needed. The new method also conserves the angular momentum for an N-body problem. The use of reversible schemes, together with a step control based on normalization of the vector field (arclength reparameterization), is demonstrated in several numerical experiments, including a double pendulum, the Kepler problem, and a three-body problem.

Symplectic Partitioned Runge–Kutta Methods for Constrained Hamiltonian Systems

This article deals with the numerical treatment of Hamiltonian systems with holonomic constraints. A class of partitioned Runge–Kutta methods, consisting of the couples of s-stage Lobatto IIIA and Lobatto IIIB methods, has been discovered to solve these problems efficiently. These methods are symplectic, preserve all un-derlying constraints, and are superconvergent with order $2s - 2$. For separable Hamiltonians of the form $H(q,p) = \frac{1}{2}p^T M^{ - 1} p + U(q)$ the Rattle algorithm based on the Verlet method was up to now the only known symplectic method preserving the constraints. In fact this method turns out to be equivalent to the 2-stage Lobatto IIIA–IIIB method of order 2. Numerical examples have been performed which illustrate the theoretical results.

A new molecular dynamics method combining the reference system propagator algorithm with a fast multipole method for simulating proteins and other complex systems

An efficient molecular dynamics (MD) algorithm is presented in this paper for biomolecular systems, which incorporates a novel variation on the fast multipole method (FMM) coupled to the reversible reference system propagator algorithm (r-RESPA). A top-down FMM is proposed which calculates multipoles recursively from the top of the box tree instead of from the bottom in Greengard’s original FMM, in an effort to be more efficient for noncubic or nonuniform systems. In addition, the use of noncubic box subdivisions of biomolecular systems is used and discussed. Reversible RESPA based on a Trotter factorization of the Liouville propagator in generating numerical integration schemes is coupled to the top-down FMM and applied to a MD study of proteins in vacuo, and is shown to be able to use a much larger time-step than the standard velocity Verlet method for a comparable level of accuracy. Furthermore, by using the FMM it becomes possible to perform MD simulations for very large biomolecules, since memory and CPU time requirements are now nearly of order of O(N) instead of O(N2). For a protein with 9513 atoms (the photosynthetic reaction center), the efficient MD algorithm leads to 20-fold reduction in CPU time for the Coulomb interaction and approximately 15-fold reduction in total CPU time over the standard velocity Verlet algorithm with a direct evaluation of Coulomb forces.

Qualitative study of the symplectic Störmer–Verlet integrator

Symplectic numerical integrators, such as the Störmer–Verlet method, are useful in preserving properties that are not preserved by conventional numerical integrators. This paper analyzes the Störmer–Verlet method as applied to the simple harmonic model, whose generalization is an important model for molecular dynamics simulations. Restricting our attention to the one -dimensional case, both the exact solution and the Störmer–Verlet solution to this model are expressed as functions of the number of time steps taken, and then both of these functions are interpreted geometrically. The paper shows the existence of an upper bound on the error from the Störmer–Verlet method, and then an example is worked to demonstrate the closeness of this bound.

Comparison of Parallel Verlet and Implicit Runge-Kutta Methods for Molecular Dynamics Integration

Comparison of Parallel Verlet and Implicit Runge-Kutta Methods for Molecular Dynamics Integration

Theory and algorithms for mixed Monte Carlo-stochastic dynamics simulations

The recently introduced mixed MC-SD method is a fundamentally new procedure which essentially eliminates the distinction between Monte Carlo and dynamics. Unlike other methods which utilize forces, Brownian motion or dynamical steps to generate new trial configurations in a Monte Carlo search, mixed MC-SD does stochastic dynamics on the cartesian space of a molecule and Monte Carlo on the torsion space of the molecule simultaneously. After each dynamical step, a random deformation of a rotatable torsion is performed and accepted or rejected according to the Metropolis criteria. The next dynamical step is performed from the most recent configuration and the velocities from the previous dynamical step. The smooth merging of Monte Carlo and dynamics requires the use of the stochastic velocity Verlet integration scheme. Here, the velocity Verlet stochastic dynamics method is derived, and the reasons why it can be joined with Metropolis Monte Carlo in a continuous fashion are explored.

A Non-Iterative Matrix Method for Constraint Molecular Dynamics Simulations

Abstract A non-iterative version of the matrix method which can be easily vectorised and parallelized is presented. The original matrix method, which is conceptually identical to the familiar SHAKE algorithm, is shown to be non-iterative when incorporated with the Verlet and leap-frog integration schemes with the same constraint error order as the (local) error order of the accompanying integration schemes. The method is checked by test simulations with n-butane and the liquid crystal moleculae 5CB (4-n-pentyl-4′-cyanobiphenyl) in which its effectiveness is demonstrated.

A Multiple-Time-Step Molecular Dynamics Algorithm for Macromolecules

We present a computationally efficient molecular dynamics algorithm designed to take advantage of the inherent separation of time scales in biomolecular systems. The algorithm is essentially a generalization of the previously introduced reversible reference system propagation algorithm (r-RESPA) which employs a Trotter factorization of the Liouville propagator to generate numerical integration schemes for molecular dynamics applications. The method is compared with thevelocity Verlet integration algorithm in a molecular dynamics (MD) simulation of the protein crambin in vucuo. The multiple-time-step algorithm is shown to be able to take a much larger time step for a comparable level of accuracy than that of the standard method, leading to a 4-5-fold reduction in the CPU time required to calculate the nonbonded forces.

Symplectic Numerical Integrators in Constrained Hamiltonian Systems

Recent work reported in the literature suggests that for the long-time integration of Hamiltonian dynamical systems one should use methods that preserve the symplectic (or canonical) structure of the flow. Here we investigate the symplecticness of numerical integrators for constrained dynamics, such as occur in molecular dynamics when bond lengths are made rigid in order to overcome stepsize limitations due to the highest frequencies. This leads to a constrained Hamiltonian system of smaller dimension. Previous work has shown that it is possible to have methods which are symplectic on the constraint manifold in phase space. Here it is shown that the very popular Verlet method with SHAKE-type constraints is equivalent to the same method with RATTLE-type constraints and that the latter is symplectic and time reversible. (This assumes that the iteration is carried to convergence.) We also demonstrate the global convergence of the Verlet scheme in the presence of SHAKE-type and RATTLE-type constraints. We perform numerical experiments to compare these methods with the second-order backward differentiation method, commonly recommended for ordinary differential equations with constraints.

Dangers of Multiple Time Step Methods

In this work we demonstrate two potential dangers of multiple time step techniques for numerical integration of differential equations. This idea of using different time steps for different interactions was proposed in 1978 for molecular dynamics but undoubtedly has other useful applications. The use of multiple time stepping with the popular Verlet method was proposed in a 1991 paper. However, the method advocated in this paper does not retain the symplectic (or canonical) property of the Verlet method, which is an abstract property satisfied by the flow of any Hamiltonian system, of which molecular dynamics is an example. Recent work reported in the literature suggests that this property is important for the long-time integration of Hamiltonian dynamical systems. We perform experiments on linear and nonlinear problems comparing symplectic and nonsymplectic multiple time stepping extensions of the Verlet method. We observe that in the nonsymplectic case either instability or dissipation becomes evident after a long integration. However, our experiments also indicate that it is quite possible to obtain an artificial "resonance" for the symplectic method that is much worse than that for the nonsymplectic methods.

The royal road to an energy-conserving predictor-corrector method

The performance of a fourth-order Störmer-Cowell predictor-corrector method is described. Provided that the corrector is iterated a fixed number of times, at large timesteps this method is justas stable with respect to energy drift as the well-known Verlet method, while being much more accurate. The key to arriving at a stable predictor-corrector method lies in the use of an explicitly time-reversible corrector, and in recognising that explicit time reversibility can only be achieved if the corrector is iterated until converged.

A transition-rate investigation by molecular dynamics with the Langevin/implicit-Euler scheme

We report results from molecular dynamics simulations for a bistable piecewise-harmonic potential. A new method for molecular dynamics—the Langevin/implicit-Euler scheme—is investigated here and compared to the common Verlet integration algorithm. The implicit scheme introduces new computational and physical features since it (1) does not restrict integration time step to a very small value, and (2) effectively damps vibrational modes ω≫ωc, where ωc is a chosen cutoff frequency. The main issue we explore in this study is how different choices of time steps and cutoff frequencies affect computed transition rates. The one-dimensional, double-well model offers a simple visual and computational opportunity for observing the two different damping forces introduced by the scheme—frictional and intrinsic—and for characterizing the dominating force at a given parameter combination. Another question we examine here is the choice of time step below which the Langevin/implicit-Euler scheme produces ‘‘correct’’ transition rates for a model potential whose energy distribution is ‘‘well-described’’ classically.

Generalized Verlet Algorithm for Efficient Molecular Dynamics Simulations with Long-range Interactions

Abstract For the purpose of molecular dynamics simulations of large biopolymers we have developed a new method to accelerate the calculation of long-range pair interactions (e.g. Coulomb interaction). The algorithm introduces distance classes to schedule updates of non-bonding interactions and to avoid unnecessary computations of interactions between particles which are far apart. To minimize the error caused by the updating schedule, the Verlet integration scheme has been modified. The results of the method are compared to those of other approximation schemes as well as to results obtained by numerical integration without approximation. For simulation of a protein with 12 637 atoms our approximation scheme yields a reduction of computer time by a factor of seven. The approximation suggested can be implemented on sequential as well as on parallel computers. We describe an implementation on a (Transputer-based) MIMD machine with a systolic ring architecture.

The Dufour and Soret coefficients of isotopic mixtures from equilibrium molecular dynamics calculations

The Soret and Dufour cross correlations functions (CCFs), as determined by equilibrium molecular dynamics calculations, are compared for soft sphere and Lennard-Jones (LJ) mixtures. The isotopic mixture was investigated to avoid separate determination of the partial enthalpy of the system. For the LJ system, potential truncation and particle number effects on the thermal diffusion coefficient were additionally studied. Main findings are: (i) Soret and Dufour CCFs are numerically equivalent, but the fluctuations in the latter are more pronounced, particularly for the LJ mixture. (ii) Using a symmetrized form of the tensor occurring in the heat current, differences between the two CCFs are noticeable, when the Stoermer–Verlet integration algorithm is applied. (iii) The particle number dependence of the CCF, and hence of the thermal diffusion coefficient, is very small.

4533522 Apparatus for the desulfurization of flue gases

We compare experimentally several techniques for combining geometric integrators with variable time steps. In particular, we study modifications of the Verlet method due to Leimkuhler and a technique for symplectic integration based on Poincaré transformations suggested by Hairer and Reich independently. We conclude that it is feasible to develop symplectic variable step size codes that, for Hamiltonian problems, are competitive with standard software. We also analyze the error growth of the new algorithms when integrating periodic orbits.