Langevin thermostat for robust configurational and kinetic sampling

TL;DR

This paper reformulates the GJF Langevin thermostat algorithm as a leap frog scheme, enhancing its robustness and accuracy in configurational and kinetic sampling for molecular dynamics, especially at larger time steps.

Contribution

It introduces a new leap frog formulation of the GJF algorithm that improves kinetic and configurational sampling accuracy and robustness in Langevin dynamics simulations.

Findings

Exact temperature estimation for harmonic oscillators at any stable time step.

Robust configurational sampling in strongly non-linear systems.

Enhanced kinetic measure accuracy with the new formulation.

Abstract

We reformulate the algorithm of Gr{\o}nbech-Jensen and Farago (GJF) for Langevin dynamics simulations at constant temperature. The GJF algorithm has become increasingly popular in molecular dynamics simulations because it provides robust (i.e., insensitive to variations in the time step) and accurate configurational sampling of the phase space with larger time steps than other Langevin thermostats. In the original derivation [Mol. Phys. {\bf 111}, 983 (2013)], the algorithm was formulated as a velocity-Verlet type integrator with an in-site velocity variable. Here, we reformulate it as a leap frog scheme with a half-step velocity variable. In contrast to the original form, the reforumlated one also provides robust and accurate estimations of kinetic measures such as the average kinetic energy. We analytically prove that the newly presented algorithm gives the exact configurational and kinetic temperatures of a harmonic oscillator for any time step smaller than the Verlet stability limit, and use computer simulations to demonstrate the configurational and kinetic robustness of the algorithm in strongly non-linear systems. This property of the new formulation of the GJF thermostat makes it very attractive for implementation in computer simulations.

Equations

$$m\ddot{r}=f-\alpha v+\beta (t),$$

$$⟨\beta (t)⟩$$

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

$$⟨E_p⟩=⟨\frac{k{\left(r^n\right)}^2}{2}⟩=\frac{k_{B}T}{2},$$

$$⟨E_k⟩=⟨\frac{m{\left(v^n\right)}^2}{2}⟩=\frac{k_{B}T}{2}\left[1-\frac{{\left({\Omega }_{0}dt\right)}^2}{4}\right].$$

$$r^{n+1}$$

$$v^{n+1}$$

$$a=\frac{1-\frac{\alpha dt}{2m}}{1+\frac{\alpha dt}{2m}}$$

$$b=\frac{1}{1+\frac{\alpha dt}{2m}}.$$

$${\beta }^{n+1}={\int }_{t_n}^{t_{n+1}}\beta (t^{\prime })d{t}^{\prime },$$

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

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

$$⟨{\left(r^{n+1}\right)}^2⟩=⟨{\left(r^n\right)}^2⟩{\left[1-\frac{b{\left({\Omega }_{0}dt\right)}^2}{4}\right]}^2+b^{2}d{t}^2⟨{\left(v^n\right)}^2⟩+\frac{b^{2}d{t}^2}{4m^2}⟨{\left({\beta }^{n+1}\right)}^2⟩+2⟨r^n{v}^n⟩\left[1-\frac{b{\left({\Omega }_{0}dt\right)}^2}{4}\right]bdt.$$

$${\left(\Omega dt\right)}^2\left(b^2-b+\frac{\alpha dt}{2m}\right)+2⟨r^n{v}^n⟩\left[1-\frac{b{\left({\Omega }_{0}dt\right)}^2}{2}\right]dt=0.$$

$$⟨r^n{v}^n⟩=0.$$

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

$$r^{n+1}$$

$$v^{n+1}$$

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

$$r^{n+1}$$

$$⟨{\left(u^{n+\frac{1}{2}}\right)}^2⟩$$

$$=⟨{\left(v^n\right)}^2⟩$$

$$⟨E_k⟩=⟨\frac{m{\left(u^{n+\frac{1}{2}}\right)}^2}{2}⟩=\frac{k_{B}T}{2},$$

$$w^{n+\frac{1}{2}}=\frac{r^{n+1}-r^n}{\sqrt{b}dt}$$

