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

Over the past several decades, discrete-time Langevin and Molecular Dynamics (MD) simulations have provided a wealth of information about the properties of nonlinear and complex systems.^1–5 While the simulations are intended to represent the continuous-time equations of motion, the inevitable temporal discretization alters not only the accuracy of the simulated trajectories but in some cases also the fundamental aspects of the system itself. Thus, an integral part of any simulation task is to explore and optimize the balance between the two conflicting objectives; namely, simulation efficiency by increasing the discrete time step and simulation accuracy by decreasing the discrete time step. For computational statistical mechanics, one of the key equations of motion to investigate is the Langevin equation,^6,7 which is the topic we are concerned with in this paper,

where m is the mass of an object with spatial (configurational) coordinate r and velocity $v=ṙ$⁠. The object is subjected to a force f and linear friction, which is represented by the non-negative constant α. The fluctuation–dissipation relationship specifies that the thermal fluctuations β can be represented through a temporally uncorrelated Gaussian variable⁸ 

where δ(t) is Dirac’s delta function, k_B is Boltzmann’s constant, and T is the thermodynamic temperature.

Many methods (thermostats) for controlling the temperature of a simulated system have been developed and most of them fall into two major categories: Deterministic (e.g., velocity rescaling⁹ or Nosé–Hoover^10–14) and stochastic (Langevin) thermostats (see the large body of work represented in Refs. 15–30). It has also been proposed to use hybrids of those two approaches in which the individual degree of freedom is largely deterministic, but with the input from a global measure of kinetic temperature, which is stochastically regulated.^31,32 Common for these methods, and all other methods until recently, is that they include at least second order time step errors in the configurational sampling statistics. As a result, reliable simulations of thermodynamic properties necessitate rather small time steps compared to a given stability limit. The difficulty with obtaining robust, large time step statistics is rooted in the discrete-time inconsistencies between a coordinate and its momentum (velocity), the latter being defined by an imperfect finite difference of the former (see Refs. 33–36 for discussions). As a result, validation and calibration of the imperfect kinetic temperature are inconsistent with the accompanying configurational temperature, leading to misinterpretations of the simulated temperature for configurational sampling. This is a particular difficulty for methods such as velocity rescaling⁹ or Nosé–Hoover,^10–12 where the observed discrete-time kinetic temperature is adjusted to a desired value. In this context, it is interesting to recall that discrete-time dynamics is governed solely by the spatial coordinate as described by the second order difference equation (the position-Verlet method) usually attributed to Størmer and Verlet^37,38 but identified much earlier by Newton (see, e.g., the recent review in Ref. 39). The later introduction of approximate finite difference velocities has traditionally been done in two ways: one leading to the on-site velocity-explicit^40,41 method and the other leading to the half-step leap-frog^42,43 method. The advantage of explicit inclusion of a velocity variable is obviously that one has direct access to the kinetic measures, but these approximate velocities induce the above-mentioned ambiguities of how to interpret those measures in discrete time (see Ref. 36 for a discussion).

Focusing directly on configurational sampling in convex potentials, Leimkuhler and Matthews⁴⁴ developed the BAOAB method, which to our knowledge is the first to correctly sample Boltzmann statistics in a harmonic potential. Soon after, the stochastic GJF algorithm⁴⁵ demonstrated correct, time step-independent configurational statistics for linear systems; that is, correct drift and configurational diffusion in addition to the Boltzmann distribution for a harmonic potential. Due to the above-mentioned discrete-time inconsistencies, it was shown that the accompanying on-site velocity coordinate exhibits statistical deviations that increase quadratically with the time step, thereby depressing kinetic measures in simulations with correct configurational sampling. By introducing a time scale revision to the BAOAB method, a very similar algorithm followed⁴⁶ with seemingly similar properties to the GJF method. A thermostat that can give correct thermodynamic response for both configurational and kinetic core quantities was demonstrated by the GJF-2GJ method,³⁵ which adopted the configurational sampling of the GJF method with the addition of a tailored true half-step velocity that responds with correct, time step-independent Maxwell–Boltzmann distribution for as long as the configurational sampling is correct. It was further demonstrated that the diffusion measured by the velocity autocorrelation function (ACF) also is correct if the appropriate Riemann sum is chosen for approximating the Green–Kubo integral.⁴⁹ The GJF-2GJ method is now available within the LAMMPS suite.^47,48 This method was recently generalized by the complete GJ set of stochastic, Verlet-type methods.³⁶ These methods all display correct, time step-independent configurational statistics, and they have accompanying true half-step velocities that also respond with correct kinetic statistics. Reference 36 further illuminates that it is not possible to find an on-site velocity that correctly samples the Maxwell–Boltzmann distribution if the configurational coordinate is statistically correct.

