[gmx-users] trouble diagnosing a simulation that's blowing up (shake not converging --> segmentation fault)

[Nathan K Houtz]

Thanks for your explanations, Dr. Lemkul.

I had already corrected a couple of the things you suggested. Gromacs won't actually let me run with Nose-Hoover and Parrinello-Rahman together (or at least, it gives a warning not to do that and stops). I'd like to run in NVT anyway, so I had set Pcoupl to 'no'. I had also increased the cutoffs and my nstlist is 20.

I have now also changed the tcoupl to v-rescale, but unfortunately that alone didn't help. So I'm not sure exactly what I can inspect to find the source of my error. I know that temperature, pressure, and total energy are all warning signs of bad things if they misbehave, but they were all fairly constant throughout the simulation. (This is for the flexible-model simulation:) The starting energy for 1700 water molecules at 120K was -1.08744e+05, and it finished at -1.06047e+05 with no significant variations on any step. Pressure and temperature were also fine. I attempted to look at the results visually in vmd but I might have done something wrong because the .trr file has some error in it and vmd crashes. But the starting geometry after minimization looks fine: none of the molecules were moved out of their position in the crystal. 

For the constrained molecules, the energy stays about the same (about -1e+5) for the first 8 timesteps and then quickly blows up to +2e15 at step 12 before it obviously fails. I get shake warnings from step 0 though. In vmd, the first 8 steps look reasonable (checking the .pdb files that gromacs outputs when there are warnings) and all the molecules seem to hold their places in the crystal, but then at step 9 suddenly some of the molecules become misshapen with very long or very short bonds. Of course it goes downhill from there, but I can't figure out what's causing the problems since it's clearly happening from the very beginning (according to the shake warnings). 

What other things could I look at to troubleshoot my problem?

Thanks for your help,
Nathan

----- Original Message -----
From: "Justin Lemkul" <jalemkul at vt.edu>
To: gmx-users at gromacs.org
Sent: Thursday, July 9, 2015 7:39:47 AM
Subject: Re: [gmx-users] trouble diagnosing a simulation that's blowing up (shake not converging --> segmentation fault)



On 7/8/15 8:06 PM, Nathan K Houtz wrote:
> Hello,
>
> I deleted the email and can't respond to my last reply directly - sorry! I
> got this response from Mark Abraham:
>
> Hi, Try doing some EM and initial equilibration with no constraints at all,
> perhaps? Mark
>
> I tried commenting out the shake commands, and got a short (5000 step)
> simulation to run just fine without blowing up. Before, I would get shake
> warnings from the first few steps and a segmentation fault around step 13 or
> 14. I would like to be able to simulate with rigid molecules, though. Why
> would the simulation work with flexible molecules but not rigid ones?
>

Flexible water allows weird geometry, which is probably coming up and causing 
your constraint algorithm to fail.  I'd inspect the outcome carefully.  Just 
because it runs doesn't mean it's right.

> Also, in the example .mdp file for tip4p water, there is the (outdated)
> option, 'unconstrained-start', which is now 'continuation'. I got errors when
> trying to make the input .tpr file when I attempted to set that option to
> 'yes'. The warning said it was because I want Gromacs to generate velocities
> to start the simulation, which is incompatible with that command. Is there
> another way I can try to start the simulation unconstrained? Or would you
> suggest another idea to fix my shake warnings?
>

What this setting says is "have the constraints already been solved 
(continuation = yes) or should mdrun constrain the starting configuration itself 
(continuation = no)."

The .mdp file has a number of weird settings.  I would never use Nose-Hoover and 
Parrinello-Rahman when generating velocities; that's likely to be very unstable. 
  See if a more forgiving thermostat and/or barostat resolves the issue, e.g. 
tcoupl = v-rescale and pcoupl = Berendsen.

Also note that the cutoffs are very short (probably because they want to use a 
small box, but note that it *does* affect the physics) and nstlist = 1 is 
totally unnecessary.  It doesn't hurt your physics, but it's a major waste of 
performance.

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================
-- 
Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.