$$r^{n+1}-r^n=b\left[u^{n+\frac{1}{2}}dt+\frac{dt}{2m}{\beta }^{n+1}\right].$$

$${\left(r^{n+1}-r^n\right)}^2=b^{2}d{t}^2\left[⟨{\left(u^{n+\frac{1}{2}}\right)}^2⟩+\frac{1}{4m^2}⟨{\left({\beta }^{n+1}\right)}^2⟩\right],$$

$${\left(r^{n+1}-r^2\right)}^2=b^{2}d{t}^2\frac{k_{B}T}{m}\left[1+\frac{\alpha dt}{2m}\right]=bd{t}^2\frac{k_{B}T}{m}.$$

$$⟨E_k⟩=⟨\frac{m{\left(w^{n+\frac{1}{2}}\right)}^2}{2}⟩=\frac{k_{B}T}{2}.$$

$$r^{n+1}=2b{r}^n-a{r}^{n-1}+\frac{bd{t}^2}{m}{f}^n+\frac{bdt}{2m}\left({\beta }^n+{\beta }^{n+1}\right),$$

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

$$r^{n+1}$$

$$T_c=\frac{1}{k_B}\frac{⟨{\left(\frac{\partial U}{\partial r}\right)}^2⟩}{⟨\frac{{\partial }^{2}U}{\partial r^2}⟩}.$$

$$T_k=\frac{2}{k_B}⟨E_k⟩,$$

Full text

I Introduction

One of the prominent approaches for conducting molecular simulations in the canonical $(N,V,T)$ ensemble is based on the idea that the statistical ensemble can be sampled by considering the dynamics of each particle in the system to be governed by Langevin equation langevin:1908

[TABLE]

where $r$ and $v=\dot{r}$ denote, respectively, the coordinate and velocity of the particle. Langevin’s equation is essentially Newton’s second law describing the motion of a particle of mass $m$ under the action of (i) a deterministic force, $f$, and two additional forces representing the interactions with a heat bath - (ii) a friction force, $-\alpha v$, where $\alpha>0$ is the friction coefficient, and (iii) Gaussian white noise with zero mean and delta-function auto-correlation risken:book :

[TABLE]

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

In computer simulations the time is discretized in intervals of $dt$, and a Langevin “thermostat” algorithm is used for discrete-time integration of Langevin’s equation of motion, yielding a sequence of coordinates $r^{n}=r(t_{n})$ and velocities $v^{n}=v(t_{n})$, where $t_{n}=ndt$. A major problem of thermostat algorithms is the discretization errors that they introduce pbs88 . These cause computed averages of thermodynamic quantities of interest to vary with the time step $dt$ - an alarming feature that raises concerns about the reliability of the simulation results. A simple test-case for the robustness of an algorithm is the one-dimensional harmonic oscillator where $f=-kr$. The average potential energy satisfies $\langle E_{p}\rangle=\langle kr^{2}/2\rangle=k_{B}T/2$, but most popular and widely-implemented algorithms, e.g., BBK bbk84 , Schneider-Stoll ss78 , and van Gunsteren-Berendsen vGB82 thermostats, exhibit systematic deviations from this result. Depending on the method of choice, the integration error of the potential energy may scale as ${\cal O}(dt)$ or ${\cal O}(dt^{2})$ wang03 . It was only several years ago that two new algorithms were introduced, by Leimkuhler and Matthews (LM) lm12 and by Grønbech-Jensen and Farago (GJF) gjf1 , that reproduce the exact harmonic oscillator potential energy for any time step within the stability limit $dt<2/\Omega_{0}$, where $\Omega_{0}=\sqrt{k/m}$ is the frequency of the oscillator.

When conducting a molecular simulations study, one is often interested in measuring the temperature of the simulated system in order to compare it to the target thermodynamic temperature. The most straightforward quantity to calculate for this purpose is the average kinetic energy per degree of freedom, $\langle E_{k}\rangle=\langle mv^{2}/2\rangle=k_{B}T/2$. Unfortunately, the discrete-time variables $r^{n}$ and $v^{n}$ are only approximations of their continuous-time counterparts. In contrast to the latters, the formers are not exactly conjugated to each other, which causes the “kinetic” and “configurational” measures of the temperature to be different. This feature is nicely captured by the harmonic oscillator test-case. As mentioned above, the LM and GJF algorithms yield the correct configurational temperature,