Common for the GJ methods is that they include a single stochastic number for each degree of freedom in each time step. This is not unreasonable, since the Langevin equation (1) has a single random sequence over a time step. However, direct integration of the Langevin equation over one time step reveals that this single sequence of randomness translates into two different, correlated noise terms in the discrete-time equations for the coordinate and its velocity¹⁶ (see also  Appendixes A and  B in this paper). Thus, it is obvious to methodically explore what benefits, if any, the inclusion of two random variables per time step can add to the quality of statistical simulations and if it is possible to go beyond the quality that is exhibited by the GJ methods. Once it is understood which methods may be available for statistically robust sampling in discrete time, we further wish to determine if one can rationally discriminate among the available methods in order to select one or more methods with particularly desirable features. The objective in constructing methods that are statistically robust against the applied time step revolves around linear analysis and properties. We notice that for nonlinear systems, which are the typical systems for which one would apply the methods, additional discrete-time complexity enters into the dynamics.⁵⁰ As these new methods offer the possibility for applying large time steps, certain important quantities may unexpectedly suffer as a result. Thus, it is important to validate and investigate any developed method against systems that represent the true complexity and nonlinearity that the method needs to handle properly. In support of such validation, recent publications have investigated the statistical properties of the GJ methods applied to molecular dynamics with direct comparisons to earlier methods. For example, a review of the statistical accuracy of the methods in Refs. 18, 29, and 35 was given in Ref. 35, comparisons between different GJ methods were shown in Ref. 36, and a comprehensive review of the configurational properties of the methods in Refs. 18, 44, and 45 was recently given in Ref. 51. Additionally, by implementing the GJF method into several established simulation suites, AMBER,⁵² LAMMPS,^47,48 GROMACS,²³ and ESPResSo,⁵³ Refs. 54 and 55 have demonstrated the beneficial features of the configurational sampling of the GJF method compared with other, widely used stochastic thermostats.

This presentation starts by reviewing the direct integration approaches¹⁶ to the Langevin equation and how that leads to a discretization with two random variables per time step. Following Ref. 36, we then analyze a general stochastic form of the discrete-time Størmer–Verlet equations including two random variables with mutual correlation. We conclude that correct configurational statistics for linear systems can only be accomplished if the two random variables are fully correlated, implying that only a single random variable per degree of freedom should be applied to simulations. This result points to the recently published complete set of GJ methods as representing all discrete-time options in the stochastic Størmer–Verlet form that can reliably provide robust statistics for relatively large time steps. We finally analyze the GJ methods for their mutual relationship and practical functionality. We observe a direct scaling relationship between the methods and verify this scaling through comprehensive simulations of a complex molecular system in equilibrium.

The seemingly most direct approach to constructing discrete-time Langevin algorithms for simulations is to analytically integrate the Langevin equation (1) over one time step Δt, from t_n to t_n+1 = t_n + Δt. This can be done in at least two ways, which we have illustrated separately in  Appendixes A and  B, respectively. Both of these approaches consist of exact integrals of the Langevin equation (1) except for symmetric trapezoidal approximations to the integral over the conservative force f in (1). With rⁿ ≈ r(t_n) and $v$ⁿ ≈ $v$(t_n) and fⁿ representing the conservative force at time t_n, the immediate results of the two approaches can be expressed as follows:

