<link rel="self" href="https://redmine.gromacs.org/issues/172.atom"/>
<link rel="alternate" href="https://redmine.gromacs.org/"/>
<id>https://redmine.gromacs.org/</id>
<icon>https://redmine.gromacs.org/favicon.ico?1402403081</icon>
<updated>2007-10-17T06:23:31Z</updated>
<author>
<name>GROMACS development</name>
</author>
<entry>
<title>GROMACS - Bug #172: incorrect handling of angle bending leads to errors for obtuse angles and crashes for angles of 180https://redmine.gromacs.org/issues/172?journal_id=10342007-10-17T06:23:31ZDavid van der Spoelspoel@xray.bmc.uu.se
<ul></ul><p>If the solution is simple then please provide code or at least equations or preferably both.</p> GROMACS - Bug #172: incorrect handling of angle bending leads to errors for obtuse angles and crashes for angles of 180https://redmine.gromacs.org/issues/172?journal_id=10352007-10-17T06:59:25ZJohn D.jchodera@zimm.compbio.ucsf.edu
<ul></ul><p>As far as my meagre brain can grasp, the angle energy terms and gradient are computed in gmxlib/bondfree.c, in a function called 'angles'. The current torsion angle is computed by a call to bond_angle(), which computes the cos of the angle between the ij and ik bond vectors, and then calls 'acos' to compute the inverse cosine, which returns the angle theta in the range [0,pi]. The energy and gradient is then computed by a call to 'harmonic', which simply computes the linearly interpolated harmonic contributes by computing the deviation from the reference angle (let's call it theta0).</p>
<p>Keeping either of these angle terms in David's .top file would cause David Mobley's vacuum calculations to segfault:</p>
<p>[ angles ]<br />; ai aj ak funct theta cth<br /> 2 1 4 1 1.7838e+02 3.7489e+02 ;<br /> 1 2 3 1 1.7799e+02 4.7196e+02 ;</p>
<p>Note that both reference angles are nearly 180 degrees. They run fine in AMBER, apparently.</p>
<p>The problem is more subtle than I first appreciated. At first, I thought it was simply that the range of the result of acos(cos_theta) would need to be shifted, but once theta opens up to pi, it is impossible to tell that it is continuing to open up larger than pi, and it looks like the angle is shrinking again.</p>
<p>AMBER seems to ensure that cos_theta remains on the domain [-0.999,+0.999] before feeding it to acos(cos_theta) by clipping. I imagine that numerical errors in the computation of cos_theta in gromacs might be causing values that exceed the defined domain to be fed to acos(), and that could be causing the observed problems, but I haven't experimented to figure this out.</p>
<p>I apologize for my earlier remarks that this bug could be fixed by just shifting the range of the output of acos() -- it looks to be more subtle, and hopefully simply a matter of ensuring the input to acos() remains in this domain.</p>
<p>It should also be noted that the derivative is formally discontinuous at theta = pi here, but I don';t think there's any way to fix that except to move to force terms that are harmonic in cos(theta) rather than theta, but that's a forcefield problem rather than a numerical stability of the implementation problem.</p>
<p>I'll upload David Mobley's example in a subsequent mail.</p>
<p>Cheers,</p>
<p>- John</p> GROMACS - Bug #172: incorrect handling of angle bending leads to errors for obtuse angles and crashes for angles of 180https://redmine.gromacs.org/issues/172?journal_id=10362007-10-17T07:08:02ZJohn D.jchodera@zimm.compbio.ucsf.edu
<ul></ul><p>Created an attachment (id=250)<br />propyne in vacuum test case from David Mobley.</p>
<p>Added propyne in vacuum test case from David Mobley. Unpack and run ./run0.sh to run, which should crash during dynamics phase. The first and last members of the [ angles ] section seem to be responsible for the crash.</p> GROMACS - Bug #172: incorrect handling of angle bending leads to errors for obtuse angles and crashes for angles of 180https://redmine.gromacs.org/issues/172?journal_id=10392007-10-17T10:01:10ZDavid van der Spoelspoel@xray.bmc.uu.se
<ul></ul><blockquote>
<p>AMBER seems to ensure that cos_theta remains on the domain [-0.999,+0.999]<br />before feeding it to acos(cos_theta) by clipping. I imagine that numerical<br />errors in the computation of cos_theta in gromacs might be causing values that<br />exceed the defined domain to be fed to acos(), and that could be causing the<br />observed problems, but I haven't experimented to figure this out.</p>
</blockquote>
<p>So the real problem is what to do when cos_theta falls outside the domain! Suggestions? What does AMBER do in these cases? Set the force to zero? In that case it probably becomes discontinuous.</p>
<p>It's not that we haven't thought of this before, we have, but for proteins it is not a problem.</p> GROMACS - Bug #172: incorrect handling of angle bending leads to errors for obtuse angles and crashes for angles of 180https://redmine.gromacs.org/issues/172?journal_id=10402007-10-17T14:40:06ZDavid Mobleydmobley@gmail.com
<ul></ul><p>I've been doing a little more testing, and, perhaps not surprisingly, this all works basically fine in double precision, presumably because the extra precision means that cos_theta never hits exactly +1 or -1. Of course, that's not really a solution if you want to be able to use single precision.</p>
<p>Maybe clipping as AMBER does would be a good solution.</p> GROMACS - Bug #172: incorrect handling of angle bending leads to errors for obtuse angles and crashes for angles of 180https://redmine.gromacs.org/issues/172?journal_id=10432007-10-17T17:25:10ZJohn D.jchodera@zimm.compbio.ucsf.edu
<ul></ul><p>Hi guys,</p>
<p>Could you simply use double-precision accumulators for the computation of cos_theta and acos(cos_theta)? If David Mobley is right and cos_theta never hits exactly +1 or -1, then this should fix the problem.</p>
<p>Clipping, as in AMBER, doesn't seem like the optimal solution, though it at least avoids the problem of feeding acos() an argument outside of its defined domain.</p>
<p>Cheers,</p>
<p>- John</p> GROMACS - Bug #172: incorrect handling of angle bending leads to errors for obtuse angles and crashes for angles of 180https://redmine.gromacs.org/issues/172?journal_id=10442007-10-17T17:28:24ZErik Lindahllindahl@sbc.su.se
<ul></ul><p>But will double precision really solve it, rather than just making it a couple of orders of magnitude more improbable?</p>
<p>Cheers,</p>
<p>Erik</p> GROMACS - Bug #172: incorrect handling of angle bending leads to errors for obtuse angles and crashes for angles of 180https://redmine.gromacs.org/issues/172?journal_id=10452007-10-17T17:38:28ZDavid Mobleydmobley@gmail.com
<ul></ul><p>Erik, I worried about that too. What I can say at this point is it seems more than a couple orders of magnitude better in double precision, in that before it was crashing reliably in single precision within a few hundred ps and so far I've run ~30 5 ns simulations in double precision with no crashes.</p>
<p>My intuition is that you are correct and this is not really a fix, but would simply make the problem a lot more rare. Maybe rare enough, though (it seems to be making it rare enough for my purposes).</p>
<p>As John points out, clipping is also less than ideal, though.</p>
<p>My suggestion is to either:<br />(a) Go with some sort of clipping because it's what AMBER does and at least prevents the problem<br />(b) Go to double precision here and put a note in the documentation in case it ever comes up again, in which case then add some (very slight) clipping in double precision.</p> GROMACS - Bug #172: incorrect handling of angle bending leads to errors for obtuse angles and crashes for angles of 180https://redmine.gromacs.org/issues/172?journal_id=10532007-10-20T07:56:24ZDavid van der Spoelspoel@xray.bmc.uu.se
<ul></ul><p>I have added a patch in routine bond_angle like this:</p>
<pre><code><strong>costh=cos_angle(r_ij,r_kj); /</strong> 25 */<br /> if (*costh == -1.0)<br /> *costh += GMX_REAL_EPS;<br /> else if (*costh == 1.0)<br /> *costh -= GMX_REAL_EPS;</code></pre>
<p>it seems to solve the problem. If you'd looked in include/vec.h you'd see that we already had tests in cos_angle for costh > 1 or costh < -1, however it is the problem that angles get exactly 180 or 0 that causes these problems.</p>
<p>Shall I commit this to the 3.3 cvs branch?</p>
<p>By the way David, you are really good at setting up "troublesome" systems: I also had to patch grompp to be able to get your topology grompp-ed, since there is a different of constraint generation in 3.3.1 and 3.3.2. Oh well, all the better.</p> GROMACS - Bug #172: incorrect handling of angle bending leads to errors for obtuse angles and crashes for angles of 180https://redmine.gromacs.org/issues/172?journal_id=10542007-10-20T08:13:15ZJohn D.jchodera@zimm.compbio.ucsf.edu
<ul></ul><p>Thanks for the quick diagnosis and fix, David!</p>
<p>You're right -- I had missed the rectification conditions to ensure that cos_theta was ensured of being in the domain [-1,+1].</p>
<p>Cheers,</p>
<p>- John</p> GROMACS - Bug #172: incorrect handling of angle bending leads to errors for obtuse angles and crashes for angles of 180https://redmine.gromacs.org/issues/172?journal_id=10562007-10-22T16:28:48ZDavid Mobleydmobley@gmail.com
<ul></ul><p>David, committing this to CVS sounds like a good idea -- it sounds like your solution should do it.</p>
<p>As far as setting up troublesome systems, I'm always happy to oblige. Can I get paid twice my usual rate for finding two problems with one Bugzilla? :)</p>
<p>David</p> GROMACS - Bug #172: incorrect handling of angle bending leads to errors for obtuse angles and crashes for angles of 180https://redmine.gromacs.org/issues/172?journal_id=10572007-10-24T10:13:26ZDavid van der Spoelspoel@xray.bmc.uu.se
<ul></ul><p>Fixed in CVS for 3.3 and 4.0. The fix is applied to all force calculations (angles dihedrals etc.)</p>