Accurate Configurational and Kinetic Statistics in Discrete-Time Langevin Systems

TL;DR

This paper introduces a new velocity definition for the GJF thermostat in Langevin systems, enabling accurate, time-step-independent kinetic and configurational statistics in discrete-time simulations, including molecular dynamics.

Contribution

A novel velocity variable for the GJF thermostat that provides correct kinetic responses and enables a new Leap-Frog algorithm for accurate statistical sampling.

Findings

The new velocity yields correct kinetic energy and fluctuations.

The algorithm is stable and accurate for nonlinear systems.

It maintains time-step independence in complex molecular simulations.

Abstract

We expand on the previously published Gr{\o}nbech-Jensen Farago (GJF) thermostat, which is a thermodynamically sound variation on the St{\o}rmer-Verlet algorithm for simulating discrete-time Langevin equations. The GJF method has been demonstrated to give robust and accurate configurational sampling of the phase space, and its applications to, e.g., Molecular Dynamics is well established. A new definition of the discrete-time velocity variable is proposed based on analytical calculations of the kinetic response of a harmonic oscillator subjected to friction and noise. The new companion velocity to the GJF method is demonstrated to yield correct and time-step-independent kinetic responses for, e.g., kinetic energy, its fluctuations, and Green-Kubo diffusion based on velocity autocorrelations. This observation allows for a new and convenient Leap-Frog algorithm, which efficiently and precisely represents statistical measures of both kinetic and configurational properties at any time step within the stability limit for the harmonic oscillator. We outline the simplicity of the algorithm and demonstrate its attractive time-step-independent features for nonlinear and complex systems through applications to a one-dimensional nonlinear oscillator and three-dimensional Molecular Dynamics.

Equations

$$m\ddot{r}$$

$$r^{n+1}$$

$$v^n$$

$$r^{n+1}$$

$$v^{n+1}$$

$$v^{n+\frac{1}{2}}$$

$$r^{n+1}$$

$$v^n$$

$$m\dot{v}+\alpha \dot{r}$$

$$⟨\beta (t)⟩$$

$$⟨\beta (t)\beta (t^{\prime })⟩$$

$$r^{n+1}$$

$$v^{n+1}$$

$$a$$

$$b$$

$${\beta }^{n+1}$$

$$⟨{\beta }^n⟩$$

$$⟨{\beta }^n{\beta }^l⟩$$

$$r^{n+1}$$

$$v^n$$

$$⟨E_p⟩$$

$$⟨E_k⟩$$

$$⟨E_k^{n+\frac{1}{2}}⟩$$

$$=$$

$$⟨(r^{n+1}{)}^2⟩$$

$$\frac{k_{B}T}{\kappa }$$

$$⟨r^n{r}^{n-1}⟩$$

$$⟨E_k^{n+\frac{1}{2}}⟩$$

$$u^{n+\frac{1}{2}}$$

$$⟨E_k^{n+\frac{1}{2}}⟩$$

$$u^{n+\frac{1}{2}}$$

$$r^{n+1}$$

$$u^{n+\frac{1}{2}}$$

Full text

Submitted to Molecular Physics November 1, 2018

Accepted for publication January 7, 2019.

I Introduction

Molecular Dynamics (MD) offers the possibility for simulating a set of $N$ interacting particles over a representative time period, and thereby producing trajectories that can be used to generate statistical and thermodynamic properties AllenTildesley ; Frenkel ; Rapaport . This is accomplished by formulating a set of Newtonian equations of motion of the form

[TABLE]

where $m$ is the mass of a particle with coordinate $r$, which is subjected to the force $f$ that represents the particle’s interaction with its surroundings. The discrete-time numerical approximation to the solution of Eq. (1) is most often acquired by methods rooted in the Størmer-Verlet (SV) Stormer_1921 ; Verlet ; Gear algorithm that evolves the coordinate $r(t_{n})=r^{n}$ at time $t_{n}$ to the next time step, $dt$ later, at $t_{n+1}=t_{n}+dt$ when the coordinate is $r^{n+1}=r(t_{n+1})$. The second order (in $dt$) discrete-time equation is

[TABLE]

where $f^{n}=f(r^{n},t_{n})$. This method is appealing for a number of reasons, including its simplicity, efficiency, stability properties, and time reversibility, which implies attractive conservation properties for the trajectory $r^{n}$ for closed systems ($f=f(r)$). A simple model system that directly illuminates the features of this algorithm is the harmonic oscillator, for which the SV method Eq. (2) can produce a discrete-time harmonic oscillator trajectory for as long as the reduced time step $dt\Omega_{0}<2$, where $\Omega_{0}$ is the natural frequency of the oscillatorVenneri ; Pastor_88 . While the SV method will reproduce a perfect harmonic oscillator, the discrete-time oscillation frequency $\Omega_{V}$ will be enhanced beyond the continuous-time frequency $\Omega_{0}$ (see, e.g., Ref. GJF1 and references therein, or the Appendix in this paper).

