- Contemporary messages sorted: [ by date ] [ by thread ] [ by subject ] [ by author ] [ by messages with attachments ]

From: David Case <case_at_scripps.edu>

Date: Thu 28 Mar 2002 17:52:38 -0800

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

Date: Thu 28 Mar 2002 17:52:38 -0800

On Thu, Mar 28, 2002, Margaret Cheung wrote:

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.]

This is just your first equation, above.

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

Custom Search