[TABLE]

but the kinetic temperature computed by these thermostats exhibits a discretization error and reads

[TABLE]

There exist other thermostats that reproduce the kinetic energy without discretization errors, but no existing algorithm has simultaneously both the correct kinetic and potential energy of the harmonic oscillator. Since the aim of computer simulation studies of molecular systems at equilibrium is phase space sampling, the velocity variable is essentially an auxiliary field and one should favor the use of algorithms like the GJF thermostat, which have been demonstrated to provide robust configurational sampling not only for the harmonic oscillator but also for non-linear molecular systems gjf2 ; gjf3 . Nevertheless, the kinetic energy constitutes a useful and a simple measure for the temperature of the system and, therefore, a question arises on whether it is possible to devise a thermostat featuring both correct potential and kinetic energies of the harmonic oscillator. Here, we show that the GJF algorithm can be reformulated with a different velocity variable which, in contrast to the one in the original formulation, exhibits no discretization errors. We use simulations of a simple toy model to demonstrate the robustness of the newly-defined velocity also in non-linear systems.

II Half-step velocity

Our starting point is the GJF algorithm, which in the velocity-Verlet formulation reads gjf1

[TABLE]

where $f^{n}=f(r^{n})$, and the damping coefficients of the algorithm are given by

[TABLE]

and

[TABLE]

The discrete-time noise,

[TABLE]

is a random Gaussian number satisfying

[TABLE]

where $\delta_{m,n}$ is Kronecker delta.

We now invoke another property of the canonical ensemble, which is the fact that $r$ and $v$ are statistically independent degrees of freedom, namely $\langle rv\rangle=0$. Let us demonstrate that the discrete-time variables in the GJF algorithm satisfy this relation. To show this, we begin by squaring Eq. (6) and taking statistical averages of all terms. Keeping also in mind that $f^{n}=-kr^{n}$, this yields

[TABLE]

Using Eqs. (4), (5), and (12), and the fact that $\left\langle\left(r^{n+1}\right)^{2}\right\rangle=\left\langle\left(r^{n}\right)^{2}\right\rangle$ in Eq. (13), yields the equality

[TABLE]

The first term on the r.h.s. of Eq. (14) vanishes because $b^{2}-b+\frac{\alpha dt}{2m}=0$, and we immediately conclude that

[TABLE]

Let us look at the “half-step” velocity variable, $u^{n+1/2}$, defined by rewriting the frictionless ($\alpha=0$) velocity-Verlet algorithm verlet67 in the following form:

[TABLE]

With the definition of $u^{n+1/2}$ by Eq. (16), the GJF equations (6)-(7) in the velocity-Verlet form, can be converted into the following leap-frog form

[TABLE]

These equations constitute a new formulation of the GJF thermostat, to be henceforth referred to as the GJF-F algorithm. The half-step velocity variable $u^{n+1/2}$ satisfies

[TABLE]

and using Eqs. (4), (5), and (15) we readily find the kinetic energy associated with $u^{n+1/2}$

[TABLE]

which is exact for any time step $dt$ (within the stability limit).

With the above derivation, it is easy to define another half-step velocity

[TABLE]

with similar properties. This velocity variable was independently identified recently by Grønbech-Jensen and Grønbech-Jensen (2GJ) 2gj . In order to prove that $w^{n+1/2}$ is a robust velocity variable, we rewrite Eq. (20) in a slightly different form

[TABLE]

Squaring both sides of Eq. (24) and taking averages, we arrive at

[TABLE]

and by using Eqs. (22) and (12) we find that

[TABLE]

From the last result we immediately conclude that for any $dt$

[TABLE]

A leap-frog scheme involving $w^{n+1/2}$ can be derived by complementing Eq. (23) with the Strømer-Verlet form of the GJF algorithm (see Eq. (11) in ref. gjf2 )

[TABLE]

which, together with the relationship $a+1=2b$, leads to the following set of equations