Given the attractive properties of the SV discrete-time trajectory $r^{n}$, it is desirable to also approximate the accompanying velocity. This, however, introduces a new set of errors since discrete-time does not offer a proper definition of a differential, which is clearly necessary for evaluating the velocity. The two most obvious approximations are given by second order central differences. One of these AllenTildesley ; Swope ; Beeman defines an on-site velocity $v^{n}$ by

[TABLE]

which, combined with Eq. (2), can be expressed as the Velocity-explicit Verlet (VV) method

[TABLE]

The other approximation AllenTildesley ; Buneman ; Hockney defines a half-step velocity $v^{n+\frac{1}{2}}$ by

[TABLE]

which, combined with Eq. (2), can be expressed as the Leap-Frog (LF) method

[TABLE]

While the trajectories of Eqs. (2), (4), and (8) are identical, the associated velocity variables are different as can be seen from the following relationship

[TABLE]

Of course, both velocity definitions limit the true velocity $v=\dot{r}$ for $dt\rightarrow 0$.

Having highlighted the attractive features of the Størmer-Verlet method above, it is important to emphasize that the discrete-time velocity variables $v^{n}$ and $v^{n+\frac{1}{2}}$ are fundamentally inconsistent with $r^{n}$; i.e., the discrete-time velocities do not precisely represent the conjugated variable of the simulated trajectory. This has previously been indicated in various ways (see, e.g., Refs. vgb_1982 ; Pastor_88 ; holian95 ; GJF1 ; GJF3 ) through harmonic oscillator analysis. Thus, while $r^{n}$, $v^{n}$, and $v^{n+\frac{1}{2}}$ are all second order approximations to the true values $r(t_{n})$ and $v(t_{n})$, the mutual inconsistency between the coordinate and its velocity results in the discrete-time momentum not being the exact conjugated variable to the simulated trajectory. The harmonic oscillator illuminates this clearly, and the appendix of this paper further emphasizes this critical issue by including the LF half-step velocity into the analysis, indicating that properties extracted from Verlet methods will generally show inconsistencies between kinetic and configurational measures in the same simulation.

The SV, VV, and LF methods allow for simulations of the constant energy, or $(N,V,E)$, micro-canonical ensemble, whereas the constant temperature $(N,V,T)$ ensemble is often more relevant. The well-known governing equation for this type of dynamics is given by the Langevin equation Langevin

[TABLE]

where linear friction is represented by the constant $\alpha\geq 0$, which is related to the thermal fluctuations $\beta$, which can be chosen to be represented by the Gaussian distribution Parisi

[TABLE]

where $k_{B}$ is Boltzmann’s constant and $T$ is the thermodynamic temperature.

Many methods for controlling the temperature of a simulated system ( thermostats ) have been developed, and most of them fall into two major categories: Deterministic (e.g., Nosé-Hoover Nose ; Hoover ) and stochastic (Langevin) thermostats SS ; BBK ; Pastor_88 ; Vanden ; ML ; Paquet . The deterministic approach includes additional degrees of freedom, which act as an energy reservoir and thereby mimic a thermal heat bath. A requirement for such method is that the temperature of a simulated system can be reliably measured in order for the method to interact properly with the heat-bath. The stochastic approach is to directly simulate the Langevin equation Eq. (10), which does not include additional degrees of freedom, but instead interacts with a heat-bath through the fluctuation-dissipation balance for $\alpha>0$. Since discrete-time tends to distort the conjugated relationship between the positional coordinate and its corresponding momentum (see, e.g., Appendix in Ref. GJF3 and the Appendix in this paper), a common problem for all methods is that kinetic and configurational measures of temperature disagree, which is a concern for both the integrity of a simulation and the extraction of self-consistent information, which may depend on configurational as well as kinetic sampling. It is therefore imperative to understand how to properly define a kinetic measure consistent with the statistics of the trajectory.

