[gmx-users] Hbond Shake and Gromacs Version Related

[Mark Abraham]

On 21/08/2012 8:45 AM, Sapna Sarupria wrote:
> Hello All,
>
> I have been having rather bizarre experience with Gromacs and was
> wondering if any one can shed some light on what is happening. I am
> running simulations of a hydrate + water system (think of hydrate as
> ice if you are not familiar with them). So basically I have a simple
> solid+liquid system. The system is minimized. The water model I am
> using is TIP4P/2005. Now when I run the simulation using grompp and
> mdrun both from Gromacs 4.5.5, my system crashes and writes out
> step*.pdb files. This I know happens when the system is exploding or
> something however, I do not see an obvious reason for this.
> However, if I compile the tpr file with older Gromacs version 4.0.5
> and then use mdrun from 4.5.5 to run the simulation, my simulation
> proceeds just fine....the way I expect it to.

Minor differences in how things are constructed can have big effects on 
numerical integration. MD is chaotic.

> Can anyone tell me what might be happening here. I have pasted my mdp
> file below.  I think this might be related to the constraint algorithm
> but I have no idea how to fix it. I have tried using Lincs instead of
> Shake and that also gives the same outcome. Runs with Gromacs 4.0.5
> (i.e. grompp using 4.0.5 and mdrun using 4.5.5) and fails with 4.5.5
> (both grompp and mdrun using 4.5.5).
>
> Please let me know if you need any further information and thanks a
> lot for your help. Any insight will be helpful.
>
> #### BEGIN MDP FILE #######
> title               =  s1-Hydrate melting simulation
> dt                  =  0.002                    ; time step
> nsteps              =  5000                    ; number of steps
> nstcomm             =  10                       ; reset c.o.m. motion
> nstxout             =  7500                     ; write coords
> nstvout             =  7500                     ; write velocities
> nstlog              =  25000                    ; print to logfile
> nstenergy           =  250                      ; print energies
> xtc_grps            =  System
> nstxtcout           =  25
> nstlist             =  10                       ; update pairlist
> ns_type             =  grid                     ; pairlist method
> coulombtype         =  PME
> rvdw                =  1.20                     ; cut-off for vdw
> rcoulomb            =  1.20                     ; cut-off for coulomb
> rlist               =  1.20                     ; cut-off for coulomb
> Tcoupl              =  Nose-Hoover
> ref_t               =  200.0
> tc-grps             =  System
> tau_t               =  0.5
> Pcoupl              =  Parrinello-Rahman
> Pcoupltype          =  semiisotropic            ; pressure geometry
> tau_p               =  1.0   1.0                ; p-coupoling time
> compressibility     =  4.5e-5  4.5e-5           ; compressibility
> ref_p               =  100.0   100.0            ; ref pressure
> DispCorr            =  EnerPres                 ; long range correction
> gen_vel             =  yes                      ; generate init. vel
> gen_temp            =  200                      ; init. temp.
> gen_seed            =  98247                    ; random seed
> constraints         =  hbonds                ; constraining bonds with H
> constraint_algorithm = shake
>

As you will see in 
http://www.gromacs.org/Documentation/Terminology/Blowing_Up, the problem 
is usually not with the constraints, but they're the first thing to show 
deviant behaviour. Your use of T and P coupling algorithms is likely 
inappropriate for equilibration given your starting point. Equilibrating 
initially with Berendsen and/or a small time step will likely move you 
to a point where you can switch to the algorithms that produce proper 
ensembles.

Mark