Langevin Dynamics
© Victor Anisimov, 2011
Time evolution of a molecular system is fully described by time-dependent Shrödinger equation
(1)
where Ψ is the wave function of the system, and Ĥ is the Hamilton operator. Eg. 1 describes dynamic changes of coordinate r of electrons and atomic nuclei R in time t as of intricately coupled quantum objects. Based on Born-Oppenheimer approximation, the motion of electrons and nuclei can be separated and the total wave function presented in the form of a product of wave functions of electrons Ψ(r) and nuclei Ψ(R)
(2)
Having done so, for many practical applications it is possible to go further and to consider the motion of nuclei according to classical mechanics described by Newton’s second law
(3)
where m is mass of particle, a(t) its acceleration, and F(t) is force acting on the particle. Eq. 3 describes a closed system with conservative forces, i.e. the total energy of the system is constant, and the work done to move a particle from one point to another is independent of the path. Essential for Eq.3 is a complimentary set of equations of motion
(4)
In theory any classical object abides to the Newton’s low of motion. However, including all interacting particles in a simulation is not always feasible and hardly even necessary. Langevin equation presents an important extension to Newton’s equation for approximate accounting of the effect of environment on molecular system
(5)
where γ is friction coefficient, v(t) is particle velocity, γ is a friction coefficient, and R(t) is Gaussian random force with zero mean and variance of 
, where kB is Boltzmann constant, T0 is thermostat temperature, δ(t-t´) is Dirac delta function. The term R(t) serves as a stochastic force responsible for random collisions of the molecular system with imaginary particles of the environment.
Due to the two extra terms in Eq. 5 vs. Eq. 3 a set of simple relations Eq. 4 no longer holds for Langevin dynamics, thus numerical solution of Eq. 5 becomes seemingly more complicated than that of Eq. 3. However, this limitation of the Langevin equation greatly outweighs the burden of explicitly treating all particles of the environment as it is necessary with Eq. 3.
Brünger-Brooks-Karplus integrator
A popular solution to Langevin equation is given by Brünger, Brooks and Karplus [1], which is known as BBK integrator after the first letters of the authors. However, in their paper the authors gave only the final equation without explaining its derivation. As a result, many readers opening this paper in the search for theoretical insights perhaps felt disappointed when they found such condensed presentation. Therefore, the following is our attempt to reconstruct the logics of the derivation in order to help the reader to grasp the idea of this very important integrator.
To find time evolution of molecular system according to Eq. 5 one may introduce the following approximation, which is implicitly assumed in the BBK integrator [1]
(6)
Strictly speaking, due to extra terms added to force in Langevin dynamics, particle acceleration is no longer a time derivative of particle velocity as it is defined in Newtonian dynamics. The conventional relation between speed and velocity (Eq. 4) will be recovered in the limit of friction coefficient γ approaching zero. Once Eq. 4 is assumed to hold, we can expand particle position in Taylor series about time t in the forward and backward directions according to Verlet algorithm [2] using time increment Δt
(7)
Summing up the corresponding terms in Eq. 7 gives the approach to compute time evolution, which is order 4 accurate in coordinate propagation
(8)
Eq. 8 contains in its right hand particle velocity v(t), which is unknown. It can be obtained from Taylor expansion (Eq. 7) by subtracting second expression from first one
(9)
Using the obtained expression for velocity, Eq. 8 can be rewritten as
(10)
This is the BBK integrator in the form it was given by Brünger, Brooks and Karplus [1]. It is order 4 accurate based on the accuracy of the original expansion (Eq. 8). Substituting velocities from Eq. 9 to Eq. 8 retains the accuracy of Eq. 8 since after multiplication of velocities by Δt2 the the uncertainty coming from Eq. 9 goes down to order 4.
The Verlet form of the BBK integrator does not propagate particle velocity. Therefore velocities at time t can be found from Eq. 9.
Last part in the implementation is discretization of stochastic force, R(t). On discretization the delta function becomes 
. From this it follows that
(11)
Velocity Verlet formulation of the BBK integrator
Browsing a literature on Langevin dynamics the reader may encounter all sorts of different equations called the BBK integrator. In reality these seemingly different equations constitute a class of Langevin dynamics integrators known as the BBK-type integrators. In their root they are all based on the BBK approximation expressed in Eq. 6. The differences start from using different techniques for coordinate and velocity propagation. As shown above, the original BBK integrator is based on the Verlet algorithm. By using leapfrog, position Verlet, and velocity Verlet algorithms one may come to different formulations of the BBK integrator.
Since Verlet algorithm does not give velocities, it may be useful to consider verlocity Verlet integrator which directly computes velocities. Typically velocities computed by integrator are more accurate than those computed from coordinate increments. More accurate velocities allow more accurate temperature determination. Based on Eq. 6 one can propagate velocities by half-time interval
(12)
After rearranging v(t) terms one obtains
(13)
In the next step, coordinate can be propagated
(14)
The coordinate propagation is order 3 accurate because it uses velocities determined for half-step ahead of the coordinate. Substituting Eq. 13 in Eq. 14 explains where order 3 comes from.
Having force computed at time t + Δt makes possible further propagation of velocities
In this equation velocity v(t+Δt) is present in both the left and right hand sides. Resolving the dependence about v(t+Δt) one obtains
(15)
Eqs. 13, 14, and 15 represent the formulation of the BBK integrator based on velocity Verlet algorithm, which is order 3 accurate in coordinate and velocity. Although the velocity Verlet integrator has the advantage of more accurately determined velocities over Verlet formulation, this integrator leads to a bit less accurately determined coordinates.
Limitations of the BBK integrator
Both Verlet and velocity Verlet integrators remain at their stated accuracy as long as Taylor expansion given by Eq. 8 holds. The latter is not fully satisfied in reality, though. In Langevin dynamics a particle can suddenly change its trajectory due to random collision with imaginary particles of the environment. Velocity as time derivative of coordinate will experience discontinuity in the point of collision. Taylor expansion cannot be used under such circumstances since it requires the existence of continuous derivatives. Therefore the actual accuracy of the BBK algorithm is somewhat lower than the estimate obtained from the Taylor expansion. The accuracy of BBK integration approaches its theoretical limit when using smaller friction coefficient γ and integration time step Δt.
Acknowledgements
The author is grateful to Prof. Hiroshi Matsukawa at Dept. of Physics and Mathematics at Aoyama Gakuin University for the stimulating discussion which helped improving this presentation.
References
[1] A. Brünger, C. L. Brooks III, M. Karplus, Stochastic boundary conditions fro molecular dynamics simulations of ST2 water. Chem. Phys. Letters, 1984, 105 (5) 495-500.
[2] L. Verlet, Computer experiments on classical fluids: I. Thermodynamical properties of Lennard-Jones molecules. Phys. Rev., 1967, 159 (1) 98-103
