Linear Analysis of Stochastic Verlet-Type Integrators for Langevin Equations

TL;DR

This work formulates a general, objective framework to evaluate stochastic Verlet-type integrators for Langevin dynamics by deriving closed-form expressions for three configurational benchmarks: diffusion on a flat surface, drift on a tilted surface, and Boltzmann sampling in a harmonic potential. By expressing integrators in a configurational form with coefficients $c_i$ and noise correlations, the authors compute normalized measures $\Gamma_{\rm diff}$, $\Gamma_{\rm drift}$, and $\Gamma_{\rm dist}$ and identify the unique $GJ_{13-20}$ set as capable of achieving unity for all three benchmarks within stability. The paper systematically analyzes twelve historical integrators, revealing that many fail to reproduce all three statistics at finite time steps, while BAOAB$_{12}$ and BBK$_{84}$ offer reliable transport properties but can incur sampling-temperature deviations; in contrast, the $GJ_{13-20}$ family achieves exact configurational statistics, providing a principled basis for selecting thermostats in complex, nonlinear systems. The findings offer a practical, problem-independent guide for high-quality thermodynamic simulations and emphasize the importance of linear-benchmark validation before applying integrators to nonlinear problems.

Abstract

We provide an analytical framework for analyzing the quality of stochastic Verlet-type integrators for simulating the Langevin equation. Focusing only on basic objective measures, we consider the ability of an integrator to correctly simulate two characteristic configurational quantities of transport, a) diffusion on a flat surface and b) drift on a tilted planar surface, as well as c) statistical sampling of a harmonic potential. For any stochastic Verlet-type integrator expressed in its configurational form, we develop closed form expressions to directly assess these three most basic quantities as a function of the applied time step. The applicability of the analysis is exemplified through twelve representative integrators developed over the past five decades, and algorithm performance is conveniently visualized through the three characteristic measures for each integrator. The GJ set of integrators stands out as the only option for correctly simulating diffusion, drift, and Boltzmann distribution in linear systems, and we therefore suggest that this general method is the one best suited for high quality thermodynamic simulations of nonlinear and complex systems, including for relatively high time steps compared to simulations with other integrators.

Figures (12)