This work will consider the stochastic Langevin approach, which was pioneered by, e.g., Refs. SS ; BBK ; Pastor_88 ; Gunsteren by including discrete-time friction and noise into the Verlet framework. As clearly analyzed in Ref. Pastor_88 , the inconsistencies between kinetic and configurational properties persist into the statistical averages obtained from the discrete-time thermostats. Consistent with the harmonic oscillator analysis of the algorithms, the configurational temperature is typically found to be higher than the measured kinetic temperature. One method was identified in Ref. Pastor_88 to provide the correct time-step independent kinetic energy for a harmonic oscillator, but the inconsistency with the configurational properties made the authors reason that this was not desirable. Instead, their analysis revealed that while another method BBK displayed significant time-step errors in energies, these errors were the same for kinetic and configurational measures. This consistency in error made the authors recommend using the stochastic BBK method BBK . Since then, much work has been devoted to the development of methods that minimize the time-step error on statistical measures (see, e.g., Refs. ML ; Goga_2012 ; vdSpoel_2005 ; GJF1 ). A major development of accomplishing correct sampling of configurational statistics for the harmonic oscillator was reported in 2012 ML . This method, however, did not provide correct configurational diffusion in discrete time. We here focus specifically on the GJF method GJF1 , since that method has been analytically shown to provide exact configurational thermodynamic response for both flat and harmonic potentials for any time step within the Verlet stability criterion. The method has further been demonstrated to give extremely robust configurational sampling in complex molecular dynamics environments of both soft and hard materials GJF1 ; GJF2 ; GJF3 ; GJF4 , and it has been implemented into the MD suite LAMMPS Plimpton ; LAMMPS-Manual as an option for a stochastic thermostat. As expected, the kinetic measures of, e.g., kinetic energy (temperature) are, however, not reliable due to the time-discretization. We here seek to address this fundamental problem by redefining the discrete-time velocity variable such that the revised GJF method can provide reliable statistical information simultaneously for configurational as well as kinetic properties.

II Discrete-Time Langevin Dynamics

We adopt the well-documented GJF method that was derived in Ref. GJF1 and further analyzed for complex applications Ref. GJF2 ; GJF4 . The method reads

[TABLE]

where

[TABLE]

The discrete-time noise is

[TABLE]

which results in an uncorrelated Gaussian random number with zero mean and a variance given by the temperature and friction coefficient:

[TABLE]

Notice that the special case $\alpha=0$ reduces the method to the VV form of the standard Størmer-Verlet algorithm, Eqs. (4) and (5). The derivation GJF1 of the GJF method was conducted by directly integrating the Langevin equation Eq. (10) over a time-step. By taking advantage of the normal distribution of the fluctuations and the exact integral of the friction term only as a function of spatial displacement, this approach maintains the fluctuation-dissipation relationship in discrete time. Thus, the known inaccuracies of velocity are avoided in the fluctuation-dissipation balance, and the configurational Boltzmann statistics remains accurate even for large time steps within the stability range (see, e.g., Eq. (22) below). As outlined in Ref. GJF2 , Eqs. (13) and (14) can be combined to give the trajectory in the SV form, where the trajectory $r^{n}$ is separated from the associated velocity $v^{n}$ by linear transformation. The result is:

[TABLE]

where the GJF velocity Eq. (14) can be expressed directly from the trajectory by note_on_typo

[TABLE]

Thus, Eqs. (LABEL:eq:gjf_sv) and (II) can be recombined to yield Eqs. (13) and (14), implying that the trajectory Eq. (LABEL:eq:gjf_sv) is identical to what results from Eqs. (13) and (14).

We will now give a new definition of a useful velocity parameter, based on a linear analysis of the GJF trajectory, Eq. (LABEL:eq:gjf_sv).

II.1 Linear Analysis, $f=-\kappa r$

The introduction of the GJF algorithm GJF1 was rooted in the remarkable exact features for the linear Langevin equation (10) for $f=-\kappa r$, $\kappa>0$, and $\alpha>0$:

[TABLE]

with the additional observations that for $f=0$ ($\kappa=\Omega_{0}=0$) and $\alpha>0$, the method yields the exact Einstein diffusion $D=k_{B}T/\alpha$, and the correct kinetic energy given by Eq. (23) for $\Omega_{0}=0$. These results are true for all $\Omega_{0}dt\leq 2$.

Equation (23) displays the features of the imperfect velocity given in Eq. (14). Thus, while Eq. (22) is always an excellent measure of the configurational sampling, (23) is an increasingly poor representation of the kinetic temperature for increasing time-step $dt$. Of course, since the error factor in the measured kinetic energy is explicitly given, one is tempted to redefine the velocity in order to obtain the correct kinetic energy (temperature). However, that will only work for the specific harmonic mode with frequency $\Omega_{0}$ and a velocity correction of this kind can therefore not be successful for more complex or nonlinear problems.

Inspired by Ref. Pastor_88 , where more than one velocity definition was applied to a given trajectory, we evaluate the kinetic energy of the GJF trajectory Eq. (LABEL:eq:gjf_sv) based on the half-step velocity definition $v^{n+\frac{1}{2}}$ given in Eq. (6):

[TABLE]

where we have used $\langle(r^{n+1})^{2}\rangle=\langle(r^{n})^{2}\rangle$.

