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
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Only plain text messages are allowed!
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
