<br><br>----- Original Message -----<br>From: Warren Gallin &lt;wgallin@ualberta.ca&gt;<br>Date: Saturday, June 5, 2010 7:02<br>Subject: Re: [gmx-users] Constraint causing system to explode<br>To: Discussion list for GROMACS users &lt;gmx-users@gromacs.org&gt;<br><br>&gt; OK, I've done as Mark suggested, replaced the type 2 constraint <br>&gt; between the amino terminal nitrogen atom and the C-terminal <br>&gt; carbon atom with a type 6 bond, r0=0.8, r1=0.8, r2=1, which acts <br>&gt; as a distance restraint.<br>&gt; <br>&gt; This has run clean without crashing for 1.6 ns of simulation on <br>&gt; a single processor, so my explosion is not originating from the <br>&gt; linkage between the ends of the peptide per se.<br>&gt; <br>&gt; In this case the distance between the the two atoms increases <br>&gt; from the initial 0.8 nm to what appears to be an equilibrium <br>&gt; oscillating around 1.5 nm, which is what I would expect - the <br>&gt; most probable distance between the ends of the peptide is much <br>&gt; larger than 0.8 nm, so it's essentially pushing the spring.<br>&gt; <br>&gt; Which leaves me with the next question:&nbsp; What is it about <br>&gt; the type 2 constraint that is making things blow up?<br><br>I suspect the issue is that coupled constraints are hard to implement reliably. If you were using constraints = all-bonds, then at every integration step you have a cycle of constraints, which might not be a happy numerical situation if there's an appreciable force every step to increase your end-to-end constraint. One can get weird resonance effects building up (think Tacoma Narrows). You might be able to play with the constraint algorithm settings to give it enough flex to work, or try constraints = h-bonds. If you're not using constraints = all-bonds then the above might still be worth trying, although I haven't got as good a rationale for it. A smaller integration step might work, also.<br>&nbsp;<br>&gt; Intuitively it seems possible that after an unusual step of <br>&gt; simulation a water molecule gets between the two ends or atoms <br>&gt; attached to the two constrained atom rotate into line between <br>&gt; the two constrained atoms, and when the constraint pulls the <br>&gt; ends back to 0.8 nm&nbsp; it creates a ridiculous force.<br>&gt; <br>&gt; Is that a valid worry?<br><br>Probably not in a single MD step - if the end-to-end repulsive force is strong enough to create a void into which a water could move later before being smashed, then the ends would have been smashing up other stuff beforehand.<br><br>&gt; But why would that be true only for a type 2 constraint between <br>&gt; the ends of the peptide, and not all the other bond constraints <br>&gt; in the rest of the peptide?&nbsp; Should I be considering the <br><br>See above.<br><br>&gt; use of a type 1 constraint (with exclusions), and if so, is that <br>&gt; formulation legitimate for determining the free energy?<br><br>The presence or absence of exclusions is not going to be part of the issue. The purpose of exclusions in the model physics is to create a particular balance between general bonded and non-bonded terms. Your desire for an end-to-end constraint isn't going to affect this.<br><br>Mark<br><br>&gt; On 2010-05-31, at 10:09 PM, Mark Abraham wrote:<br>&gt; <br>&gt; &gt; ----- Original Message -----<br>&gt; &gt; From: Warren Gallin &lt;wgallin@ualberta.ca&gt;<br>&gt; &gt; Date: Tuesday, June 1, 2010 13:38<br>&gt; &gt; Subject: Re: [gmx-users] Constraint causing system to explode<br>&gt; &gt; To: Discussion list for GROMACS users &lt;gmx-users@gromacs.org&gt;<br>&gt; &gt; <br>&gt; &gt;&gt; Chris,<br>&gt; &gt;&gt; <br>&gt; &gt;&gt;&nbsp;        As you suggested, I ran it on one node, and it still <br>&gt; blew up, <br>&gt; &gt;&gt; so it appears to be some oddity with the constraint that I am <br>&gt; &gt;&gt; imposing.<br>&gt; &gt; <br>&gt; &gt; Does it happen with a distance restraint, instead? (Is that a <br>&gt; better physical model?) Is there any possible correlation with <br>&gt; one atom crossing a PBC boundary, such that a buggy <br>&gt; implementation would then be horribly broken on the next step?<br>&gt; &gt; <br>&gt; &gt; Mark<br>&gt; &gt; <br>&gt; &gt;&gt;&nbsp;        When I look at the dgdl.xvg file from the run, there <br>&gt; appear to <br>&gt; &gt;&gt; be a number of transient spikes that are a couple orders of <br>&gt; &gt;&gt; magnitude larger than the typical variation during the run, <br>&gt; and <br>&gt; &gt;&gt; I am getting errors on an angle changing by more than 30 <br>&gt; degrees <br>&gt; &gt;&gt; as a LINCS warning during those short spikes.&nbsp; The final <br>&gt; &gt;&gt; system blow-up starts with the same situtation, which <br>&gt; &gt;&gt; degenerates into a full blow system explosion, with bond <br>&gt; lengths <br>&gt; &gt;&gt; going to crzy high values.<br>&gt; &gt;&gt; <br>&gt; &gt;&gt;&nbsp;        So my question boils down to this:&nbsp; What is it <br>&gt; about a <br>&gt; &gt;&gt; type 2 constraint that will generate these ridiculously large <br>&gt; &gt;&gt; forces on relatively rare occasions?&nbsp; Or perhaps, what <br>&gt; &gt;&gt; error am I making in setting up the simulation that would <br>&gt; cause <br>&gt; &gt;&gt; only the constrained system to manifest this behavior?<br>&gt; &gt;&gt; <br>&gt; &gt;&gt; Warren Gallin<br>&gt; &gt;&gt; <br>&gt; &gt;&gt; <br>&gt; &gt;&gt; On 2010-05-27, at 11:22 AM, chris.neale@utoronto.ca wrote:<br>&gt; &gt;&gt; <br>&gt; &gt;&gt;&gt; Dear Warren:<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt; I don't have your answer, but I'll point out that when you <br>&gt; &gt;&gt; ask: "Is it possible that this is a problem that arises <br>&gt; because <br>&gt; &gt;&gt; of domain decomposition over multiple nodes" that you are <br>&gt; &gt;&gt; probably the person in the best position to address this. <br>&gt; 300ps <br>&gt; &gt;&gt; should not take too long to simulate so why not try it on a <br>&gt; &gt;&gt; single node and also on multiple nodes with mdrun -pd and <br>&gt; report back?<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt; Somebody may have your answer, but the system blowing up <br>&gt; &gt;&gt; question is so common that it is probably faster for you to <br>&gt; rule <br>&gt; &gt;&gt; out some things first.<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt; Chris.<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt; -- original message --<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt;&nbsp;&nbsp;        I am looking at the the free energy profile of <br>&gt; end-to <br>&gt; &gt;&gt; end distances of various peptides, but I am consistently <br>&gt; getting <br>&gt; &gt;&gt; a system blow-up when running simulations with that distance <br>&gt; &gt;&gt; constrained by a type 2 constraint.<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt;&nbsp;&nbsp;        I run a simulation of the unconstrained peptide <br>&gt; in a box <br>&gt; &gt;&gt; of tip4p water, Na+ and Cl- ions, and it runs with no problem.<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt;&nbsp;&nbsp;        Then I grab a frame of that simulation in which <br>&gt; the end-<br>&gt; &gt;&gt; to-end distance is 0.8 nm (full frame including water) as a <br>&gt; .gro <br>&gt; &gt;&gt; file.&nbsp; Then I add a type 2 constraint between the N-<br>&gt; &gt;&gt; terminal nitrogen atom and the C-terminal carboxyl carbon, <br>&gt; &gt;&gt; create a new .tpr file using the revised topology and the <br>&gt; &gt;&gt; already equilibrated frame as starting files and a .mdp file <br>&gt; &gt;&gt; that now has the free_energy set to on, and then launch mdrun.<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt;&nbsp;&nbsp;        About 189 ps into the simulation I start getting <br>&gt; &gt;&gt; warnings as follows, ultimately leading to blow-up and the <br>&gt; run <br>&gt; &gt;&gt; failing (fragment of error file output shown at end of message).<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt;&nbsp;&nbsp;        I am obviously missing something about how the <br>&gt; &gt;&gt; constraint is handled.&nbsp; Is it possible that this is a <br>&gt; &gt;&gt; problem that arises because of domain decomposition over <br>&gt; &gt;&gt; multiple nodes, ir is there something more basic that needs <br>&gt; to <br>&gt; &gt;&gt; be dealt with when imposing a type 2 constraint?<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt; Warren Gallin<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt; Step 94636, time 189.272 (ps)&nbsp; LINCS WARNING<br>&gt; &gt;&gt;&gt; relative constraint deviation after LINCS:<br>&gt; &gt;&gt;&gt; rms 0.041028, max 0.161215 (between atoms 217 and 219)<br>&gt; &gt;&gt;&gt; bonds that rotated more than 30 degrees:<br>&gt; &gt;&gt;&gt; atom 1 atom 2&nbsp; angle&nbsp; previous, current, <br>&gt; constraint length<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt; Step 94637, time 189.274 (ps)&nbsp; LINCS WARNING<br>&gt; &gt;&gt;&gt; relative constraint deviation after LINCS:<br>&gt; &gt;&gt;&gt; rms 0.018884, max 0.100333 (between atoms 217 and 218)<br>&gt; &gt;&gt;&gt; bonds that rotated more than 30 degrees:<br>&gt; &gt;&gt;&gt; atom 1 atom 2&nbsp; angle&nbsp; previous, current, <br>&gt; constraint length<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt; &lt;SNIP&gt;<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt;&gt; --<br>&gt; &gt;&gt;&gt; gmx-users mailing list&nbsp;&nbsp;&nbsp; gmx-users@gromacs.org<br>&gt; &gt;&gt;&gt; http://lists.gromacs.org/mailman/listinfo/gmx-users<br>&gt; &gt;&gt;&gt; Please search the archive at http://www.gromacs.org/search <br>&gt; &gt;&gt; before posting!<br>&gt; &gt;&gt;&gt; Please don't post (un)subscribe requests to the list. Use the<br>&gt; &gt;&gt;&gt; www interface or send it to gmx-users-request@gromacs.org.<br>&gt; &gt;&gt;&gt; Can't post? Read http://www.gromacs.org/mailing_lists/users.php<br>&gt; &gt;&gt;&gt; <br>&gt; &gt;&gt; <br>&gt; &gt;&gt; -- <br>&gt; &gt;&gt; gmx-users mailing list&nbsp;&nbsp;&nbsp; gmx-users@gromacs.org<br>&gt; &gt;&gt; http://lists.gromacs.org/mailman/listinfo/gmx-users<br>&gt; &gt;&gt; Please search the archive at http://www.gromacs.org/search <br>&gt; &gt;&gt; before posting!<br>&gt; &gt;&gt; Please don't post (un)subscribe requests to the list. Use the <br>&gt; &gt;&gt; www interface or send it to gmx-users-request@gromacs.org.<br>&gt; &gt;&gt; Can't post? Read http://www.gromacs.org/mailing_lists/users.php<br>&gt; &gt; --<br>&gt; &gt; gmx-users mailing list&nbsp;&nbsp;&nbsp; gmx-users@gromacs.org<br>&gt; &gt; http://lists.gromacs.org/mailman/listinfo/gmx-users<br>&gt; &gt; Please search the archive at http://www.gromacs.org/search <br>&gt; before posting!<br>&gt; &gt; Please don't post (un)subscribe requests to the list. Use the<br>&gt; &gt; www interface or send it to gmx-users-request@gromacs.org.<br>&gt; &gt; Can't post? Read http://www.gromacs.org/mailing_lists/users.php<br>&gt; &gt; <br>&gt; <br>&gt; -- <br>&gt; gmx-users mailing list&nbsp;&nbsp;&nbsp; gmx-users@gromacs.org<br>&gt; http://lists.gromacs.org/mailman/listinfo/gmx-users<br>&gt; Please search the archive at http://www.gromacs.org/search <br>&gt; before posting!<br>&gt; Please don't post (un)subscribe requests to the list. Use the <br>&gt; www interface or send it to gmx-users-request@gromacs.org.<br>&gt; Can't post? Read http://www.gromacs.org/mailing_lists/users.php