For the harmonic oscillator, where $f=-\kappa r$, $\langle(r^{n})^{2}\rangle$ is known for the GJF trajectory to be given by Eq. (22), $\langle r^{n+1}r^{n}\rangle$ can be found from Eq. (LABEL:eq:gjf_sv) as follows: Inserting $f=-\kappa r$ into Eq. (LABEL:eq:gjf_sv) yields the autocorrelation

[TABLE]

Since $\langle r^{n}\beta^{n}\rangle=\langle(\beta^{n})^{2}\rangle bdt/2m$ (see Eq. (LABEL:eq:gjf_sv) for $\langle r^{n+1}\beta^{n+1}\rangle$), we have from Eqs. (19) and (22) that

[TABLE]

which, after some reduction, yields

[TABLE]

Finally, inserting Eqs. (22) and (28) into Eq. (25) gives

[TABLE]

This is a useful result since it is independent of $\Omega_{0}dt$. Thus, for any dynamical mode, we can for the GJF method propose the following re-assigned definition of discrete-time velocity:

[TABLE]

For the noisy harmonic oscillator, this 2GJ definition of velocity results in the average kinetic energy

[TABLE]

Thus, we now have self-consistency between time-step independent configurational and kinetic statistics for the noisy harmonic oscillator, and we have the correct statistical values of both potential and kinetic energies.

II.2 Specific Verlet-type algorithms for using the 2GJ velocity

The linear analysis above suggests that the Langevin equation with constant friction coefficient $\alpha$, and a general force of the kind $f=f(r)$, can be thermodynamically analyzed with some improvements by incorporating the re-assigned velocity $u^{n+\frac{1}{2}}$ into the GJF algorithm.

Direct applications of the accompanying velocity $u^{n+\frac{1}{2}}$ is straightforward for the GJF method in the SV form Eq. (LABEL:eq:gjf_sv), since statistics of kinetic measures are easily conducted by $u^{n+\frac{1}{2}}$ from Eq. (30). Similarly, the application of $u^{n+\frac{1}{2}}$ is straightforward for the GJF equations in the VV form Eqs. (13) and (14), where the on-site GJF velocity $v^{n}$ is used for facilitating the calculation of $r^{n}$, and $u^{n+\frac{1}{2}}$ is used for calculating the kinetic properties.

Since the new measure of velocity is defined at the half-step, we can write an LF version of the GJF trajectory to directly take advantage of the 2GJ definition. Combining Eq. (30) with the SV expression, Eq. (LABEL:eq:gjf_sv), we obtain the GJF-2GJ algorithm

[TABLE]

or, equivalently

[TABLE]

where Eq. (35) is simply the 2GJ velocity definition from Eq. (30). Thus, one can express the method we are proposing here in many different convenient ways. However, the result is always the same; the hybrid between the GJF trajectory and the re-assigned 2GJ half-step velocity presented above. For $\alpha=0$, the method reduces to the standard LF form given in Eqs. (7) and (8).

II.3 Diffusion for $f=0$

As mentioned above, the GJF trajectory has been shown to give the correct diffusion as calculated from the configurational Einstein definition of the diffusion coefficient $D$ (for $f=0$):

[TABLE]

This is the true diffusion coefficient, since the measure considers the actual square displacement as a function of discrete time. The corresponding continuous-time kinetic Green-Kubo GreenKubo expression for the same diffusion coefficient, calculated from

[TABLE]

where $v(t)$ is the continuous-time velocity.

Using the 2GJ velocity $u^{n+\frac{1}{2}}$ Eq. (32) for $f=0$ we have the thermodynamically limiting form

[TABLE]

for $|a|^{q}\rightarrow 0$, which represents the thermodynamic limit (infinite volume) for $t\rightarrow\infty$, where all inertial information about the initial condition $u^{\frac{1}{2}}$ has been lost. The velocity auto-correlation for equilibrated initial conditions ($|a|^{q}\ll 1$) is then

[TABLE]

There are two discrete-time approximations when adapting the Green-Kubo expression (37) to the discrete-time velocity correlation function. One is the accuracy of the velocity variable, which we have emphasized above. The other is the discretization of the $ds$ integral, since a non-zero $dt$-discretization allows the use of any Riemann sum for as long as that sum converges to the correct result when $dt\rightarrow 0$. Thus, in discrete-time there are a multiple of different results for the Green-Kubo calculation of diffusion for each definition of the velocity. Choosing the right-Riemann sum for approximating the integral Eq. (37) we get

[TABLE]

