Benedict Leimkuhler

Automatic integration of Euler-Lagrange equations with constraints☆

Abstract Numerical difficulties in the integration of Euler-Lagrange and similar equations are discussed. A technique for reducing their index from three to two is introduced and it is shown that variable-order, variable-step BDF methods converge for these index two problems. The practical application of this reduction in a numerical setting is examined.

Numerical solution of differential-algebraic equations for constrained mechanical motion

SummaryThe two most popular formulations of the equations of constrained mechanical motion, thedescriptor andstate-space forms, each have severe practical limitations. In this paper, we discuss and relate some proposed reformulations of the equations which have improved numerical properties.

Symplectic splitting methods for rigid body molecular dynamics

Rigid body molecular models possess symplectic structure and time-reversal symmetry. Standard numerical integration methods destroy both properties, introducing nonphysical dynamical behavior such as numerically induced dissipative states and drift in the energy during long term simulations. This article describes the construction, implementation, and practical application of fast explicit symplectic-reversible integrators for multiple rigid body molecular simulations. These methods use a reduction to Euler equations for the free rigid body, together with a symplectic splitting technique. In every time step, the orientational dynamics of each rigid body is integrated by a sequence of planar rotations. Besides preserving the symplectic and reversible structures of the flow, this scheme accurately conserves the total angular momentum of a system of interacting rigid bodies. Excellent energy conservation can be obtained relative to traditional methods, especially in long-time simulations. The method is imple...

Algorithms for constrained molecular dynamics

In molecular dynamics simulations, the fastest components of the potential field impose severe restrictions on the stability and hence the speed of computational methods. One possibility for treating this problem is to replace the fastest components with algebraic length constraints. In this article the resulting systems of mixed differential and algebraic equations are studied. Commonly used discretization schemes for constrained Hamiltonian systems are discussed. The form of the nonlinear equations is examined in detail and used to give convergence results for the traditional nonlinear solution technique SHAKE iteration and for a modification based on successive overrelaxation (SOR). A simple adaptive algorithm for finding the optimal relaxation parameter is presented. Alternative direct methods using sparse matrix techniques are discussed. Numerical results are given for the new techniques, which have been implemented in the molecular modeling software package CHARMM and show as much as twofold improvement over SHAKE iteration.

Approximation methods for the consistent initialization of differential-algebraic equations

The algebraic constraints in a system of differential-algebraic equations (DAEs) impose a consistency requirement on the initial values that can be difficult to satisfy. In this thesis the consistency requirement is characterized by a system of equations. An approximation method is introduced for these equations, and the numerical solution of the resulting system is analyzed for certain important classes of DAEs. Finally, some numerical experiments are described.

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.

Rational Construction of Stochastic Numerical Methods for Molecular Sampling

In this article, we focus on the sampling of the configurational Gibbs-Boltzmann distribution, that is, the calculation of averages of functions of the position coordinates of a molecular

Integration Methods for Molecular Dynamics

N

The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics

-body system modelled at constant temperature. We show how a formal series expansion of the invariant measure of a Langevin dynamics numerical method can be obtained in a straightforward way using the Baker-Campbell-Hausdorff lemma. We then compare Langevin dynamics integrators in terms of their invariant distributions and demonstrate a superconvergence property (4th order accuracy where only 2nd order would be expected) of one method in the high friction limit; this method, moreover, can be reduced to a simple modification of the Euler-Maruyama method for Brownian dynamics involving a non-Markovian (coloured noise) random process. In the Brownian dynamics case, 2nd order accuracy of the invariant density is achieved. All methods considered are efficient for molecular applications (requiring one force evaluation per timestep) and of a simple form. In fully resolved (long run) molecular dynamics simulations, for our favoured method, we observe up to two orders of magnitude improvement in configurational sampling accuracy for given stepsize with no evident reduction in the size of the largest usable timestep compared to common alternative methods.

Robust and efficient configurational molecular sampling via Langevin dynamics.

Classical molecular dynamics simulation of a macromolecule requires the use of an efficient time-stepping scheme that can faithfully approximate the dynamics over many thousands of timesteps. Because these problems are highly nonlinear, accurate approximation of a particular solution trajectory on meaningful time intervals is neither obtainable nor desired, but some restrictions, such as symplecticness, can be imposed on the discretization which tend to imply good long term behavior. The presence of a variety of types and strengths of interatom potentials in standard molecular models places severe restrictions on the timestep for numerical integration used in explicit integration schemes, so much recent research has concentrated on the search for alternatives that possess (1) proper dynamical properties, and (2) a relative insensitivity to the fastest components of the dynamics. We survey several recent approaches.