[gmx-developers] Scaling velocities in Velocity Verlet

[David van der Spoel]

Hi,

I'm implementing an algorithm where I change the potential energy Epot 
on time t by an amount DE. In order to conserve total energy I would 
like to change the kinetic energy Ekin by an amount -DE. This is done by 
scaling the velocities at time t by lambda = sqrt(1-DE/Ekin).

The problem is that the change in Ekin comes one step too late. I tried 
to deduce the algorithm from the code (4.5 code base) and deduced the 
following MD algorithm:

v(t) = v(t-dt/2) + a(t) dt/2

v'(t) = lambda v(t)   <======= my addition

v(t+dt/2) = v'(t) + a(t) dt/2

r(t+dt) = r(t) + v(t+dt/2) dt

However, the Ekin that is printed is the uncorrected one even though I 
scale both ekind->ekin and ekind->tcstat[x].ekinf .

Is the algorithm correct as I write it above? Is the Ekin(t) already 
stored at the beginning of the time step?


-- 
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se