which clearly shows the correct Green-Kubo diffusion coefficient for the 2GJ velocity definition for any time step $dt$. Notice that other approximations of the integral may give different results for nonzero $dt$. For example, the left-Riemann sum yields $D_{lR}=k_{B}T(\frac{1}{\alpha b}+\frac{dt}{m})$ and the trapezoidal sum gives $D_{tr}=k_{B}T/\alpha b$. While all possible Riemann approximations give the correct limit for $dt\rightarrow 0$, only the right Riemann sum is independent of $dt$. We observe that if the GJF trajectory is paired with the traditional LF half-step definition of velocity, as given by Eq. (6), the correct time-step independent result for the Green-Kubo expression for diffusion is obtained for the trapezoidal approximation to the integral Eq. (37). However, that choice of defining the velocity would yield the incorrect ($dt$ dependent) result Eq. (29) for the average kinetic energy in a harmonic potential. Thus, for the 2GJ velocity definition Eq. (30), we define the right-Riemann sum to be the proper discrete-time approximation to the Green-Kubo integral, and we have then established consistency between kinetic and configurational measures, both for transport ($f=0$) and for the sampling of distributions and energies when $f$ represents a linear Hooke’s force.

We note parenthetically that the case of $f=constant$ is entirely equivalent to the special case $f=0$ above, if one adds a constant velocity drift term $f/\alpha$ to the solution.

III Numerical Simulations

The above analysis and results show that the proposed algorithm gives the correct configurational and kinetic statistical sampling of a noisy harmonic oscillator regardless of the applied time step. It also shows that configurational and kinetic measures are correct and consistent for any time step when the potential is flat. In order to verify the usefulness of the 2GJ velocity definition as a companion to the GJF trajectory in nonlinear and complex systems, we show the results of characteristic simulations both for a one-dimensional nonlinear oscillator and for representative Molecular Dynamics examples.

III.1 One-dimensional nonlinear oscillator

We consider a one-dimensional oscillator with the deterministic force $f$ given by the potential energy surface $E_{p}(r)$

[TABLE]

where $E_{p}$ is chosen to be a nonlinear convex potential

[TABLE]

for $|r|<r_{0}$. The characteristic energy and distance are denoted $E_{0}$ and $r_{0}$, respectively. We define the characteristic time scale from the inverse of the characteristic frequency $\omega_{0}=\sqrt{E_{0}/m}/r_{0}$, which yields the small amplitude oscillation frequency $\Omega_{0}=\sqrt{E^{\prime\prime}_{p}(0)/m}=\sqrt{60}\,\omega_{0}$, and the characteristic velocity $v_{0}=r_{0}\omega_{0}$.

Simulations are conducted for thermodynamic temperatures $k_{B}T=0.1E_{0}$ and $k_{B}T=0.5E_{0}$, and for friction coefficients $\alpha=1\,m\omega_{0}$ and $\alpha=5\,m\omega_{0}$. Statistical averages are acquired over 1010 time steps. The statistical measures are as follows: Configurational and kinetic temperatures are defined as

[TABLE]

where $T_{c}$ is the configurational temperature Hirschfelder ; Rickayzen and $T_{k}$ is the kinetic temperature from Eq. (31). We further evaluate the variances $\sigma^{2}_{p}$ and $\sigma_{k}^{2}$ of the potential and kinetic energies, respectively:

[TABLE]

Figure 1 shows the simulation results of these measures as a function of the reduced time step $\Omega_{0}dt$ for the four different parameter combinations of $T$ and $\alpha$ given above. The markers show the statistical data, while the horizontal lines indicate the correct theoretical results for $T_{c}=T$ and $\langle E_{k}\rangle=\frac{1}{2}k_{B}T_{k}=\frac{1}{2}k_{B}T$, and $\sigma_{k}=\frac{1}{\sqrt{2}}k_{B}T$. The inserted table indicates the symbols that represent the different parameter combinations and quantities. The arrows indicate the approximate stability limits for the four simulated parameter combinations. These limits are found by starting a simulation for small $\Omega_{0}dt$, then conducting 1010 time steps before increasing $\Omega_{0}dt$ until the simulated particle displays instability by exceeding the bounds of the potential; i.e., the system is defined unstable when $|r|>r_{0}$. The stability limit results are: $\Omega_{0}dt<0.215$ ($k_{B}T=0.5E_{0}$ and $\alpha=0.5m\omega_{0}$), $\Omega_{0}dt<0.255$ ($k_{B}T=0.5E_{0}$ and $\alpha=5m\omega_{0}$), $\Omega_{0}dt<0.405$ ($k_{B}T=0.1E_{0}$ and $\alpha=0.5m\omega_{0}$), $\Omega_{0}dt<0.525$ ($k_{B}T=0.1E_{0}$ and $\alpha=5m\omega_{0}$). The different stability limits are due to the nonlinearity, which has a diverging curvature for increasing $|r|$. What is clear from the figure is that all statistical properties are nearly independent of the time step for as long as the stability limit has not been exceeded. Additionally, the known exact statistical values (indicated by horizontal lines and on the right legends) are very well matched by the simulations. None of the acquired data shows any deviating trends for increasing $\Omega_{0}dt$, which indicates that even the measures of potential energy $\langle E_{p}\rangle$ and its fluctuations $\sigma_{p}$ are correctly sampled, as has previously been demonstrated GJF2 ; GJF4 . We notice that simulations conducted for very low temperatures must result in a stability range $\Omega_{0}dt<2$, since the configurational sampling would effectively be happening in the harmonic limit.

