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

[Justin Lemkul]

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
>
> 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" <jalemkul at vt.edu>
> To: gmx-users at 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" <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

==================================================