CHARMM Development Project Forums User Discussion & Questions General Chemistry Discussions Calculation of velocities in Langevin Dynamics

Dear advanced users,

Recently I am studying the source code of CHARMM (c36b2) and get a little confused of the calculation of velocities in Langevin Dynamics.

Based on the equations provided in reference [1], one can readily figure out the algorithm of the coordinate propagation.
However, for the calculation of velocities, the code provided for the LEAP method is :

------------------dynlng.src--------------------
gam=timfac*fbeta(i)*delta
gamma(i+natom3)=half*sqrt(one+gam*half)/delta

------------------dynamc.src--------------------
fact=gamma(i+natom3)
vx(i)=(xnew(i)+xold(i))*fact


Thus, there seems to be an additional coefficient, sqrt(one+gam*half), comparing to typical Verlet methods (as in reference [2]).

So I was wondering if there is any theoretical foundation of this coefficient.
Any suggestions would be greatly appreciated.


References:
[1] A. Brünger, C. L. Brooks III, M. Karplus, Stochastic boundary conditions for 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

Langevin dynamics perturbs the velocities, based on a collision frequency (FBETA), which is related to viscosity, and the temperature. It is distinct from Verlet in that energy is not conserved, and it approaches Brownian dynamics behavior in the limit of a large FBETA value (high viscosity). It can be used to regulate T, but a real thermostat is preferred; with Langevin dynamics, the time scale is somewhat arbitrary, and requires careful calibration for comparison to experimental time values.

There's a good treatment of T regulation methods in the book

Computer Simulation of Liquids
Oxford Science Publications
M. P. Allen, D. J. Tildesley
ISBN-13: 978-0198556459 ISBN-10: 0198556454

A brief explanation of Langevin dynamics can be found here

Thank you for the reply!

So, does this mean that the additional coefficient, sqrt(one+gam*half), is somewhat a correction factor for the original Verlet methods when applied in a perturbed system?

For example, accodring to the "Running Molecular Dynamics" section of "dynamc.doc", the Langevin dynamics algorithm presently in CHARMM was intented
to be used primarily with the "Stochastic Boundary Molecular Dynamics" method and consequently has been restricted to an algorithm which is valid only
for the case FBETA*TIMESTEP<<1.0. The algorithm itself reduces to the Verlet algorithm when FBETA is zero and consequently may be used to do regular dynamics.

Thus, to calculate the velocties, instead of using a velocity-Verlet propagator, an original-Verlet-like method is applied. However, since this is a perturbed system, a correction factor is needed and this
factor only works well when FBETA*TIMESTEP<<1.0.

Is this understanding correct? Thank you for the time.

I should note that the LEAPfrog integrator is used, rather the original coordinate Verlet; there are implications for how the velocities are calculated.

There are also velocity verlet integrators available in CHARMM, but I'm not sure they support Langevin dynamics at present.