Encouraged by the robust statistical averages exemplified in Figure 1, we explore if the sampling distributions of configurational coordinate $r$ and the 2GJ velocity $u$ match our expectations. We exemplify the results by showing data for $\alpha=5m\omega_{0}$ and $k_{B}T=0.5E_{0}$ for $\Omega_{0}dt=0.25$; i.e., the highest time step for which the system is stable. For this extreme time step we produce the configurational and kinetic density distributions $\rho_{c}(r)$ and $\rho_{k}(u)$ over the 1010 simulated time steps. Ideally, these distributions relate to the potential and kinetic energies through the Maxwell and Maxwell-Boltzmann distributions such that

[TABLE]

from where we can generate the effective potential and kinetic surfaces (potentials of mean force):

[TABLE]

The constants $C_{c}$ and $C_{k}$ are determined such that, e.g., $E_{\rm PMF}(0)=E_{p}(0)$ and $U_{\rm MB}(0)=E_{k}(0)=0$. If the sampled discrete-time statistics is correct and consistent with the Boltzmann and Maxwell-Boltzmann distributions, we have the ideal comparisons: $U_{\rm PMF}(r)=E_{p}(r)$ and $U_{\rm MB}(u)=E_{k}(u)=\frac{1}{2}mu^{2}$.

Figure 2 shows the effective potentials $U_{\rm PMF}(r)$ and $U_{\rm MB}(u)$ for the parameters and procedures described above. Fig. 2a shows the simulation result $U_{\rm MB}(u)$ as a solid curve with the correct function $\frac{1}{2}mu^{2}$ shown as a dashed curve. It is obvious from this comparison that the 2GJ velocity is not only giving the correct, time-step independent statistical averages shown in Figure 1, but also has the correct quadratic distribution even for this extreme value of the time step $\Omega_{0}dt$. Figure 2b shows the comparable potentials of mean force $U_{\rm PMF}(r)$ as a solid curve with a comparison to the exact behavior $E_{p}(r)$ shown with a dashed curve. Again, this comparison is nearly flawless even for this very high time step. Included in the figure are the simulation results of a harmonic oscillator with the same parameters as for the nonlinear oscillator described above. We have chosen an oscillator with $\kappa=E_{p}^{\prime\prime}(0)$, and the results are shown for $\Omega_{0}dt=1.95$, which is very close to the stability limit for the harmonic oscillator. These results are shown as solid curves along with dotted curves that represent the exact results $\frac{1}{2}mu^{2}$ and $\frac{1}{2}\kappa r^{2}$. As expected from the derivation and analysis of the previous sections, we observe perfect agreement for the harmonic oscillator.

We conclude that the combination of the 2GJ velocity with the GJF trajectory yields very robust statistical results for both kinetic and configurational properties at any time step within the stability range, even for highly nonlinear systems.

III.2 Molecular Dynamics

We wish to further demonstrate the proposed algorithm for use in Molecular Dynamics by showing characteristic simulation results for both solid and liquid states of a simple, one-component system in the $(N,V,T)$ ensemble. We adopt the splined, short range Lennard-Jones potential described in Ref. GJF3 :

[TABLE]

where $r$ is a three dimensional coordinate between any two particles. The parameters are given by

[TABLE]

The potential $E_{p}(r)$ has minimum $-E_{0}$ at $|r|=r_{0}$, and is smoothly splined between the inflection point of the Lennard-Jones potential and zero with continuity through the second derivative at $|r|=r_{s}$, and continuity through the third derivative at $|r|=r_{c}$. We simulate $N=864$ particles of identical mass $m$ in a cubic box of side-length $L$, and with periodic boundary conditions. Initial conditions are chosen to be the Face Centered Cubic (FCC) structure, which is allowed to evolve for a long transient time before acquiring statistical data in order to reach a thermodynamically representative ensemble for a given applied temperature and a reduced time step $\omega_{0}dt$, where $\omega_{0}^{2}mr_{0}^{2}=E_{0}$.

We conduct simulations for two characteristically different temperatures, representing a crystalline solid at a set volume of $V=L^{3}=(8.51442r_{0})^{3}$ for $k_{B}T=0.3E_{0}$, and a liquid at a set volume of $V=L^{3}=(9.378812r_{0})^{3}$ for $k_{B}T=0.7E_{0}$. Each simulation makes statistical averages over a time span of $\omega_{0}\Delta t=2\times 10^{5}$.

