On Thu, Mar 28, 2002, Margaret Cheung wrote:
>
>
> If I understand the leap-frog correctly, the expression for the new
> velocities, based on the old velocities:
>
> v(t+dt/2) = v(t-dt/2) + dt* f(t)/m
>
> The update of the new coordiates is:
>
> r(t+dt)=r(t)+dt*v(t+dt/2)
>
> Which implied that kinetic energy (obtained at time t+dt/2) and potential
> energy (obtained at time t+dt) can not be obtained
> at the same time.
>
For "strict" Verlet, this is true. There are various approaches to getting
the KE at the same time as you have the PE (from the force() routine).
Amber 4,5,6 used the following ansatz:
KE(t) approx= [ KE(t-dt/2) + KE(t+dt/2) ]/2
i.e. averaging the kinetic energy at the two points where we know the
velocities. These values are called EOLD3 and ENEW in the code.
Amber 7 uses a slightly better algorithm: do a "half-step" update to
propagate the velocities from t-dt/2 to t; then use the V(t) values to
compute the KE.
Note that the method of estimating the KE has no effect on the trajectory;
this is why I say that how to get the KE goes "beyond" the Verlet
integrator.
[Note that all of the above becomes more complicated when SHAKE is turned on.]
> However, in the code (runmd.f)
>
> velocity is obtained by
> V(I3) = (V(I3) + F(I3)*wfac)
This is just your first equation, above.
>
> and the coordinates are later obtained by
> X(I3) = X(I3)+V(I3)*DTX
This is your second equation, above.
...hope this helps....dac
--
==================================================================
David A. Case | e-mail: case_at_scripps.edu
Dept. of Molecular Biology, TPC15 | fax: +1-858-784-8896
The Scripps Research Institute | phone: +1-858-784-9768
10550 N. Torrey Pines Rd. | home page:
La Jolla CA 92037 USA | http://www.scripps.edu/case
==================================================================
Received on Thu Mar 28 2002 - 17:52:38 PST