• The result of approach A is given in Eqs. (A16a) and (A16c),

where $c̃2$ and $c̃3$ are given in Eqs. (A10) and (A11), respectively. The noise terms $d̃rβrn+1$ and $βAvn+1$ are given in Eqs. (A12) and (A4) (see  Appendix A for details on parameters and derivations).

• B.

  The result of approach B is given in Eqs. (B10a) and (B10b),

where the coefficient $d̃Bv$ is given in Eq. (B5) and where the noise terms $d̃rβrn+1$ and $d̃BvβBvn+1$ are given in Eqs. (A12) and (B4). Approach B coincides with that of Ref. 16 (see  Appendix B for details on parameters and derivations).

The two sets of Eqs. (A16) and (B10) [or, equivalently, Eqs. (3) and (4)] are identical except for the appearance of the noise term in the second equation of each set. However, by reviewing the definitions of the noise terms in Eqs. (A4), (A12), and (B4), we observe that

Thus, the two sets of equations, describing methods A and B, are identical. Curiously, implementation of this direct integration method (we will refer to the direct integration method in the form given in  Appendix B) shows that it does not reproduce time step-independent statistics if the conservative force represents a potential with non-zero curvature. However, the method and its derivation, which is correct in its fluctuation–dissipation relationship, allude to the rationale for considering the inclusion of two different noise terms for each time step, as these two noise terms represent different convolutions, one in velocity and the other in position, of the integrated Langevin noise over the time step Δt. Thus, the following analysis generalizes the functional forms of these equations in order to investigate the possibilities for creating algorithms that can take advantage of two different noise terms per time step. This analysis will illuminate why the direct integration approach fails and which possibilities for useful algorithms exist.

Inspired by the two direct integrations, A and B, outlined above, we write the general expressions in the form offered by  Appendix B,¹⁶ 

where c₂ is the pivotal one-time step velocity attenuation factor, and the other functional parameters are to be determined. The two stochastic terms are given by

where the Gaussian random variables $σrn$ and $σvn$ are independent across different time steps and for the same step possess a joint probability distribution with covariance so that

where δ_n,ℓ is Kronecker’s delta function. Note that ζ is the correlation between the two noise terms of the same time step, with |ζ| ≤ 1. The choice |ζ| = 1 corresponds to the complete set of GJ methods,³⁶ which has a single, shared noise term per time step. As promoted in Refs. 35 and 36, the discrete-time velocity is an inherently ambiguous quantity, which for any given configurational coordinate rⁿ can be defined with different features in mind. Thus, in order to design and evaluate the statistics of a method, it is convenient to first express this method without the direct representation of velocity, thereby focusing exclusively on configurational statistics before the possibility of correct kinetic statistics is evaluated. We therefore rewrite Eq. (6) in the form

where we recognize the one time step attenuation factor c₂ from the GJ set of methods.³⁶ The initial relationships between parameters are

$ForT=0̲$⁠, we write the velocity $v$ⁿ in its general, central difference form

where γ₁ + γ₂ + γ₃ = 0. Inserting Eq. (9) into Eq. (11) yields

Comparing this expression with Eq. (6b) implies that γ₃$c3′$ = 0. It is further reasonable to require that the method can correctly mimic drift for T = 0 and fⁿ = f = const. For fⁿ = f = const and T = 0, we write the constant drift velocity $v$_d = f/α

 Appendix C investigates the special case γ₃ = 0, which allows $c3′$ ≠ 0, and concludes that $c3′$ ≠ 0 cannot lead to a meaningful method that supports a correct and time step-independent Boltzmann distribution. Moreover, in the case of an explicit velocity, it can be shown that lim_α→0$c3′$ = 0 is a necessary condition for the conservative method to be symplectic. Thus, we will proceed with the condition $c3′$ = 0 and

