Rational Construction of Stochastic Numerical Methods for Molecular Sampling

TL;DR

The paper tackles sampling the configurational Gibbs-Boltzmann distribution for molecular systems at constant temperature by deriving a formal invariant-measure expansion for Langevin-splitting integrators via the Baker-Campbell-Hausdorff lemma. It demonstrates a superconvergence phenomenon, yielding effectively 4th-order accuracy in configurational sampling in the high-friction limit for the BAOAB-type method, and presents a simple Brownian-dynamics limit with a coloured-noise modification. The authors provide explicit expressions for the invariant-density corrections, compare several integrators (BAOAB, ABOBA, SPV, BBK), and show that BAOAB delivers substantial accuracy gains with a single force evaluation per step, along with stability at larger timesteps. Numerical experiments on oscillator models and small molecular clusters confirm up to two orders of magnitude improvement in configurational sampling accuracy and illustrate the method’s robustness across large friction values. Overall, the work offers a practical, high-accuracy framework for molecular sampling that leverages BCH-based invariant-measure expansions and targeted integrator design to achieve efficient, accurate configurational statistics.

Abstract

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 $N$-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.

Figures (5)

• Figure 1: The graphs show the comparison of four different Langevin dynamics methods when applied using different stepsizes. The configurational distribution errors are plotted against the stepsize in a log-log scale. Here $k_BT=1$. The simulation time was fixed for all runs at $t= 5 \times10^7$, and five runs were averaged to further reduce sampling errors. At left, $\gamma=1$, at right $\gamma=50$. The graphs are entirely in keeping with the theory presented in the article.
• Figure 2: Computed distributions of the 1D model problem are compared for $\gamma=20$. The three curves for each of the four methods show the results for three different stepsizes: $\delta t=0.1$ (circles), $\delta t=0.2$ (crosses), $\delta t=0.3$ (dashed), compared to the dark, solid curve representing the exact distribution. The graph shows that the generalized GLA methods are superior to SPV and BBK in the moderate $\gamma$ regime.
• Figure 3: The diagrams illustrate the distributions of interatomic distances, $G(r)$, for Morse (left) and Lennard-Jones (right) clusters. They also show the choice of bins used in calculating the numerical distributions.
• Figure 4: The radial distribution errors are plotted in log-log scale against stepsize, demonstrating the first order decay of the error in the case of Euler-Maruyama and the second-order behavior of the BAOAB limit method (in modified timestep $h$). Left: Morse potential; Right: Lennard-Jones. For Morse we used a temperature of $k_BT=0.1$, a fixed time interval of $t=4\times 10^6$, with stepsizes ranging from $0.0075$ to $0.0225$. For Lennard-Jones the temperature was $k_B T=0.2$, $t=2.5\times 10^5$ and stepsizes ranged from $.001$ to $0.0022$. In order to drive the variance of the results down, a large number of runs were necessary: for the Morse simulation the error is computed using the average histogram computed from 200 independent runs, while the Lennard-Jones simulation used 1000 independent runs. Around two orders of magnitude of improvement are observed in the accurate regime, but, perhaps even more important, the BAOAB limit method is usable at substantially larger stepsizes than Euler-Maruyama.
• Figure 5: The diagrams illustrate the performance of the different algorithms (labelled) by showing the error in the computed radial distribution functions as a function of both $\gamma$ and the number of timesteps (samples) taken, in the case of the 7 atom Lennard-Jones model. The graphs address the potential concern that the larger values of $\gamma$ needed to give the superconvergence property may reduce the rate of convergence to equilibrium (it does not, in the case of BAOAB). The same stepsize of $\delta t=0.044$ was used for all these simulations.