The GJF and GJF-2GJ results shown in Figures 3-8 are labeled as follows: For purely configurational measures, shown in Figs. 3 and 6, the results are labeled GJF, since the GJF and GJF-2GJ methods differ only in the velocity. For purely kinetic measures, shown in Figs. 4 and 7, GJF results using the on-site GJF velocity in Eqs. (14) are labeled GJF-v, while GJF results for the half-step 2GJ velocity of Eq. (30) are labeled GJF-u. Finally, for composite measures of kinetic and configurational statistics, shown in Figs. 5 and 8, results using the original GJF method Eqs. (13) and (14) are labeled $GJF$, while results using the GJF trajectory with the 2GJ velocity (e.g., Eqs. (32) and (33)) are labeled GJF-2GJ.

In addition to conducting MD simulations using the GJF and GJF-2GJ methods that we propose in this paper, we have chosen to include a select few of well studied and well-used traditional methods BBK ; Loncharich_2004 ; Pastor_88 that exemplify the inconsistencies between configurational and kinetic properties in discrete time. There are many, more recent methods, which notably includes the key advance shown in Ref. ML , where the proper configurational Boltzmann distribution was first obtained. Given that the GJF method includes this essential feature, as well as other advances, we let the GJF trajectory represent this class of attractive methods. Our main objective here is to show that the GJF-2GJ combination gives very attractive results for MD in both configurational and kinetic sampling.

Apart from the GJF and GJF-2GJ results of this paper, we will thus compare to four other sets of results. These are generated by the well-characterized methods given in Ref. Pastor_88 : Equation (2.20) in Ref. Pastor_88 with the on-site velocity Eq. (3) (BBK-v); Equation (2.20) in Ref. Pastor_88 with the half-step velocity Eq. (6) (BBK-u); Equation (2.21) in Ref. Pastor_88 with the on-site velocity Eq. (3) (PBS-v); and Equation (2.21) in Ref. Pastor_88 with the half-step velocity Eq. (6) (PBS-u).

Figures 3-5 show results for $\alpha=1\,m\omega_{0}$ as a function of the reduced time step $\omega_{0}dt$ for the entire ranges of stability. Figure 3a confirms the attractive configurational sampling properties of the GJF method by displaying only very modest deviation from time-step independent behavior in average potential energies and their fluctuations compared to the traditional BBK and PBS methods. As expected, the liquid state (left panels (b) and (d)) at $k_{B}T=0.7E_{0}$ exhibits the most deviation (about 2-3% near the stability range) since this state is a disordered ensemble with relatively slow diffusion, which makes the statistical sampling more challenging than for a crystalline state (right panels (a) and (c)) at $k_{B}T=0.3E_{0}$.

Figure 4 displays the attractive kinetic sampling properties of the GJF with the 2GJ velocity ($\blacksquare$). The time-step independence of the kinetic measures is striking for both kinetic energy averages and their fluctuations, even in the liquid state at $k_{B}T=0.7E_{0}$. We notice that the PBS-v and, to a lesser degree, the BBK-v results also display very good time-step independence. However, the comparable configurational properties from Figure 3a reveals the internal inconsistencies of those methods. As pointed out in Ref. Pastor_88 , the better use of BBK and PBS is to use BBK-u and PBS-u since these deviations in kinetic properties are somewhat consistent with the deviations in the BBK and PBS configurational measures.

Figure 5 investigates the relationships between the kinetic and configurational statistics with (a) and (b) exploring the difference $\langle E_{k}\rangle-\langle E_{p}\rangle$, and (c) and (d) showing the system pressure ${\cal P}$ calculated from

[TABLE]

where $f_{i}$ is the total force on the particle with coordinate $r_{i}$, and where $T_{k}$ is calculated from Eq. (47). Again, we observe that the GJF-2GJ combination ($\blacksquare$) is superior in its behavior as the time step is increased, with only minor deviations in both statistical energy difference and the pressure. We observe the GJF result ($\bullet$) deviating significantly for the energy differences, due to the known depression of the velocity magnitude, and this translates into the visible depression of the calculated pressure, although the GJF pressure is still better than the ones calculated from the BBK and PBS methods.

Figures 6-8 show the same kinds of results as in Figures 3-5, but for a friction coefficient of $\alpha=10\,m\omega_{0}$. The overall observation from these figures is, as expected, much the same for GJF and GJF-2GJ, since these methods are shown in the earlier sections of this paper to have deviations independent of the friction coefficient. In fact, it is clear that the higher friction allows the system to be in more tight contact with the heat-bath, which results in better statistical averages. Only the fluctuation of the potential energy of the liquid state observed in Figure 6d shows any appreciable deviation (2-3%). The kinetic measures using the 2GJ velocity are as independent of the time step as one can expect. In contrast, we observe that the BBK and PBS results generally show significant deviations as the time step is increased. The expected exception are the PBS-v kinetic results seen in Figure 7, but these very good kinetic properties are not matched by reliable configurational properties, as observed in Figure 6.