in conjunction with the GJ set of methods.³⁶ 

Using these observations, we can hereby write the equation for the trajectory rⁿ in Eq. (9) as

where c₅ = c₂c₆. Given the one time step attenuation factor, c₂, there remain three parameters (c₄, c₆, and ζ) in order to determine the noise terms such that the method correctly represents basic thermodynamic measures. Notice that the GJ set of methods³⁶ is the same as Eq. (16) for c₄ = c₁c₃, c₆ = c₃, and ζ = 1 (i.e., $βrn=βvn$⁠), where ζ expresses the correlation given in Eq. (8b). The functional parameter c₂ can be chosen among decaying functions for αΔt/m > 0 and functions that satisfy $c2→1−αΔtm$ for αΔt/m → 0,³⁶ and this is the parameter that distinguishes the GJ methods from each other [see Eq. (53) below for key examples]. The resulting method can only be meaningful for |c₂| < 1, since c₂ is the velocity attenuation factor over the time Δt. The two core statistical measures that we initially wish to reproduce are diffusion, as defined by the configurational Einstein definition, for a flat potential (fⁿ = 0), and the Boltzmann distribution for a harmonic potential (fⁿ = −κrⁿ with κ > 0).

For fⁿ = 0, we follow Refs. 35 and 36 and define the velocity $vn+12$⁠,

such that

Combining the velocity Eq. (17) with Eq. (16) gives

which, when inserted into Eq. (18), produces the diffusive displacement

The configurational diffusion is thus

Thus, basic diffusion requires that the noise parameters satisfy the condition

Notice that this condition is satisfied for the direct integration method,¹⁶ described in  Appendixes A and  B, and summarized in Eqs. (3)–(5).

For $fn=−κrn=−Ω02mrn$⁠, where $Ω0=κ/m$ is the natural frequency of the harmonic oscillator with spring constant κ > 0, we adopt the methodology of Ref. 36 to find the autocorrelation ⟨rⁿrⁿ⟩ from the linearized version of Eq. (16). In order to simplify the visual impression of the expressions, we define

and

such that the linearized Eq. (16) becomes

Multiplying Eq. (27) with rⁿ⁻¹, rⁿ, and rⁿ⁺¹, respectively, and then making the statistical averages of the resulting equations, we get

which can also be written

The autocorrelation of immediate interest is then directly given by

With

Equation (30) can be written

Thus, we must require both of the following conditions to be satisfied:

We notice that adding twice Eq. (34b) to Eq. (34a) yields precisely Eq. (24). Thus, the requirement Eq. (24) for correct free particle diffusion is redundant with the set Eq. (34) that enforces the correct Boltzmann distribution for the harmonic potential. It is noteworthy to retroactively recognize that this redundancy between the statistics of a confined and a free particle is rooted in c₃ from Eq. (15), which in turn follows from the basic expectation that constant deterministic drift is correct in discrete time.

Since the three conditions, given in Eqs. (24) and (34), determining the three parameters of the stochastic terms in Eq. (16), are redundant such that, e.g., Eq. (24) can be discarded, we investigate only the two remaining oscillator conditions in Eq. (34), which can be rewritten as

We here recognize the single random variable GJ result (c₆ = c₃, c₄ = c₁c₃, ζ = 1) as a solution to the two equations. Defining the two variables, z₆ and z₄,

we write Eq. (35) in the simple form

The task of identifying the solutions (c₄, c₆, ζ) now becomes one of finding the intersections between the two curves in Eq. (37). As illustrated in Fig. 1, it is obvious that the curve described by Eq. (37a) is a unit circle and that the curve described by Eq. (37b) is outside the unit circle for both small and large positive values of z₆. Inserting Eq. (37b) into Eq. (37a) gives