• Figure 1: Normalized discrete-time diffusion [$\Gamma_{\rm diff}$, Eqs. (\ref{['eq:DE_discrete']}) and (\ref{['eq:SS_G_diff']})], drift [$\Gamma_{\rm drift}$, Eqs. (\ref{['eq:Drift_discrete']}) and (\ref{['eq:SS_G_drift']})], and Boltzmann configurational temperature [$\Gamma_{\rm dist}$, Eqs. (\ref{['eq:Dist_discrete']}) and (\ref{['eq:SS_G_dist']})] for the SS$_{78}$ integrator of Ref. SS, as a function (a) of $\gamma\Delta{t}$ for $\Gamma_{\rm diff}$ (solid) and $\Gamma_{\rm drift}$ (dashed), and (b) of $\Omega_0\Delta{t}$ for $\Gamma_{\rm dist}$, the latter for select values of normalized damping $\gamma/\Omega_0$ as indicated on the figure. Vertical arrows in (b) indicate the stability limit of the time step as given by Eq. (\ref{['eq:Stability_eq_R-']}). Discrete-time quantities are normalized to the correct continuous-time quantities, Eq. (\ref{['eq:Basic_truths']}).
• Figure 2: Normalized discrete-time diffusion [$\Gamma_{\rm diff}$, Eqs. (\ref{['eq:DE_discrete']}) and (\ref{['eq:EB_III_G_diff']})], drift [$\Gamma_{\rm drift}$, Eqs. (\ref{['eq:Drift_discrete']}) and (\ref{['eq:EB_III_G_drift']})], and configurational temperature [$\Gamma_{\rm dist}$, Eqs. (\ref{['eq:Dist_discrete']}) and (\ref{['eq:EB_III_G_dist']})] for the EB$_{80}$ integrator of Ref. Ermak1980, as a function (a) of $\gamma\Delta{t}$ for $\Gamma_{\rm diff}$ (solid) and $\Gamma_{\rm drift}$ (dashed), and (b) of $\Omega_0\Delta{t}$ for $\Gamma_{\rm dist}$, the latter for select values of normalized damping $\gamma/\Omega_0$ as indicated on the figure. Vertical arrows in (b) indicate the stability limit of the time step as given by Eqs. (\ref{['eq:Stability_eq_C']}) and (\ref{['eq:Stability_eq_R-']}). Discrete-time quantities are normalized to the correct continuous-time quantities, Eq. (\ref{['eq:Basic_truths']}).
• Figure 3: Normalized discrete-time diffusion [$\Gamma_{\rm diff}$, Eqs. (\ref{['eq:DE_discrete']}) and (\ref{['eq:REB_III_G_diff']})], drift [$\Gamma_{\rm drift}$, Eqs. (\ref{['eq:Drift_discrete']}) and (\ref{['eq:REB_III_G_drift']})], and Boltzmann configurational temperature [$\Gamma_{\rm dist}$, Eqs. (\ref{['eq:Dist_discrete']}) and (\ref{['eq:REB_III_G_dist']})] for the MPA$_{80-82}$ integrator of Ref. Allen80Allen82, as a function (a) of $\gamma\Delta{t}$ for $\Gamma_{\rm diff}$ (solid) and $\Gamma_{\rm drift}$ (dashed), and (b) of $\Omega_0\Delta{t}$ for $\Gamma_{\rm dist}$, the latter for select values of normalized damping $\gamma/\Omega_0$ as indicated on the figure. Vertical arrow in (b) indicates the stability limit, $\Omega_0\Delta{t}<2$, of the time step as given by Eq. (\ref{['eq:Stability_eq_R-']}). Discrete-time quantities are normalized to the correct continuous-time quantities, Eq. (\ref{['eq:Basic_truths']}).
• Figure 4: Normalized discrete-time diffusion [$\Gamma_{\rm diff}$, Eqs. (\ref{['eq:DE_discrete']}) and (\ref{['eq:vg82_G_diff']})], drift [$\Gamma_{\rm drift}$, Eqs. (\ref{['eq:Drift_discrete']}) and (\ref{['eq:vg82_G_drift']})], and Boltzmann configurational temperature [$\Gamma_{\rm dist}$, Eqs. (\ref{['eq:Dist_discrete']}) and (\ref{['eq:vg82_G_dist']})] for the vGB$_{82}$ integrator of Ref. vGB82, as a function (a) of $\gamma\Delta{t}$ for $\Gamma_{\rm diff}$ (solid) and $\Gamma_{\rm drift}$ (dashed), and (b) of $\Omega_0\Delta{t}$ for $\Gamma_{\rm dist}$, the latter for select values of normalized damping $\gamma/\Omega_0$ as indicated on the figure. Vertical arrows in (b) indicate the stability limit of the time step as given by Eq. (\ref{['eq:Stability_eq_R-']}). Discrete-time quantities are normalized to the correct continuous-time quantities, Eq. (\ref{['eq:Basic_truths']}).
• Figure 5: Normalized discrete-time diffusion [$\Gamma_{\rm diff}$, Eqs. (\ref{['eq:DE_discrete']}) and (\ref{['eq:BBK_G_diff']})], drift [$\Gamma_{\rm drift}$, Eqs. (\ref{['eq:Drift_discrete']}) and (\ref{['eq:BBK_G_drift']})], and Boltzmann configurational temperature [$\Gamma_{\rm dist}$, Eqs. (\ref{['eq:Dist_discrete']}) and (\ref{['eq:BBK_G_dist']})] for the BBK$_{84}$ integrator of Ref. BBK, as a function (a) of $\gamma\Delta{t}$ for $\Gamma_{\rm diff}$ (solid) and $\Gamma_{\rm drift}$ (dashed), and (b) of $\Omega_0\Delta{t}$ for $\Gamma_{\rm dist}$, the latter for any value of normalized damping $\gamma/\Omega_0>0$. The stability limit of the time step is given by Eq. (\ref{['eq:Stability_eq_R-']}) to be $\Omega_0\Delta{t}<2$. Discrete-time quantities are normalized to the correct continuous-time quantities, Eq. (\ref{['eq:Basic_truths']}).