On 7/10/15 7:21 PM, Justin Lemkul wrote:
On 7/10/15 7:15 PM, Nathan K Houtz wrote:
Actually, I think I found the problem. When i looked at the pdb files in vmd
the first time I missed it, but a colleague had a hunch and I found that he
was right! I think that shake is constraining atoms to the wrong molecules.
Here's a screenshot I took where you can see what I
mean:http://imgur.com/1bgThDv
Sorry, this doesn't make sense to me. Constraints are applied per-moleculetype;
they cannot be applied between molecules. I can't tell anything from this
image.
I'm still not sure how to fix it though. The molecule numbering in my .gro
file is correct, and there are only 4 atoms per molecule (or 3 plus the
dummy). How can I tell shake to look at each molecule individually? Or have I
made a mistake somewhere? Here are my most recent files:
.gro file: (only part of it -- there is a character limit. this is the first
and last ten molecules with the box dimensions at the bottom)
http://textuploader.com/e898
topology file:
http://textuploader.com/e89g
minimization:
http://textuploader.com/e8rc
.mdp file:
http://textuploader.com/e896
One final thing after I took another look on a hunch - constraining all angles
is very unstable, and doing it with SHAKE is (I think) even worse. You've got
[constraints] manually defined in the topology to keep the geometry rigid; you
should set "constraints = none" to just use those constraints specified in the
topology.
-Justin
log file from a failed simulation:
http://textuploader.com/e80y
Thanks for your help! If you do still want the full .gro file to run the
simulation, I can look for a different file host.
I really would like a nice tarball with everything so that I can play with it a
bit. That's the only way I'm going to have any chance of figuring it out.
I would also dispute your previous assertion that the temperature is fine - it
goes out of control immediately :)
-Justin
Nathan
----- Original Message -----
From: "Justin Lemkul" <jalem...@vt.edu>
To: gmx-us...@gromacs.org
Sent: Friday, July 10, 2015 5:06:34 PM
Subject: Re: [gmx-users] trouble diagnosing a simulation that's blowing up
(shake not converging --> segmentation fault)
On 7/9/15 9:23 PM, Nathan K Houtz wrote:
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?
Not sure, but at least the failure happens fast. Upload all your input files
somewhere and I'll take a few minutes to see if I can spot anything.
-Justin
Thanks for your help,
Nathan
----- Original Message -----
From: "Justin Lemkul" <jalem...@vt.edu>
To: gmx-us...@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
jalem...@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-requ...@gromacs.org.