While the left-hand side of this equation is non-negative, the right-hand side is non-positive. Thus, the intersection can only happen if

i.e., any discrete-time method of the form Eq. (6), or equivalently Eq. (16), that is required to satisfy the basic thermodynamic requirements of Einstein diffusion and Boltzmann distribution in harmonic potentials must have noise correlation |ζ| = 1. Notice that the two cases ζ = ±1 simply imply the sign of z₄z₆, which is a distinction with no significance for the resulting algorithm. We therefore choose to set ζ = 1, which implies z₄z₆ > 0. It follows that there does not exist any method that can be expressed in the form Eq. (6) for which two different random variables can accomplish the required basic thermodynamic measures for the coordinate rⁿ. The unique solution ζ = 1 is the one for which c₆ = c₃ and c₄ = c₁c₃. This unique possibility coincides with the complete GJ set of statistically sound methods,³⁶ which has one random variable per time step.

Having learned from the previous sections that no stochastic Størmer–Verlet-based thermostat can yield time step-independent statistics with two different random variables per time step, we conclude that the complete GJ set of stochastic thermostats,³⁶ using one stochastic variable per time step, is the only possible avenue for statistically accurate simulations of the Langevin equation for Størmer–Verlet type algorithms. This set is given by a free functional parameter, the one time step velocity attenuation factor c₂, which is a function of the reduced variable αΔt/m. Using the result from Sec. III C, the method can be written from Eq. (16) with ζ = 1 (⁠$βrn=βvn=βn$⁠), c₄ = c₁c₃, and c₅ = c₂c₆ = c₂c₃,

The unique on-site $v$ⁿ and half-step $un+12$ velocities for this method were found in Ref. 36 to be

respectively, where on-site and half-step velocities are defined³⁶ such that

in accordance with the expectations from statistical mechanics.⁷ While it is not possible to find an on-site velocity that can provide correct, time step-independent sampling of the kinetics, the unique half-step velocity was found to give

for the noisy harmonic oscillator, regardless of Δt.^35,36 The combination of Eqs. (42) and (43) yields the specific velocity-explicit GJ method³⁶ 

This is the resulting form that determines the functional parameters in the general expression (6) that we set out to investigate in this work. The resulting noise parameters that should be compared to the parameters $d̃r$⁠, $d̃Bv$⁠, and ζ_B for the direct integration method reviewed in  Appendix B (Fig. 8) are d_r = c₃, $dv=c1c3$⁠, and ζ = 1. It is noticeable that these noise parameters are completely different from the ones arising from direct integration of the noise (see  Appendix B). Furthermore, while direct integration points to an exponential velocity attenuation factor $c̃2$ in line with the GJ-II method, it also points to a parameter relationship [d₂ = d₃, d₆ = c₂, and d₇ = 1 in Eq. (6)], which is aligned with the GJ-I method where c₁ = c₃ [see Eq. (53) below for relevant c₂ functions or Ref. 36 for details on the methods]. Thus, the result of direct integration, leading to two different noise variables of each time step, is curiously inconsistent with basic statistical properties that are found in the GJ set.

Since the GJ methods are the only possibilities that will provide correct, time step-independent statistics, we seek to determine if there are any of those methods that are especially advantageous. As indicated in Ref. 36, the stability and harmonic period of these methods are the same if time and damping are scaled with the quantity $c3/c1$ [see Eqs. (56) and (70) in Ref. 36]. More generally, even beyond the linear limit, Eq. (48) can be written

where

and

Consequently, the GJ set of methods can be mapped into a simpler form with revised parameters $Δt̃$ and $α̃$⁠. The scaled form (49) has the same linear stability limit $Ω0Δt̃<2$ for any value of $α̃$⁠. We denote a method with a friction-independent stability limit as temporally unscaled. Within the GJ set of methods, GJ-I is the only such method,³⁶ and the BAOAB method⁴⁴ is also temporally unscaled with the same stability limit as GJ-I.

