The challenge of stochastic Størmer–Verlet thermostats generating correct statistics

Abstract

In light of the recently developed complete GJ set of single random variable stochastic, discrete-time Størmer–Verlet algorithms for statistically accurate simulations of Langevin equations [N. Grønbech-Jensen, Mol. Phys. 118, e1662506 (2020)], we investigate two outstanding questions: (1) Are there any algorithmic or statistical benefits from including multiple random variables per time step and (2) are there objective reasons for using one or more methods from the available set of statistically correct algorithms? To address the first question, we assume a general form for the discrete-time equations with two random variables and then follow the systematic, brute-force GJ methodology by enforcing correct thermodynamics in linear systems. It is concluded that correct configurational Boltzmann sampling of a particle in a harmonic potential implies correct configurational free-particle diffusion and that these requirements only can be accomplished if the two random variables per time step are identical. We consequently submit that the GJ set represents all possible stochastic Størmer–Verlet methods that can reproduce time step-independent statistics of linear systems. The second question is thus addressed within the GJ set. Based on numerical simulations of complex molecular systems, as well as on analytic considerations, we analyze apparent friction-induced differences in the stability of the methods. We attribute these differences to an inherent, friction-dependent discrete-time scaling, which depends on the specific method. We suggest that the method with the simplest interpretation of temporal scaling, the GJ-I/GJF-2GJ method, be preferred for statistical applications.

Bringing discrete-time Langevin splitting methods into agreement with thermodynamics

“…It is worth clarifying that the use of the second, independent noise contribution β n+ 1 2 within the time step is exclusively for the half-step velocity, and it is therefore not in conflict with the result of Ref. [34], stating that there can only be one stochastic variable per time step for the configurational coordinate r n to be statistically correct. We notice that this new half step velocity coincides with the previously identified statistically correct half-step velocity for the GJ-III method, where c 3 = 1 [29] (see also Eq.…”

“…The configurational sampling is, however, not limited by the additivity of the O operation, since Ref. [34] has shown that no stochastic Verlet-type integrator can produce correct Boltzmann statistics with more than one independent stochastic number contributing to the configurational equation per time step. If O were not at the pivotal (i.e., central) position, but instead used as two symmetrically positioned 1 2 ∆t operations, then the configurational variable r n+1 would be subject to two separate stochastic contributions, σ n + and σ n+1 − .…”

“…Numerous direct simulations of both simple nonlinear and complex molecular dynamics systems have since confirmed these attractive features (see, e.g., Refs. [28][29][30][31][32][33][34]) and two variants of this method [26,28] are included in the public distribution of the molecular dynamics simulation suite LAMMPS [35,36].…”

“…Most recently, Ref. [34] has explicitly shown that a stochastic Verlet-type method cannot produce correct time-step-independent configurational statistics if more than one stochastic number is applied to each degree of freedom per time step. This crucial result underlines that the set of GJ methods is, indeed, a complete and exclusive set of methods that have the above mentioned desired statistical properties.…”

“…We note that the GJ methods have previously been extensively tested on nonlinear and classical MD systems; see, e.g., Refs. [28,30,31,33,34].…”