{'fraction_non_alphanumeric': 0.03713714890116347, 'fraction_numerical': 0.012104830179809614, 'mean_word_length': 4.8585197934595525, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Application of the G-JF Discrete-Time Thermostat for Fast and Accurate Molecular Simulations

A new Langevin-Verlet thermostat that preserves the fluctuation-dissipation relationship for discrete time steps, is applied to molecular modeling and tested against several popular suites (AMBER, GROMACS, LAMMPS) using a small molecule as an example that can be easily simulated by all three packages. Contrary to existing methods, the new thermostat exhibits no detectable changes in the sampling statistics as the time step is varied in the entire numerical stability range. The simple form of the method, which we express in the three common forms (Velocity-Explicit, Stormer-Verlet, and Leap-Frog), allows for easy implementation within existing molecular simulation packages to achieve faster and more accurate results with no cost in either computing time or programming complexity.


I. INTRODUCTION
Recently, a new stochastic thermostat [1], based on an exact implementation of the fluctuation-dissipation relationship in discrete time, was presented as integrated into the well-known and widely used Verlet formalism. It was analytically demonstrated that the method provides exact thermodynamic response for both flat and harmonic potentials for any time step within the Verlet stability criteria. This unique thermodynamic feature of the approach, as implemented into the Verlet framework, makes it attractive for a wide range of applications, where it is desired to execute efficient and accurate simulations. Here, we demonstrate that the G-JF method of Ref. [1] is not limited to simple linear cases, but also extends its usefulness to complex nonlinear systems.
Langevin dynamic simulations constitute an appealing approach for simulations of physical systems in contact with a thermodynamic heat bath. A very popular class of such systems is molecular dynamics (MD) [2] of classical particle ensembles. The method is based on numerical integration of the Langevin equation where r is the coordinate, v =ṙ is its velocity, m is its mass, f is a net deterministic force acting on r. The friction constant α > 0 and the noise β are connected by the fluctuation-dissipation relationship [3] β(t) = 0 (2) β(t)β(t ′ ) = 2αk B T δ(t − t ′ ) with k B and T being Boltzmann's constant and the thermodynamic temperature, respectively. Given the appli-cability of this equation to a wide spectrum of physical systems with phenomenological dissipation and thermal noise, there have been decades of development in numerical methods for solving such equations. Specifically, the perfection of the most important thermodynamic properties of discrete-time numerical simulations have been of particular interest. We here point to the vast literature through Refs. [1,[4][5][6][7][8][9][10][11] and references therein. The most commonly sought after properties in stochastic simulations have been i) diffusion of a particle in a flat potential, ii) transport on a linear ramp potential (which can be mapped directly onto the diffusive behavior in a flat potential), and iii) Boltzmann sampling in harmonic potentials [12]. Since the G-JF method [1] exhibits all these features thermodynamically correct in discrete time, we here wish to provide additional expressions for practical implementation of the method. We also provide a simple, yet representative, example of how the method performs for a nonlinear and complex system in comparison to several widely used contemporary molecular dynamics simulation suites.

II. THE THREE VERLET EXPRESSIONS
Verlet integrators are usually expressed in one of three forms [13]: A) velocity-explicit Verlet (VE), which advances the trajectory one time step based on the coordinate and its conjugate variable (here the velocity), B) Störmer-Verlet (SV), which uses coordinates at two consecutive time steps to advance time, and C) leap-frog (LF), which advances the trajectory based on the coordinate and its conjugate variable, the latter being rep-resented at half-integer time steps relative to the coordinate. The three typical Verlet formulations produce identical trajectories, and, thus, are different expressions of the exact same method. Consequently, applications of the Verlet method are commonly expressed in any of the available forms. Since all three identical methods are frequently used in the literature, we start here with the recently published G-JF thermostat that was derived in the natural VE form, and re-express the algorithm in the two other popular forms. The previous harmonic oscillator analysis [1] of the method applies to all three variants since they produce identical results. Specifically, the three following formulations will result in the correct fluctuation-dissipation relationship and thus, the correct thermodynamic response in linear systems.

A. Velocity Explicit G-JF
We take our starting point with the VE G-JF expressions that were derived and analyzed in [1]. Denoting discrete time variables by the integer time step superscript, such that, e.g., r n = r(t n ), the algorithm for advancing r n and v n one time step of dt reads and where is a standard Gaussian random number that satisfies Setting the initial conditions (r 0 , v 0 ), Eqs. (4)-(8) can be directly used to generate the trajectory (r n , v n ) from which the dynamical and statistical information can be derived.

