Reading forces from trr when restarting a simulation

GROMACS version: 2021.2
GROMACS modification: No
Hello Everyone,
Is it possible to let mdrun read the forces, previously stored into trr trajectory file, by restarting the simulation? I tried to find something on that topic in the documentation but I see not clear answer there. Can forces be incorporated into the input files?

Thank you very much in advance for your answer!

Veselin

Hi,
to extract force from trr file you can use gmx traj https://manual.gromacs.org/current/onlinehelp/gmx-traj.html

Structure input file (like gro file) can include information on position and velocity for each particles.

Best regards
Alessandra

You didn’t even read the question properly.

The question is how to organize the simulation in a way that GROMACS’ mdrun can read the forces from the TRR when doing restarts and continue the simulation starting with that forces (without re-computing the forces at the beginning). I need to force the simulator reading the forces from the last frame of the previously generated TRR. Of course, one can change the velocities inside the input to reflect the forces. Another way is to just impose acceleration on certain atoms. The problem with imposing acceleration is that it is implemented upon each step of the integration and one need to perform short simulation with one step of integration, then to remove the acceleration, and to run another simulation using the product of the previous one as in input. What I need is an easy way to create a perturbation at the beginning and measure the relaxation time of that perturbation, without having two simulations. Even answer like “changing the velocities in the input will do it” would be fine.

Please keep a polite tone on the forum. People answering are in no way obliged to answer any questions here. People answering often have to guess what a user might want to do and therefore might give, useful or not useful, suggestions for thing one could do.

mdrun does not read forces from I don’t see the point in mdrun reading forces. Since the potential is a state function, the forces can be recomputed from the coordinates read from the checkpoint file.

Hello Hess,

And thank you for the answer… So if the mdrun does not read the forces, then does it read the velocities from the last frame of a previously generated TRR during the restarts? What if you specify in mdp file “ge-vel = no” but you have the velocity vectors specified in the input TRR file?

Thank you in advance for your answer

Force and velocities are not directly related. mdrun read all data that defines the state point, including the coordinates and velocities, from the checkpoint file and uses those to set the state point to continue the simulation from.

Hello Hess, and thank you for the very clear answer. I had that specific part of the C++ source code gone over and followed the way the mdrun reads the input. Also, even if one supplies very certain the velocity vectors through the input, the v-rescale process will “normalize” them to match the distribution related to the target temperature quite fast. Hence, to achieve my goal, the only option is to compute the acceleration vectors based on the force vectors that are supposed to create the initial perturbation, to define “accelerate group” (one for each heavy protein atom), and to run short simulation to propagate the perturbation. Then one have to remove the accelerated groups in the input topology and continue the MD integration unperturbed to catch the relaxation. I checked that already and the relaxation times observed are in good agreement with the expected values.

I think that solves the problem.

Thank you.

I don’t know how long you need to simulate. If it is short, you can try to decrease the verlet-buffer tolerance to get better energy conservation and run without a thermostat.

Hello Berk and thank you very much for the proposition. So far our goal is to employ the MD integrator of GROMACS for (i) EVB (where upon we need some 50 steps of integration at dt = 2 fs and we don’t need ideal conservation), and (ii) for comparing relaxation time of induced mechanical “spring” effects across single helix (where the conservation is quite important in tracing the “tail” of the relaxation). For the second protocol, we can vary the buffer until the distribution of the potential energy variation in the tail (recorded in the TRR) matches the one predicted by an external method (based on a stochastic protocol that generates samples by means of Metropolis-Hastings criterion). The tests I have done so far show that decreasing the buffer by 3 to 5% brings the distribution derived based on the TRR/EDR analysis reasonably close to the predicted one (verified based on F-test at p=0.85).

Thank you very much again for the proposition!