[gmx-developers] Scaling velocities in Velocity Verlet

[Berk Hess]

Hi,

With which integrator did you try this?
Leap-frog doesn't have v at whole time steps.

Cheers,

Berk

On 12/3/11 9:22 , David van der Spoel wrote:
> 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?
>
>