B. Störmer-Verlet G-JF
We here start by rewriting Eqs. (4) and (5) for n → n − 1: As outlined in Ref. [1], Eq. (10) can be inserted into Eq. (4) in order to replace v n , whereafter Eq. (9) is used to replace the resulting v n−1 . This yields which is the SV formulation of the G-JF method. Unlike the VE expressions, the SV equation does not contain direct information about the velocity and is therefore not directly applicable for natural initial conditions (r 0 , v 0 ). The self-consistent approach for starting this procedure from (r 0 , v 0 ) is to apply Eq. (9) for n = 0, then apply Eq. (11) for all subsequent time steps n > 0.
In order to calculate both the complete dynamical trajectory and important thermodynamic quantities, one needs the velocity v n explicitly expressed as well. Consistent with the method, we replace v n−1 in Eq. (10) by inserting Eq. (9), such that Thus, Eqs. (11), (12), and (13) constitute the identical SV form of the VE expressions.

C. Leap-Frog G-JF
This version of the Verlet method for Langevin equations comes with some flexibility in how the method is expressed. We take the starting point with the SV form given in Eq. (11) and introduce a reasonable definition of the half-step velocity We now use Eq. (14) to replace r n+1 in Eq. (11) to obtain in which we can again apply Eq. (14) and arrive at This equation is the half-time step velocity propagator, which is complemented by Eq. (14) to yield a LF G-JF method As was the case for the SV formulation of the method, this LF representation does not trivially incorporate the natural initial conditions (r 0 , v 0 ). We therefore apply Eq. (14) for n = 0 in combination with Eq. (12), resulting in v 1 2 to be used, along with Eq. (12), before applying Eqs. (16) and (17) for n > 0. The proper integer-step velocity can be found from combining Eqs. (16) and (4), resulting in where we can then, again, use Eq. (14) to obtain Thus, Eqs. (16), (17), and (20) constitute a consistent LF method, where initial conditions (r 0 , v 0 ) are applied through Eqs. (12) and (18). Notice that while the given SV method is uniquely connected to the VE expressions in both the coordinate {r n } and its evaluated velocity {v n }, the development of LF is not unique, even if LF is constructed to produce identical trajectories (r n , v n ) to those of VE and SV. The reason is the aforementioned somewhat ambiguous choice in defining the half-step velocity v n+ 1 2 in Eq. (14), and the subsequent reconstruction of the integer-time velocity v n in Eq. (20). Thus, one can exercise some freedom of choice in the velocity equations when using the LF expressions.
For example, a sensible alternative to the half-step velocity in Eq. (14) may be where we use the symbol u for the revised definition of the half-step velocity. Following the procedure starting from Eq. (14), we can develop the following alternate LF formulation, which also yields identical trajectories (r n , v n ). The two equations for half-step velocity and integer-step position become with initiating half-step velocity u 1 2 and conversion to integer-step velocity v n written respectively. As mentioned above, all LF formulations of the thermostat yield identical trajectories for (r n , v n ) and therefore constitute the same method regardless of the specific definition of the half-step velocity. Choosing which one to use is entirely a matter of convenience. The defined half-step velocity is simply an auxiliary variable, which need not have any specific useful physical interpretation for the evaluation of thermodynamic quantities.

III. TESTING THE METHOD
We exemplify the applicability of the G-JF method in the context of a simple biomolecular simulation and we make comparison to results of the widely-used molecular dynamics codes AMBER [14], LAMMPS [16,17], and GROMACS [7]. These packages employ different commonly used stochastic thermostats. More comprehensive discussions of other thermostats can be found in recent references [1,4,5]. The purpose of the following simulations is not to present new or ground-breaking results in biomolecular science; instead, we choose a simple and well-understood representative model system that can illuminate, through simulations within well-established MD packages, some key features of the new algorithm as implemented into one of the codes that is particularly amenable to revisions.

