<div dir="ltr">Hi Pavel,<div><br></div><div>Molecular Dynamics is a chaotic process, so no matter what precision you use the results will diverge - it will just take a bit longer in double precision.</div><div><br></div><div>The reason why the method works is that we rely on sampling phase space (with realistic dynamics) rather than making an infinite-precision prediction of how a particle/atom will move - remember that we anyway won't know the initial conditions (in particular atom velocities) perfectly. For this reason, single precision is virtually always accurate enough, apart from a few places where we do summation of forces. One good way to test this is to check the extent to which the total energy is conserved.</div><div><br></div><div>For the dihedral formula, IIRC the reason we changed it is that we wanted to be able to use so-called "coarse-grained" force fields. For normal atomic force fields the four atoms can *never* be on a straight line (since the angle potentials are too large and will prevent that), but for coarse-grained force fields the dihedral angles won't correspond directly to four bonded atoms, and then they can all be very close to linear once in a blue moon.</div><div><br></div><div>In that case, the scalar product between the two vectors will be *very* close to 1 or -1. Both the way you calculate the scalar product (where you might add a very small term to a large one) and the behavior of the acos() function (whose derivative approaches infinity at -1 and +1) can lead to very large numerical errors.</div><div><br></div><div>Cheers,</div><div><br></div><div>Erik</div><div><br></div><div><br></div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Sun, Oct 29, 2017 at 12:00 AM, Pavel Panchekha <span dir="ltr"><<a href="mailto:me@pavpanchekha.com" target="_blank">me@pavpanchekha.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_default" style="font-size:small">Berk, thank you for the reply. The quality of the GROMACS source code is clear from reading it.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Thank you for explaining that the sign is not used except in computing the angle. Do you know how much error in the angle is acceptable? I am trying some variations on the dihedral angle computation and seeing answers that differ by 1e-8. Could that small a difference possibly matter?<br></div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">I've noticed that earlier versions of GROMACS used a slightly different formula to compute <u>dih_angle</u> (using <u>acos</u> instead of <u>atan2</u>). Did the switch cause any known changes in behavior?<br></div></div><div class="gmail_extra"><br clear="all"><div><div class="m_-5706029853271738848gmail_signature" data-smartmail="gmail_signature"><div dir="ltr">—Pavel Panchekha</div></div></div>
<br><div class="gmail_quote"><div><div class="h5">On Sat, Oct 28, 2017 at 1:06 AM, Berk Hess <span dir="ltr"><<a href="mailto:hess@kth.se" target="_blank">hess@kth.se</a>></span> wrote:<br></div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div class="h5">
<div text="#000000" bgcolor="#FFFFFF">
<div class="m_-5706029853271738848m_3495107603991124129moz-cite-prefix">Hi,<br>
<br>
In molecular dynamics it's rather likely that you hit corner
cases, so we pay a lot of attention to numerical issues when
writing and reviewing code. Therefore it is unlikely that there
are cases that can be improved, but you are welcome to look.<br>
<br>
It is good that you mention the sign in dihedral angle. The return
value is actually never used. The sign is only used to set the
sign of the angle, which switches at angle=0, so there is no
stability issue.<br>
<br>
Cheers,<br>
<br>
Berk<div><div class="m_-5706029853271738848h5"><br>
<br>
On 10/28/2017 07:16 AM, Pavel Panchekha wrote:<br>
</div></div></div><div><div class="m_-5706029853271738848h5">
<blockquote type="cite">
<div dir="ltr">
<div class="gmail_default" style="font-size:small">Hello,</div>
<div class="gmail_default" style="font-size:small"><br>
</div>
<div class="gmail_default" style="font-size:small">I'm a student
at the University of Washington working on numerical accuracy,
and I've been looking into GROMACS as an example of a beloved
and important large numerical program. Browsing through the
code base, I was struck by the computation of the dihedral
angle in <u>bonded.c</u>, at line 1569 (in GROMACS 5.1.4).</div>
<div class="gmail_default" style="font-size:small"><br>
</div>
<div class="gmail_default" style="font-size:small">I wonder if
it is possible for the <u>sign</u> variable to get the wrong
value, if the input points are close to one another or have
some other reason to experience cancellation when doing the
cross product. I think this is fixable using an adaptive
algorithm. Could changing the sign of the angle computed by <u>dih_angle</u>
affect GROMACS's results?<br>
</div>
<div class="gmail_default" style="font-size:small"><br>
</div>
<div class="gmail_default" style="font-size:small">More broadly,
has anyone ever had numerical accuracy issues with GROMACS?
Would there be interest in reviewing a patch that improved the
accuracy (for edge case inputs) for the dihedral angle
computation?<br>
</div>
<div class="gmail_default" style="font-size:small"><br clear="all">
</div>
<div>
<div class="m_-5706029853271738848m_3495107603991124129m_3853295731210088818gmail_signature" data-smartmail="gmail_signature">
<div dir="ltr">—Pavel Panchekha</div>
</div>
</div>
</div>
<br>
<fieldset class="m_-5706029853271738848m_3495107603991124129mimeAttachmentHeader"></fieldset>
<br>
</blockquote>
<br>
</div></div></div>
<br></div></div><span class="HOEnZb"><font color="#888888">--<br>
Gromacs Developers mailing list<br>
<br>
* Please search the archive at <a href="http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List" rel="noreferrer" target="_blank">http://www.gromacs.org/Support<wbr>/Mailing_Lists/GMX-developers_<wbr>List</a> before posting!<br>
<br>
* Can't post? Read <a href="http://www.gromacs.org/Support/Mailing_Lists" rel="noreferrer" target="_blank">http://www.gromacs.org/Support<wbr>/Mailing_Lists</a><br>
<br>
* For (un)subscribe requests visit<br>
<a href="https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers" rel="noreferrer" target="_blank">https://maillist.sys.kth.se/ma<wbr>ilman/listinfo/gromacs.org_gmx<wbr>-developers</a> or send a mail to <a href="mailto:gmx-developers-request@gromacs.org" target="_blank">gmx-developers-request@gromacs<wbr>.org</a>.<br></font></span></blockquote></div><br></div>
<br>--<br>
Gromacs Developers mailing list<br>
<br>
* Please search the archive at <a href="http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List" rel="noreferrer" target="_blank">http://www.gromacs.org/<wbr>Support/Mailing_Lists/GMX-<wbr>developers_List</a> before posting!<br>
<br>
* Can't post? Read <a href="http://www.gromacs.org/Support/Mailing_Lists" rel="noreferrer" target="_blank">http://www.gromacs.org/<wbr>Support/Mailing_Lists</a><br>
<br>
* For (un)subscribe requests visit<br>
<a href="https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers" rel="noreferrer" target="_blank">https://maillist.sys.kth.se/<wbr>mailman/listinfo/gromacs.org_<wbr>gmx-developers</a> or send a mail to <a href="mailto:gmx-developers-request@gromacs.org">gmx-developers-request@<wbr>gromacs.org</a>.<br></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>Erik Lindahl <<a href="mailto:erik.lindahl@dbb.su.se" target="_blank">erik.lindahl@dbb.su.se</a>></div><div>Professor of Biophysics, Dept. Biochemistry & Biophysics, Stockholm University</div><div>Science for Life Laboratory, Box 1031, 17121 Solna, Sweden</div></div></div></div></div></div></div></div></div>
</div>