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 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.
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.
(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.
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.
Article metrics loading...
Full text loading...