CHARMM Development Project Forums User Discussion & Questions Molecular Dynamics Velocity trajectories and reversing dynamics

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, 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?

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.

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
, 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

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?

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.

Brilliant!

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

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



Does that, in fact, follow?

Many thanks!

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.

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.