[gmx-users] Constraint causing system to explode

[Warren Gallin]

OK, I've done as Mark suggested, replaced the type 2 constraint between the amino terminal nitrogen atom and the C-terminal carbon atom with a type 6 bond, r0=0.8, r1=0.8, r2=1, which acts as a distance restraint.

This has run clean without crashing for 1.6 ns of simulation on a single processor, so my explosion is not originating from the linkage between the ends of the peptide per se.

In this case the distance between the the two atoms increases from the initial 0.8 nm to what appears to be an equilibrium oscillating around 1.5 nm, which is what I would expect - the most probable distance between the ends of the peptide is much larger than 0.8 nm, so it's essentially pushing the spring.

Which leaves me with the next question:  What is it about the type 2 constraint that is making things blow up?

Intuitively it seems possible that after an unusual step of simulation a water molecule gets between the two ends or atoms attached to the two constrained atom rotate into line between the two constrained atoms, and when the constraint pulls the ends back to 0.8 nm  it creates a ridiculous force.

Is that a valid worry?

But why would that be true only for a type 2 constraint between the ends of the peptide, and not all the other bond constraints in the rest of the peptide?  Should I be considering the use of a type 1 constraint (with exclusions), and if so, is that formulation legitimate for determining the free energy?

Any suggestions on how to proceed?  All ideas welcome.

Warren Gallin

On 2010-05-31, at 10:09 PM, Mark Abraham wrote:

> ----- Original Message -----
> From: Warren Gallin <wgallin at ualberta.ca>
> Date: Tuesday, June 1, 2010 13:38
> Subject: Re: [gmx-users] Constraint causing system to explode
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> 
>> Chris,
>> 
>> 	As you suggested, I ran it on one node, and it still blew up, 
>> so it appears to be some oddity with the constraint that I am 
>> imposing.
> 
> Does it happen with a distance restraint, instead? (Is that a better physical model?) Is there any possible correlation with one atom crossing a PBC boundary, such that a buggy implementation would then be horribly broken on the next step?
> 
> Mark
> 
>> 	When I look at the dgdl.xvg file from the run, there appear to 
>> be a number of transient spikes that are a couple orders of 
>> magnitude larger than the typical variation during the run, and 
>> I am getting errors on an angle changing by more than 30 degrees 
>> as a LINCS warning during those short spikes.  The final 
>> system blow-up starts with the same situtation, which 
>> degenerates into a full blow system explosion, with bond lengths 
>> going to crzy high values.
>> 
>> 	So my question boils down to this:  What is it about a 
>> type 2 constraint that will generate these ridiculously large 
>> forces on relatively rare occasions?  Or perhaps, what 
>> error am I making in setting up the simulation that would cause 
>> only the constrained system to manifest this behavior?
>> 
>> Warren Gallin
>> 
>> 
>> On 2010-05-27, at 11:22 AM, chris.neale at utoronto.ca wrote:
>> 
>>> Dear Warren:
>>> 
>>> I don't have your answer, but I'll point out that when you 
>> ask: "Is it possible that this is a problem that arises because 
>> of domain decomposition over multiple nodes" that you are 
>> probably the person in the best position to address this. 300ps 
>> should not take too long to simulate so why not try it on a 
>> single node and also on multiple nodes with mdrun -pd and report back?
>>> 
>>> Somebody may have your answer, but the system blowing up 
>> question is so common that it is probably faster for you to rule 
>> out some things first.
>>> 
>>> Chris.
>>> 
>>> -- original message --
>>> 
>>>  	I am looking at the the free energy profile of end-to 
>> end distances of various peptides, but I am consistently getting 
>> a system blow-up when running simulations with that distance 
>> constrained by a type 2 constraint.
>>> 
>>>  	I run a simulation of the unconstrained peptide in a box 
>> of tip4p water, Na+ and Cl- ions, and it runs with no problem.
>>> 
>>>  	Then I grab a frame of that simulation in which the end-
>> to-end distance is 0.8 nm (full frame including water) as a .gro 
>> file.  Then I add a type 2 constraint between the N-
>> terminal nitrogen atom and the C-terminal carboxyl carbon, 
>> create a new .tpr file using the revised topology and the 
>> already equilibrated frame as starting files and a .mdp file 
>> that now has the free_energy set to on, and then launch mdrun.
>>> 
>>>  	About 189 ps into the simulation I start getting 
>> warnings as follows, ultimately leading to blow-up and the run 
>> failing (fragment of error file output shown at end of message).
>>> 
>>>  	I am obviously missing something about how the 
>> constraint is handled.  Is it possible that this is a 
>> problem that arises because of domain decomposition over 
>> multiple nodes, ir is there something more basic that needs to 
>> be dealt with when imposing a type 2 constraint?
>>> 
>>> Warren Gallin
>>> 
>>> Step 94636, time 189.272 (ps)  LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 0.041028, max 0.161215 (between atoms 217 and 219)
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2  angle  previous, current, constraint length
>>> 
>>> Step 94637, time 189.274 (ps)  LINCS WARNING
>>> relative constraint deviation after LINCS:
>>> rms 0.018884, max 0.100333 (between atoms 217 and 218)
>>> bonds that rotated more than 30 degrees:
>>> atom 1 atom 2  angle  previous, current, constraint length
>>> 
>>> <SNIP>
>>> 
>>> --
>>> gmx-users mailing list    gmx-users at gromacs.org
>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>> Please search the archive at http://www.gromacs.org/search 
>> before posting!
>>> Please don't post (un)subscribe requests to the list. Use the
>>> www interface or send it to gmx-users-request at gromacs.org.
>>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>>> 
>> 
>> -- 
>> gmx-users mailing list    gmx-users at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> Please search the archive at http://www.gromacs.org/search 
>> before posting!
>> Please don't post (un)subscribe requests to the list. Use the 
>> www interface or send it to gmx-users-request at gromacs.org.
>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
> --
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>