^{1}and Charles Matthews

^{1,a)}

### Abstract

A wide variety of numerical methods are evaluated and compared for solving the stochastic differential equations encountered in molecular dynamics. The methods are based on the application of deterministic impulses, drifts, and Brownian motions in some combination. The Baker-Campbell-Hausdorff expansion is used to study sampling accuracy following recent work by the authors, which allows determination of the stepsize-dependent bias in configurational averaging. For harmonic oscillators, configurational averaging is exact for certain schemes, which may result in improved performance in the modelling of biomolecules where bond stretches play a prominent role. For general systems, an optimal method can be identified that has very low bias compared to alternatives. In simulations of the alanine dipeptide reported here (both solvated and unsolvated), higher accuracy is obtained without loss of computational efficiency, while allowing large timestep, and with no impairment of the conformational exploration rate (the effective diffusion rate observed in simulation). The optimal scheme is a uniformly better performing algorithm for molecular sampling, with overall efficiency improvements of 25% or more in practical timestep size achievable in vacuum, and with reductions in the error of configurational averages of a factor of ten or more attainable in solvated simulations at large timestep.

We thank David Hardy (University of Illinois) for his support with the modification of the NAMD package. We also appreciate the support of the Lorentz Center (Leiden, NL) and the programme on “Modelling the Dynamics of Complex Molecular Systems” which supported the authors and provided valuable interactions during the preparation of the article. This work has made use of the resources provided by the Edinburgh Compute and Data Facility (http://www.ecdf.ed.ac.uk/). The ECDF is partially supported by the eDIKT initiative (http://www.edikt.org.uk). We further acknowledge the support of the Engineering and Physical Sciences Research Council which has funded this work as part of the Numerical Algorithms and Intelligent Software Centre under Grant No. EP/G036136/1.

I. INTRODUCTION

II. BACKGROUND

III. PERFORMANCE OF LANGEVIN ALGORITHMS APPLIED TO THE HARMONIC OSCILLATOR

IV. ERROR ANALYSIS FOR GENERAL SYSTEMS

V. NUMERICAL RESULTS

A. One-dimensional double well

B. Alanine dipeptide

1. Unsolvated

2. Solvated

VI. CONCLUSION

### Key Topics

- Friction
- 18.0
- Brownian dynamics
- 16.0
- Molecular dynamics
- 14.0
- Phase space methods
- 12.0
- Peptides
- 10.0

## Figures

The computed potential energy distribution is shown for the method of Brünger, Brooks, and Karplus applied to a single alanine dipeptide protein at 300 K in a vacuum using the CHARMM22 forcefield; the energy distribution becomes distorted as the step size increases.

The computed potential energy distribution is shown for the method of Brünger, Brooks, and Karplus applied to a single alanine dipeptide protein at 300 K in a vacuum using the CHARMM22 forcefield; the energy distribution becomes distorted as the step size increases.

The configurational density function is shown for a planar two-basin potential, computed using two different Langevin dynamics methods at various (stable) timesteps increasing from left to right. The color indicates the value of the computed probability density: from high (red) to low (blue) over the unit square. The methods have the same computational cost (in terms of force evaluations), but give different results at large timesteps. Details on this computation are given in the supplementary material. 13

The configurational density function is shown for a planar two-basin potential, computed using two different Langevin dynamics methods at various (stable) timesteps increasing from left to right. The color indicates the value of the computed probability density: from high (red) to low (blue) over the unit square. The methods have the same computational cost (in terms of force evaluations), but give different results at large timesteps. Details on this computation are given in the supplementary material. 13

The numerically computed average error in ⟨q 2⟩, using the 1D harmonic oscillator with a given Langevin dynamics method. Computation was fixed at 107 total force evaluations, with M = γ = β = K = 1. The results are averaged over 2000 independent repeat runs, with error bars included to give the standard deviation in these results. The exactness property of the ABOBA and BAOAB schemes result in the error decreasing as step size increases due to sampling error.

The numerically computed average error in ⟨q 2⟩, using the 1D harmonic oscillator with a given Langevin dynamics method. Computation was fixed at 107 total force evaluations, with M = γ = β = K = 1. The results are averaged over 2000 independent repeat runs, with error bars included to give the standard deviation in these results. The exactness property of the ABOBA and BAOAB schemes result in the error decreasing as step size increases due to sampling error.

The BAOAB scheme is shown to significantly reduce configurational discretization errors, which distort averages in Langevin dynamics, for the one-dimensional double-well system with potential energy function U(q) = (q 2 − 1)2 + q. Errors are computed from the average of 500 independent trajectories, with 109 total iterations per trajectory. The step sizes tested began at δt = 0.2 and were increased by 5% incrementally until all schemes became unstable. The relative error in the temperature (computed by averaging the momenta) is given in (a), with the scheme of Bussi/Parrinello giving a high order relationship with the step size. This is contrasted in (c) where the temperature is calculated using the instantaneous system position; here it is instead the BAOAB scheme that gives a high order result. The error in configurational distribution (calculated as a root mean squared deviance in the histogram bins) is shown in (b), with a sample computed distribution given in (d) for δt = 0.25. The exact distribution is shown as a dashed line, with the inset magnifying the density of the deepest well.

The BAOAB scheme is shown to significantly reduce configurational discretization errors, which distort averages in Langevin dynamics, for the one-dimensional double-well system with potential energy function U(q) = (q 2 − 1)2 + q. Errors are computed from the average of 500 independent trajectories, with 109 total iterations per trajectory. The step sizes tested began at δt = 0.2 and were increased by 5% incrementally until all schemes became unstable. The relative error in the temperature (computed by averaging the momenta) is given in (a), with the scheme of Bussi/Parrinello giving a high order relationship with the step size. This is contrasted in (c) where the temperature is calculated using the instantaneous system position; here it is instead the BAOAB scheme that gives a high order result. The error in configurational distribution (calculated as a root mean squared deviance in the histogram bins) is shown in (b), with a sample computed distribution given in (d) for δt = 0.25. The exact distribution is shown as a dashed line, with the inset magnifying the density of the deepest well.

Results from 2.5 ns simulations of alanine dipeptide in vacuum at the given step size (horizontal) and friction (vertical). Pixels are colored according to relative errors for each simulation, with white pixels indicating instability.

Results from 2.5 ns simulations of alanine dipeptide in vacuum at the given step size (horizontal) and friction (vertical). Pixels are colored according to relative errors for each simulation, with white pixels indicating instability.

The number of barrier recrossings (top) and the diffusion coefficients (bottom) are shown for each simulation, the latter in m2/s and computed by integrating the velocity autocorrelation function over an interval of 1 ps. As expected, the computed coefficients do not vary significantly between methods, though changing the friction above 1/ps has a substantial effect.

The number of barrier recrossings (top) and the diffusion coefficients (bottom) are shown for each simulation, the latter in m2/s and computed by integrating the velocity autocorrelation function over an interval of 1 ps. As expected, the computed coefficients do not vary significantly between methods, though changing the friction above 1/ps has a substantial effect.

(a) Configurational temperature error. (b) Average total potential energy error. Numerical results from 5 ns runs of alanine dipeptide solvated in a 10A sphere of TIP3P water are shown, using the given algorithms. Errors are computed against a baseline solution averaged from ten 5 ns simulations using the SPV scheme at δt = 0.5 fs. Results from a single run are shown for each scheme except in the case of the BAOAB method. The BAOAB scheme shows an order of magnitude improvement in the error in computed average total potential energy; because of the small absolute errors, we exhibit the means and standard deviations from 10 runs for each step size used. The lack of an observed trend line for BAOAB suggests that the discretization error is being dominated by the sampling error.

(a) Configurational temperature error. (b) Average total potential energy error. Numerical results from 5 ns runs of alanine dipeptide solvated in a 10A sphere of TIP3P water are shown, using the given algorithms. Errors are computed against a baseline solution averaged from ten 5 ns simulations using the SPV scheme at δt = 0.5 fs. Results from a single run are shown for each scheme except in the case of the BAOAB method. The BAOAB scheme shows an order of magnitude improvement in the error in computed average total potential energy; because of the small absolute errors, we exhibit the means and standard deviations from 10 runs for each step size used. The lack of an observed trend line for BAOAB suggests that the discretization error is being dominated by the sampling error.

## Tables

The expected long-time computed average of q 2 using each Langevin dynamics scheme, for the 1D harmonic oscillator U(q) = Kq 2/2. For brevity, some results are shown as leading order series in δt.

The expected long-time computed average of q 2 using each Langevin dynamics scheme, for the 1D harmonic oscillator U(q) = Kq 2/2. For brevity, some results are shown as leading order series in δt.

Numerical results for ten 2.5 ns simulations of unsolvated alanine dipeptide, with friction set to γ = 1/ps. The mean and standard deviations of all simulations are given. The baseline comparison run was completed by averaging ten 2.5 ns simulations using δt = 0.25 fs with the SPV scheme. Sampling error will play a large role in the determination of these averages, but it is clear that the BAOAB scheme outperforms the others by a significant margin. The number of observed recrossings is obtained by counting the number of times the central dihedral angles in the alanine dipeptide model hop between their two configurations.

Numerical results for ten 2.5 ns simulations of unsolvated alanine dipeptide, with friction set to γ = 1/ps. The mean and standard deviations of all simulations are given. The baseline comparison run was completed by averaging ten 2.5 ns simulations using δt = 0.25 fs with the SPV scheme. Sampling error will play a large role in the determination of these averages, but it is clear that the BAOAB scheme outperforms the others by a significant margin. The number of observed recrossings is obtained by counting the number of times the central dihedral angles in the alanine dipeptide model hop between their two configurations.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content