Three methods, GJ-I (GJF and GJ-I are the same for configurational sampling),^35,36,45 GJ-II (VRORV and GJ-II are the same for configurational sampling),^36,46 and GJ-III,³⁶ are exemplified in Fig. 2, where Fig. 2(a) shows the defining velocity attenuation factor c₂ for each method,

Figure 2(b) shows the corresponding scaling of time (or frequency³⁶), which in turn scales the stability limit. Notice that because αΔt/m appears explicitly in the parameter c₂, this mapping does not imply that the GJ set can be represented as a single method by scaling time and friction by the same factor. Nevertheless, consistent with the GJ method being temporally unscaled, the form of Eq. (49) coincides with the GJ-I³⁶ (GJF-2GJ³⁵) since this method is characterized by c₁ = c₃.

Thus, we submit that two of the GJ methods stand out for different purposes. The first, for thermodynamic simulations, is the GJ-I (GJF-2GJ) method, which provides correct thermodynamics in both configurational and kinetic sampling if one uses the half-step velocity in Eq. (44) for kinetic measures. This method allows for simple comparisons between simulations with different friction, since the reference stability limit is fixed for all such values. The second is the GJ-III method, which is characterized by c₃ = 1 such that not only kinetic statistics is accurate but also drift and ballistic velocity measures from the half-step velocity Eq. (44). While this method has a friction-dependent stability limit, an interest in dynamics necessitates relatively small time steps for temporal accuracy. Thus, for studying dynamics, GJ-III may be preferred due to its accurate half-step velocity representation of drift and ballistic behavior.

Coarse-grained (CG) simulations of molecular systems often lead to much larger values of αΔt/m than those that would be seen in the corresponding atomistic simulations. This larger regime of αΔt/m is precisely the area of parameter space where the GJ methods differ significantly from each other (cf. Fig. 2), and where they differ from other methods as well. For this reason, we sought to numerically study the select GJ methods, GJ-I, GJ-II, and GJ-III, using a CG system. For each of the considered methods, a CG polyethylene melt consisting of 128 coarse-grained C₄₈H₉₈ hydrocarbon chains, as pictured in Fig. 3, was simulated for 200 ns at 450 K after an initial NPT equilibration of 10 ns. Accordingly, we chose α/m = 0.1 fs⁻¹ (a large value, by atomistic standards) and used a wide range of choices for the simulation time step Δt, where this range differed for each method due to the method-dependent time scaling. To coarse-grain the polyethylene, the established Shinoda–DeVane–Klein (SDK) coarse-graining model in Ref. 56 was used, resulting in a new CG system with 2048 CG particles and a CG potential energy consisting of bond, angle, and Lennard-Jones force terms. Details on the actual force field parameters can be found in Refs. 51 and 56.

We studied both the translational diffusion rate and the rotational diffusion, or relaxation, time scale by computing the time auto-correlation function (ACF) of the end-to-end vector of each CG polymer chain,

where

is the angle between the end-to-end vector R of a given chain at time t and time t + τ. The brackets ⟨·⟩ indicate averaging among the chains and over trajectory frames t.

