Trajectory continuation with old gromacs

GROMACS version: 4.5.5
GROMACS modification: No

Hi everybody,

I would like to continue a trajectory of gromacs from the last frame, therefore I made the state.cpt file of it. Furthermore, I made the same trajectory, this time with the full length I desired in order to understand the correct way to lengthen the first mentioned: in that case, indeed, the prosecution will be the same of the one I see in the full one. Some numbers: 1st traj stops at frame 14000, whereas the full one goes to 40000.

For several reasons, I have to extend the first trajectory via gromacs 4.5.5 (both trajectories have been generated with it, of course) and, looking at the gromacs website and also at various forums, I am doing the steps that will be listed in the following. I please ask you what is wrong, since the first frame of the new trajectory generated (i.e. the prosectuon of the 1st, to be compared to the full length one) is not the same of the 14001 of the full trajectory I run, neither the last one of the 1st trajectory (which could have been a possibility, indeed). Here are the instructions I follow:

  1. grompp -c template_system_structure.gro -f mdpfile.mdp -t state.cpt -p system.top -n system.ndx -o system.tpr

  2. mdrun -s system_new.tpr -t precedent_traj.trr -append

Indeed, I’ve seen that if I use the .cpt file to generate the system_new.tpr, I don’t need to state the .cpt file again in the second instruction. Moreover, I stress that I’m doing these two commands because others do not work in gromacs 4.5.5, so please don’t tell “try doing this other way ….” because I really, really have tried the whole of them.

Of course, I changed the .mdp file taking care in the prosecution, taking care of the subsequent parameters:

init_step = 10400000 (and also tried 10401000. I have dt = 0.001 ps, so that’s why the 3 additional zeros)

gen_vel = no

continuation = yes

And, just to be sure, I tried tinit = 0, 10400 and 10401 (which are correct values and do not have to be of the same order of init_step, since they are referred to the frame, not the step. Therefore 10400000 would be an enormous error). The rest of the mdp is identic to the first I used for the run of the 1st trajectory.

I tried all the combinations of init_step and tinit. I am sure the state.cpt is the last frame of the 1st traj (I checked via check -f state.cpt). Anyone knowing how to exactly set this continuation of the run? Thanks in advance!

Here is the link of one of the .mdp I tried for the prosecution: mdout.mdp - Google Drive

Best Regards,
Jacopo

I would need some more information to help you. Do you mean that the two trajectories (1st to frame 14000 and the second to 40000) were generated from different TPR files originally? In that case, were the starting velocities the same? I suspect your problem is not that extending of the first trajectory is not working properly, but rather that the trajectories were different already before that.

You say that the first frame of the extended trajectory is not the same as in the full trajectory (of corresponding time stamp). Is the last frame of the first trajectory (before extending) the same as in the full trajectory, of the same time. If those two frames were identical the problem would indeed be from when you extend the simulations (as you write). Otherwise it can be of multiple reasons, i.e., different starting temperatures or slight deviations over the large number of steps. Small deviations are expected, especially if running in single precision for very many steps. Running in an NVE ensemble would reduce some of the sources of deviations, but is not necessarily of interest for most purposes.

Dear MagnusL,

I really I appreciate you dedicated some time to me, thanks. Before opening the thread about my issure, I made sure the two original TPR files were the same, indeed the trajectories are identical up to the frame 10400 at which I stopped the shortest. Apart from the changes in the new .mdp file I explicited in my previous message and the instructions I listed to update the .tpr for the prosecution, I did not change at all the original .mdp. Hope the error lies in the fact that I should change something more in the .mdp file I shared to you to ensure a corret propagation, as if it is not I really can’t figure out what might be the problem.

Best Regards,
Jacopo

I see. Then it seems to be related to the continuation, indeed.

You said that the first trajectory stops at frame 14000, from the mdp file and the rest of your information, I presume that you mean frame 10400, right? With the settings in the mdp file it will only continue to frame 20000. Did you use the same number of steps for the full trajectory? If so, how did you get 40000 frames?

However, that still does not explain the difference in frames after restarting. Your init_step setting looks fine (assuming that you are starting from step 1040000, i.e., frame 10400), so I don’t have any good explanation right now.

Dear MagnusL,

Good morning and kindly thanks again. I apologyze for the typo: I always intended 10400. I see your point about the .mdp file going at maximum to step 20000. Indeed, I took a chunk of the trajectory that is being lengthened (from 10401 to 15000) to notice that the coincidence with the long traj is not respected, already. My intention was to verify this mismatch before going up to the full length again. I confirm, the init_step is fine.

Good news is that I finally managed to solve the problem. Even if that could be redundant, I ried to put the flag: “-cpi state.cpt” in the 2. instruction I reported, although I expected that the “-t state.cpt” flag in instruction 1. would have made it unnecessary. I’m looking at the lengthening right now, and finally see that for the first 1200 frames of continuation it is coincident to the corresponding part of the full trajectory. Therefore, the .mdp file generated from the continuation - which owns the differences I reported in my first message, was ok. Instead, something was still missing in command 2. Hope you found our discussion “trajectory-productive” :)

Thanks a lot for your efforts!

Best Regards,
Jacopo D’Ascenzi

Good to hear that it works now. I’m also a bit surprised that you would need to specify the checkpoint as input both to grompp and mdrun. But I see that it says so in the grompp help text as well: “Actually preserving the ensemble (if possible) still requires passing the checkpoint file to gmx mdrun -cpi.”

I see, I am sorry and I wish I had looked at the help text more carefully than I did, since I missed the line “You then supply the old checkpoint file directly to mdrun with -cpi” it returns and found it counterintuitive, given the “-t” flag in instruction 1