IV Discussion

We have presented, analyzed, and tested a new method for modeling Langevin equations in discrete time. The method is based on the trajectory of the GJF thermostat, which has been shown to possess exact configurational sampling properties in linear systems for any choice of time step within the stability limit. We have here devised a definition for a companion half-step velocity, which possesses equally good statistical sampling properties for kinetic measures. The combination of the two parts, the GJF trajectory and the 2GJ velocity, can be utilized in any of the formulated GJF forms, SV, VV, or LF. However, we have formulated a specific LF formulation that directly incorporates the 2GJ velocity into the GJF trajectory calculation, such that this new GJF-2GJ method conveniently displays all the attractive features of both.

Apart from the exact, time-step-independent results obtained for both configurational and kinetic sampling for linear systems, the method has been tested for a nonlinear oscillator subjected to friction and noise, and the results have confirmed near-perfect statistical behavior for as long as the trajectory is stable. Similarly, we have tested the method for MD $(N,V,T)$ ensembles in both crystalline and liquid states for different friction coefficients. Again, the promising features of the method are confirmed for both configurational and kinetic properties as well as for combinations of the two, as exemplified in Figs. 5 and 8 by calculating differences between kinetic and potential energies as well as the pressure as a function of reduced time step with near perfect time-step independence of the results.

Given the simplicity of the method, along with the very familiar Størmer-Verlet form of the expressions, we submit that this is a very convenient method for ensuring the best possible statistical results from Langevin and MD systems. Linear analysis of other methods BBK ; Skeel_2003 ; eb_1980 ; vgb_1982 ; Skeel_2002 ; rc_2003 ; Melchionna_2007 ; Bussi_2007 ; thalmann_2007 ; Mishra_1996 ; ML have demonstrated good second and third order accuracy in the obtained thermodynamic quantities, but we are not aware of any other method that has been demonstrated to give the exact, time-step independent thermodynamic response for both position and velocity. Thus, the method allows for more accurate and more efficient acquisition of reliable statistical data at no additional computational cost compared to any other method we are aware of. The GJF integrator is currently available within the LAMMPS simulation package. Including the 2GJ velocity in simulations seems straightforward.

Appendix A The SV Method Applied to The Deterministic Harmonic Oscillator

The harmonic oscillator is given by Eq. (1) with $f=-\kappa r$, where $\kappa>0$ is the spring constant, and the natural frequency of the oscillator is therefore $\Omega_{0}=\sqrt{\kappa/m}$. The discrete time (complex) solution to, e.g., Eq. (2) yields GJF1

[TABLE]

where the discrete-time frequency $\Omega_{V}>\Omega_{0}$ is given by

[TABLE]

Thus, it is obvious that for as long as $\Omega_{0}dt\leq 2$ (the stability limit) the trajectory $r^{n}$ is a perfect harmonic oscillator with a frequency $\Omega_{V}$ that quadratically approximates $\Omega_{0}$ for small $\Omega_{0}dt$ and becomes $\Omega_{V}\rightarrow\frac{\pi}{2}\Omega_{0}$ for $\Omega_{0}dt\rightarrow 2$.

Following the appendix of Ref. GJF3 , inserting the harmonic solution Eq. (63) into the on-site approximation to the velocity Eq. (3) yields

[TABLE]

Similarly, inserting Eq. (63) into the half-step approximation Eq. (6) to the velocity yields

[TABLE]

A consistent velocity should read $v^{m}=\pm i\Omega_{V}r^{m}$, which means that both the VV velocity and the half-step LF velocity are inconsistent with the trajectory. It is further clear that both velocity definitions yield values that are depressed in magnitude for the harmonic oscillator (convex potential), and we can see that the half-step LF velocity is, as expected, considerably more accurate than the on-site VV velocity. The latter observation is particularly important if one wishes to calculate, e.g., an average kinetic energy (for assessing, e.g., a thermodynamic temperature).

Implied by this short review of the Verlet method applied to the harmonic oscillator is that the discrete-time velocity is not precisely the associated velocity of the simulated trajectory, and that one may benefit from being creative in choosing which discrete-time velocity one applies to certain kinetic measures, such as kinetic energy. Specifically, if one uses SV or VV expressions to approximate a trajectory, it is preferential to apply the half-step LF velocity for kinetic averages, leaving the VV variable $v^{n}$ as a companion variable for evaluating $r^{n}$ in the VV algorithm.