In Fig. 4, we display the linear diffusion coefficient as a function of time step, both scaled and non-scaled in order to relate the time step to the stability limit of each method. These data points are obtained from an average of 20 (GJ-II and GJ-III) and 100 (GJ-I and BAOAB) different simulations. In Fig. 5, we present the temporal ACF of P₂(θ) for the different methods with a range of time steps. These curves are each generated as an average from 20 different simulations, and the resulting standard deviation of the averaged curves is estimated to be around 2 × 10⁻³. In both figures, a comparison to the BAOAB method is made as a reference to the importance of the time scale revision that produced the VRORV method from BAOAB. Linear diffusion is measured by the configurational mean-squared distance over time. The diffusive properties observed in Figs. 4–6 are both time step-independent and, overall, the same for the three studied GJ methods. The stability ranges, as indicated by the data and the horizontal arrows, are apparently vastly different. As seen in Fig. 4(b), this difference in stability, however, is accurately compensated by the above-mentioned time scaling factor $c3/c1$⁠, which scales the results of GJ-II and GJ-III almost into concert with those of GJ-I and BAOAB. This scaling is analytically shown in Ref. 36 to enter the stability range. The physical interpretation for the apparent differences between stability of the different methods is therefore mostly an illusion, since the scaling can equally well be interpreted as a change in oscillation frequency, as also shown in Ref. 36. Thus, a method of this kind will always have a stability range in time step that is uniquely linked to the period of motion. It follows that the reason for the GJ-II method to have a stability range significantly larger than GJ-I is that the GJ-II method has slowed its dynamics relative to that of GJ-I, and the mobility per time step (relative to the stability limit) has therefore not changed. This is also the time scaling that links BAOAB and VRORV. Since these GJ results exhibit overall time step independent diffusion results, linear diffusion seems, in these simulations, to be dominated by the mechanism of particles in a viscous medium, which is a feature that the GJ methods are guaranteed to reproduce regardless of the time scale of the dynamics.

For rotational diffusion, we observe in Figs. 5(a)–5(c) that the temporal P₂ correlation is seemingly independent of both time step and the specific GJ method used, even if different time steps can be applied to the different methods. In contrast, Fig. 5(d) shows that the same simulation conducted with the BAOAB method,⁴⁴ which is not specifically designed to yield the correct free particle diffusivity, exhibits a visible time step dependence also in this measure of diffusion. Note here that as the time step is increased, the rotational relaxation time scale decreases; that is, the system’s dynamics speeds up. This is consistent with earlier findings (also displayed in Fig. 4) on the increased linear diffusivity of BAOAB observed for the regime of larger αΔt/m in Refs. 46 and 51. Emphasizing this observation, we notice that the correlation shown in Fig. 5 seems to be dominated by a single time scale for about τ > 10 ns. Extracting this decay rate for data by fitting ⟨P₂(θ(τ))⟩ to

for τ > 10 ns, we obtain Fig. 6 in which the characteristic correlation decay time τ₀ is shown as a function of Δt for the different methods reviewed in this work. As in Fig. 4, where linear diffusion is shown, we display the data for both unscaled and scaled time steps in order to illuminate the similarities between the methods. In concert with Fig. 4, the GJ methods exhibit nearly time step independent decay time of the rotational order across all three methods, while the BAOAB method shows that the orientational order is lost at an increasing rate as the time step is increased. This is consistent with the linear diffusion properties.

Having constructed the set of GJ methods to preserve the Einstein diffusion relation in the case of a free particle, we are not surprised that GJ-I/GJF, which is the GJ method that does not scale the time step as a function of αΔt/m, exhibits translational diffusivity independent of αΔt/m. This was also confirmed previously in Ref. 51. However, given the complexity of the polyethylene system and the highly non-linear interactions between the chains, it is not obvious that the relaxation time scales of the slowest motions should remain unaffected by the implicit time scaling occurring with the GJ-II and GJ-III methods. Instead, Figs. 4–6 show that these relaxation times are controlled only by the friction parameter α/m, indicating that the observed diffusion should be understood as described as an effective viscous medium.

We have conducted a thorough methodical investigation of the possibilities for creating discrete-time, stochastic Størmer–Verlet algorithms that correctly represent the equilibrium statistics of Langevin equations with conservative force fields regardless of the discrete time step. An inspection of the direct integration result, where only the integral over the conservative force is approximated, suggests that a discrete-time algorithm may benefit from two different noise terms per time step, since the noise term in the continuous Langevin equation is integrated differently for the kinetic and configurational coordinates. Following the outline of the analysis that led to the GJ set of methods,³⁶ we conduct a comprehensive analysis of a completely general form of a stochastic Størmer–Verlet expression that includes different noise terms. This is investigated for its ability to yield numerical methods that respond with correct statistical measures of diffusion and Boltzmann distributions when simulating linear (harmonic) systems. The somewhat surprising, albeit definite, result is that this can only be accomplished if the two noise terms are identical. Consequently, the only methods of this kind that can accomplish the basic statistical objectives are the previously identified GJ set of methods, which includes a single noise term per time step,³⁶ and to which discrete-time velocity variables that yield comparable correct kinetic statistics can be associated.

