Previous Thread
Next Thread
Print Thread
Joined: Feb 2010
Posts: 29
Ascaris Offline OP
Forum Member
OP Offline
Forum Member
Joined: Feb 2010
Posts: 29
I've looked in the documentation, and I've looked in the forums (and charmmtutorial.org, for that matter) but I haven't been able to find the answer so I thought I'd ask here.

When starting dynamics using, for example,
Code:
DYNA LEAP VERLET STARt -
 NSAVC 1 NSAVV 1 ... 
from what time step are the first velocities that are written to the VEL trajectory?

Based on my understanding of the leapfrog integrator, it is likely v(t0 + 0.5h), where h is the timestep ... in other words, a half (time) step ahead of the starting position coordinates and a half step behind the first set of position coordinates written to the DCD. Is that correct?

If that is not correct, please tell me what time step the first velocities written to the VEL trajectory are from, when using the LEAP VERLET options as specified above. (And read no further.)

If indeed the first set of velocities written in the .vel trajectory are the velocities at time t=t0+0.5h, then it would seem that to reverse dynamics, using the negative of those velocities (i.e., -v(t0+0.5h)) with the first set of coordinates from the dcd (i.e., x(t0+h)) and not the starting coordinates (i.e., x(t0)) would be the way to extend the trajectory backwards in time. The first step in the backwards direction would bring you back to x(t0). However, it is necessary to use x(t0+h) as the starting point for integration backwards in time because the earliest velocities that you know are from t=t0+0.5h. Subsequent steps would bring you into the t < t0 range. Does that follow?


Intestinal nematode
Bronx, NY
Joined: Dec 2005
Posts: 1,535
Forum Member
Offline
Forum Member
Joined: Dec 2005
Posts: 1,535
Sounds about right, but I don't think it will get you very far; I'd gamble that the chaotic nature of MD simulations (ie. the amplification of very small errors over time) will quickly make the reverse trajectory diverge from the forward one.

Joined: Feb 2010
Posts: 29
Ascaris Offline OP
Forum Member
OP Offline
Forum Member
Joined: Feb 2010
Posts: 29
Assuming for a second that one might have a reason to run a trajectory forwards and backwards in time from a starting position (and not just doubling back over a trajectory), I'd like to clarify things a bit further.

Returning to the scenario I outlined above, but going into slightly more detail, if one's start command corresponded to
Code:
DYNA LEAP VERLET STARt -
 NSAVC 1 NSAVV 1 -
 IUNCRD 30 IUNVEL 40 -
 IASORS 1 IASVEL 1
, would the initial velocities taken from the gaussian distribution (specified by IASVEL .gt. 0) correspond to v(t0) or to v(t=t0+0.5h)? I suppose that is a moot point because, regardless, we are saying that the first velocities written to the VEL trajectory (unit 40) are v(t=t0+0.5h). Am I understanding things correctly?

On the other hand, if one uses velocity verlet, such as
Code:
DYNA VVER STAR -  
 NSAVC 1 NSAVV 1 -
 IUNCRD 30 IUNVEL 40 -
 IASORS 1 IASVEL 1

then the velocities are assigned at t0 and the algorithm is propagated such that the first velocities written to unit 40 are v(t0+h). Similarly, the first coordinates written to unit 30 are x(t0+h). In this case, the initial velocities, assigned from the Gaussian distribution, v(t0), are not recorded. Therefore, to propagate backwards, one would have to use the negative of the first velocities from unit 40 (i.e., –v(t0+h)) along with x(t0+h) to propagate backwards in time. Does that also follow?


Intestinal nematode
Bronx, NY
Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
The velocities written for LEAP VERLET are on-step, computed as the linear interpolation of the half-step velocities. This is not typical for leapfrog integrators, but is a CHARMM feature added to facilitate the calculation of the virial for constant P using the extended Hamiltonian method.


Rick Venable
computational chemist

Joined: Feb 2010
Posts: 29
Ascaris Offline OP
Forum Member
OP Offline
Forum Member
Joined: Feb 2010
Posts: 29
Brilliant!

I had a feeling you would know. The only question was whether you or lennart would respond first. smile

And, just to confirm, my previous reasoning about vverl, therefore, applies here (because the leap velocities are also on-step). Namely,

Originally Posted By: Ascaris
In this case, the initial velocities assigned from the Gaussian distribution are not recorded. Therefore, to propagate backwards, one would have to use the negative of the first velocities from unit 40 (i.e., –v(t0+h)) along with x(t0+h) to propagate backwards in time. Does that also follow?


Does that, in fact, follow?

Many thanks!


Intestinal nematode
Bronx, NY
Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
I suspected but didn't know for sure about the on-step velocities, but I knew the right person to ask.

Using the negative velocities does seem logical.


Rick Venable
computational chemist

Joined: Dec 2005
Posts: 1,535
Forum Member
Offline
Forum Member
Joined: Dec 2005
Posts: 1,535
I'm still extrememly curious what the purpose of this excercise is; I would think the "reverse" trajectory will in practice be nothing but another "forward" trajectory with a different (and equally valid) set of initial velocities.


Moderated by  BRBrooks, lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.5
(Release build 20201027)
Responsive Width:

PHP: 7.3.31-1~deb10u1 Page Time: 0.012s Queries: 28 (0.008s) Memory: 0.7624 MB (Peak: 0.8287 MB) Data Comp: Off Server Time: 2022-09-27 16:45:27 UTC
Valid HTML 5 and Valid CSS