[TABLE]

This scheme was presented in ref. 2gj and was termed the GJF-2GJ algorithm.

III Simulations of a non-linear model

To test the robustness of the new velocity variables beyond the harmonic oscillator test case, we consider the non-linear model presented in ref. gjf3 of a particle moving in a one-dimensional potential $U(r)=kr^{2}/2-\cos(r-\xi)$, with $k=1/40$, and $\xi=3/4\pi$ [see fig. 1(a)]. In ref. gjf3 , this model provided a demonstration for the superiority of the GJF algorithm over classical popular thermostats (BBK, SS, vGB) in configurational sampling. This was done by computing the configurational temperature defined by

[TABLE]

We now wish to also measure the kinetic temperature

[TABLE]

and explore the dependence of this quantity of the simulation time step $dt$. Our simulation results for the dependence of the configurational and kinetic temperatures on $dt$ are summarized in figs. 1(b)-(d). In the simulations we set $m=1$ and $T=1$ (the thermodynamic temperature), and use three different values of $\alpha$: $\alpha=0.1$ [fig. 1(b)], $\alpha=1$ [1(c)], $\alpha=10$ [1(d)]. Based on the simulation results for this model in ref. gjf3 , we restrict the simulations to the range $0<dt\leq 1$ at which the GJF algorithm exhibits accurate configurational sampling. For $dt>1$, discrepancies between $T_{c}$ and $T$ in the low friction simulations become noticeable (relative error $>6\%$). We measure the configurational temperature, $T_{c}$ [denoted by black circles in figs. 1(b)-(d)], and three kinetic temperatures $T_{k}^{v}$ (red squared), $T_{k}^{u}$ (green diamonds), $T_{k}^{w}$ (blue triangles) corresponding, respectively, to the in-site velocity $v^{n}$ defined in the original velocity-Velret GJF algorithm [Eq. (7)], and the two half-step velocities $u^{n+1/2}$ [Eq. (19)] and $w^{n+1/2}$ [Eq. (29)] introduced in the GJF-F and GJF-2GJ leap-frog formulations of the GJF algorithm. The simulation results in figs. 1(b)-(d) clearly demonstrate the difference between the in-site and half-step velocity variables. While the kinetic temperature associated with the former tends to decrease with $dt$, the kinetic energy of the latters remains extremely close to the average thermodynamic kinetic energy ($0.99<T_{k}/T<1.015$) for any time step within the range simulated herein. These results corroborate the intuition from the harmonic oscillator analysis that the half-step discrete-time velocity variables $u^{n+1/2}$ and $w^{n+1/2}$ are robust to time step variations also in non-linear systems. This conclusion agrees with the recent findings reported in ref. 2gj , where the robustness of the half-step velocity $w^{n+1/2}$ was demonstrated in simulations of three-dimensional Lennard-Jones systems.

IV Summary

We have introduced the GJF-F algorithm, Eqs. (19)-(20), which is a new formulation of the GJF algorithm for Langevin dynamics simulations. In this formulation, the GJF thermostat is represented as a leap-frog scheme with half-step velocity $u^{n+1/2}$. In contrast to the in-site velocity variable $v^{n}$ appearing in the original GJF algorithm, the half-step velocity in the new GJF-F algorithm exhibits robustness to time step variations when applied to the harmonic oscillator problem. Computer simulations demonstrate that this feature of the half-step velocity is also observed in strongly non-linear systems. Thus, the newly-presented method allows for both accurate configurational and kinetic sampling of canonical ensembles. This makes the method very attractive for implementation in computer simulations. On the one hand, it generates the same trajectories, $\{r^{n}\}$, like the GJF algorithm and thus provides high quality configurational sampling with larger time steps compared to other popular Langevin thermostats. On the other hand, it also provides robust kinetic sampling, which offers a convenient way to assessing the temperature of the simulated system via the average kinetic energy.

Acknowledgments: I thank Niels Grønbech-Jensen for stimulating discussions, especially related to the differences between on-site and half-step velocities. This work was supported by the Israel Science Foundation (ISF) through Grant No. 991/17.