A. Simulation Details
We performed classical molecular dynamics simulations of alanine dipeptide (illustrated in Figure 1a), a small and well-studied biomolecule. Intramolecular interactions among the solute atoms, including bond, angle, dihedral, and non-bonded energies, were modeled with an AMBER classical all-atom molecular mechanical treatment [14], coupled with the recent ff12 parameter set [15], while the extra-molecular environment was treated as vacuum.
This model system was simulated at a target temperature of 300 K with four integration schemes for comparison: i) The BBK thermostat [10], which is implemented and expressed within AMBER 12 (ntt =3) in the leapfrog Verlet formulation [8]; ii) the method of Schneider & Stoll [11], as implemented in the LAMMPS simulation software [16] (with the combination of the nve and langevin fixes); iii) a variation of the method of van Gunsteren & Berendsen [9,18], implemented in GROMACS [7] (as the sd integrator [19]), which includes velocity rescaling; and finally iv) the G-JF thermostat [1].
The latter was implemented into AMBER through small modifications to the AMBER 12 source code [15]. Since the underlying LF Verlet integrator used in the BBK method is consistent with Eqs. (16), (17), and (20), revisions were only necessary with regard to the fluctuation terms. The memory framework was modified to include space for the correlating noise term from the pre-vious time step, corresponding to β n in Eqs. (16) and (20), to be cached for use in subsequent integration steps. A Langevin damping coefficient of 10 ps −1 was used throughout all the simulations. Simulation input data for LAMMPS and GROMACS runs were generated directly from AMBER-formatted coordinate and topology files using the free tools amber2lammps.py (distributed with the LAMMPS source code), and amb2gmx.pl [20], respectively.
Data were collected from sets of ten independent 100 ps simulations at several values of the integration time step, from 0.5-3.1 fs in increments of 0.1 fs, and from 3.1-3.2 fs in increments of 0.01 fs. Each simulation was initiated from the same energy-minimized starting structure of the peptide. For each independent simulation, initial atomic velocities were assigned from a Maxwell-Boltzmann distribution appropriate to the target temperature and the psuedo-random number chain was uniquely seeded.
First and second moments of the total potential energy distribution were chosen as measures of the stability and accuracy with which each integration method could reproduce the statistical-mechanical properties of the physical system. Average energy and its fluctuation were computed based on samples at all integration steps. The total number of integration steps N was varied based on the time step dt, such that N = ⌊100 ps/dt⌋.

B. Results
Figures 1b and c summarize the results and shows the average total potential energy and its standard deviation as a function of the simulated time step for all the integration methods. As expected, all methods give very similar results for both the average and standard deviation of the energies for as long as the time steps are smaller than about 1 fs.
However, for increasing time steps the unmodified contemporary codes start deviating from the expected statistical values found for small dt. The two strongest deviations are found for the BBK method, implemented in AMBER 12, and the thermostat of Schneider & Stoll, implemented into LAMMPS. The same deviating behavior, albeit seemingly with only one third deviation, is found for the thermostat of van Gunsteren & Berendsen, implemented into GROMACS. In contrast, the G-JF method, implemented into AMBER 12 as described above, shows the inherent feature of preserving the fluctuation-dissipation relationship for any time step. In fact, for these simulations, the G-JF method can give statistical estimates at all time steps up to the stability limit that are indistinguishable from those at small time steps. The stability limit is here identified by adiabatically increasing the time step of a simulation until simulations produce anomalously high velocities, indicating that the Verlet integrator is no longer capable of producing a meaningful trajectory for the simulation. The existing stochastic thermostats seem stable up to 3.0 fs for this simulation, while the G-JF is stable for up to a slightly higher value 3.1 fs.
We close this section by arguing for the use of potential energy (and its fluctuations) as a measure for the ability of the different thermostats to correctly sample the phase space corresponding to the target temperature. A seemingly more straightforward test would be to consider the kinetic energy. However, as shown in Ref. [1] for the harmonic oscillator case, the computed velocity v n , and thereby the kinetic energy, is increasingly depressed for increasing time step in any Verlet scheme, regardless of the inclusion of a thermostat. Thus, when evaluating a thermostat, measures of velocity and kinetic energy may not be appropriate for determining the quality of statistical sampling. Further, this observation may hint at statistical sampling problems arising from thermostats that involve "velocity rescaling" as a way to obtain a desired temperature.

IV. DISCUSSION AND CONCLUSION
The newly developed G-JF stochastic Verlet thermostat has been applied to simulations of a small, yet nontrivial and nonlinear, system representative of many applications in molecular modeling. It has previously been analytically demonstrated that for linear systems the new method yields exact statistical behavior of diffusion and Boltzmann distributions for any time step leading to stable dynamics [1]. The simulations presented here indicate that these attractive and robust statistical features are likely to remain in other complex and nonlinear systems. Specifically, our present simulations demonstrate that the statistical behavior of potential energy remains sound for any time step up to the limit where the atomic trajectories suddenly diverge. In contrast, available contemporary molecular dynamics codes with other stochastic thermostats show deviating behavior for increasing time step, indicating that the interpretation of thermodynamic data from those algorithms must be done with caution and small time steps. This seems to be the case for the three popular molecular dynamics codes that we have investigated here, and is likely to be true also for other available MD simulation codes. In order to facilitate comparison between the stochastic thermostats, we use the same AMBER force fields in all four sets of simulations. We emphasize that simulations with the G-JF thermostat have been completed by implementing the new simple algorithm into an existing available code (AMBER 12), to ensure direct comparison between the thermostats without any other differences in parameters or simulation details. Thus, based on the preliminary tests and analyses, we suggest that the algorithm presented here and in Ref. [1] be implemented into existing molecular dynamics codes for further use and evaluation. In order to facilitate such revisions, we have explicitly provided the algorithm in all three commonly used formulations of the Verlet method.