Given that recent developments of methods, starting with the BAOAB method,⁴⁴ have been succeeding in correctly describing the Boltzmann distribution for the configurational coordinate in convex potentials, we have conducted simulations to investigate the more subtle question of time scaling in diffusive behavior in dense, nonlinear systems, exemplified by a coarse-grained polyethylene system. Unlike the BAOAB method, the GJ set also provides the correct diffusion constant as derived from the configurational Einstein definition. However, this diffusion is for a free particle in a trivially flat energy landscape. The dense polyethylene system offers a complex and nonlinear environment in which the accuracy of diffusion of a free particle may not always translate into a reliably calculated diffusion constant by such methods. The results for both linear and rotational diffusion show, however, that the three characteristic GJ methods reproduce the same time step independent diffusion constant throughout their entire stability ranges, leading to the conclusion that the correct free particle diffusion properties of the GJ methods for these simulations translate into robust diffusion properties for more complex systems. In contrast, the BAOAB method exhibits significant deviations as the time step is increased for both linear and rotational diffusion, consistent with this method not having specifically ensured free particle diffusion in its design. It is noted that this discrepancy is purely in the interpretation of the time scale, as illuminated by the VRORV method,⁴⁶ which intersects in its configurational sampling with the GJ-II method. The VRORV method was derived by introducing a time scaling of the A and B steps of the BAOAB method, specifically to address free particle diffusion, and the success of this revision is visible in Figs. 4–6. The introduced time scaling of the methods has some important consequences that have been alluded to in Ref. 36 and that can be seen in Fig. 4 as an apparent difference in stability range of the different methods (the two extremes of the displayed methods being GJ-II and GJ-III). It is important to realize that these differences are not due to actual differences in stability, since this is solely linked to the time scaling $c3/c1$⁠, which affects the dynamics of the simulated behavior. For example, as described in Ref. 36, a damped harmonic oscillator simulated by these methods will have a scaled frequency (as a function of αΔt/m) that ensures all these methods having the same stability limit as measured relative to the period of motion. It is therefore important that one does not infer sampling efficiency or stability directly from Δt but consider the time scaling inherent to each method such that $c3/c1Δt$ is the relevant time step for stability and sampling efficiency. In light of this discussion, we point to BAOAB and GJ-I/GJF, since neither of these are subject to this scaling. Given that GJ-I/GJF-2GJ is the temporally unscaled GJ method (⁠$c3/c1=1$⁠) as a function of friction, we propose that this method generally provides the results with the simplest interpretation of the simulated time scale, especially when using the large time steps offered by these methods.

As mentioned earlier in this paper, the GJ-I/GJF-2GJ method³⁵ is available for use in the LAMMPS simulation suite.^47,48

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This research includes calculations carried out on HPC resources supported, in part, by the National Science Foundation through major research instrumentation Grant No. 1625061 and the U.S. Army Research Laboratory under Contract No. W911NF-16-2-0189, the latter also supporting the work of G.F. The work of B.S. was supported by the National Science Foundation (Grant Nos. DMS–1719640 and DMS–2012271).

We here start with the direct integration of the Langevin equation (1),

where β(s) is given by Eq. (2). This immediately becomes

which for t = t_n+1 = t_n + Δt can be written

with

such that $⟨βAvn⟩=0$ and

Equation (A2) can be written

where $μ(t)=exp(αm(t−tn))$ is the beneficial choice. Integration yields the coordinate r(t),

which can be rewritten