[gmx-developers] Dihedral angle numerical accuracy

[Erik Lindahl]

Hi Pavel,

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.

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.

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.

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.

Cheers,

Erik




On Sun, Oct 29, 2017 at 12:00 AM, Pavel Panchekha <me at pavpanchekha.com>
wrote:

> Berk, thank you for the reply. The quality of the GROMACS source code is
> clear from reading it.
>
> 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?
>
> I've noticed that earlier versions of GROMACS used a slightly different
> formula to compute *dih_angle* (using *acos* instead of *atan2*). Did the
> switch cause any known changes in behavior?
>
> —Pavel Panchekha
>
> On Sat, Oct 28, 2017 at 1:06 AM, Berk Hess <hess at kth.se> wrote:
>
>> Hi,
>>
>> 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.
>>
>> 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.
>>
>> Cheers,
>>
>> Berk
>>
>>
>> On 10/28/2017 07:16 AM, Pavel Panchekha wrote:
>>
>> Hello,
>>
>> 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 *bonded.c*, at line
>> 1569 (in GROMACS 5.1.4).
>>
>> I wonder if it is possible for the *sign* 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 *dih_angle* affect GROMACS's results?
>>
>> 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?
>>
>> —Pavel Panchekha
>>
>>
>>
>>
>> --
>> Gromacs Developers mailing list
>>
>> * Please search the archive at http://www.gromacs.org/Support
>> /Mailing_Lists/GMX-developers_List before posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers
>> or send a mail to gmx-developers-request at gromacs.org.
>>
>
>
> --
> Gromacs Developers mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-developers_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers
> or send a mail to gmx-developers-request at gromacs.org.
>



-- 
Erik Lindahl <erik.lindahl at dbb.su.se>
Professor of Biophysics, Dept. Biochemistry & Biophysics, Stockholm
University
Science for Life Laboratory, Box 1031, 17121 Solna, Sweden
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20171029/e880aea6/attachment.html>