[gmx-users] rerun forces

[Mark Abraham]

Please include the whole .mdp file. The simplest explanations for what you
describe are that you are attempting to restart with an inappropriate
combination of inputs, or that your .mdp file specifies generation of new
velocities upon restart. There's a few other comments below.

Mark

> Hi all,
> I have a little problem with the rerun of the trajectories to get out
> forces acting on atoms.
> To calculate the forces I followed the instructions I read in the
> mailing list.

If you have a reference, please give a link to it - maybe they were wrong :-)

> I changed the mdp file, as below, to get the forces every frame, without
> re-calculating the neighbor list....
>
> -----------------
> mdp ORIGINAL
> ....
> ;  OUTPUT CONTROL
> nstcheckpoint       =  1000     ; checkpointing helps you continue after
> crashes

If this option exists, it isn't documented!

> nstxtcout           =  1000              ; output frequency and
> precisionfor xtc file

This is only the output frequency, not precision

> nstxout             =  1000           ; write coords
> nstvout             =  1000
> nstfout             =  1000
> nstlog              =  1000     ; print to logfile
> nstenergy           =  1000            ; print energies
> ; NEIGHBOR SEARCHING
> nstlist             =  15                       ; update pairlist
> ns_type             =  grid                     ; pairlist method
> pbc                 =  xyz
> rlist               =  1.2                      ; cut-off for ns
>
> -----------------------
> mdp RERUN
> ....
> ;  OUTPUT CONTROL
> nstcheckpoint       =  1000000000     ; checkpointing helps you continue
> after crashes
> nstxtcout           =  1000000000              ; output frequency and
> precisionfor xtc file
> nstxout             =  1000000000           ; write coords
> nstvout             =  1000000000
> nstfout             =  1
> nstlog              =  1000000000     ; print to logfile
> nstenergy           =  1000000000     ; print energies
> ;  NEIGHBOR SEARCHING
> nstlist             =  0                       ; update pairlist
> ns_type             =  grid                     ; pairlist method
> pbc                 =  xyz
> rlist               =  1.2                      ; cut-off for ns
>
> I obtain forces of the order of 10e+30 Kj mol-1 nm-1 (are these the
> correct units, to divide by f~139 to get forces in (e/nm)^2 ?), that
> sounds quite strange.

There's a section on units in the start of the gromacs manual. I suggest
you have a read of it. I can't think where forces this large are coming
from without more information.

> In fact I run a new short (100ps) simulation in which forces are printed
> out, and I obtain much smaller forces, between 10e+0 and 10e+5.

They sound more correct.

> Then I tried putting _nstlist = 1_, updating the forces each step during
> the rerun (yes?) of the trajectory, and I obtain values of the same
> order of magnitude, but with different values of the forces.
> I would expect to find the same forces at least for the frame number 0.
> What's wrong?  I need to do some imaging of the trajectory, or to leave
> unchanged also the nslist flag? In my case, the results do not vary by
> changing neither the flag nstfout from 1 to 1000, neither the nstlist
> from 1 to 1000.

Restarts (and reruns of restarts) require that you supply the previous
.trr and .edr files to grompp so that all the coordinates, box parameters,
etc. are correct. You will also want unconstrained_start = yes and gen_vel
= no to avoid performing geometric constraints before the first step
(since they were applied after the last step of the previous simulation)
and to use the velocities in the .tpr file rather than invent new ones.

Mark