[gmx-developers] Dihedral angle numerical accuracy
Magnus Lundborg
magnus.lundborg at scilifelab.se
Mon Oct 30 08:50:44 CET 2017
Hi,
I'll just chime in and say that if there are triple bonds four atoms
will, with a rather high possibility, be on a straight line in atomic
force fields as well.
Cheers,
Magnus
On 2017-10-29 09:18, Erik Lindahl wrote:
> 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
> <mailto: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
> <mailto: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
> <http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List>
> before posting!
>
> * Can't post? Read
> http://www.gromacs.org/Support/Mailing_Lists
> <http://www.gromacs.org/Support/Mailing_Lists>
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers
> <https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers>
> or send a mail to gmx-developers-request at gromacs.org
> <mailto: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
> <http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List>
> before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> <http://www.gromacs.org/Support/Mailing_Lists>
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers
> <https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers>
> or send a mail to gmx-developers-request at gromacs.org
> <mailto:gmx-developers-request at gromacs.org>.
>
>
>
>
> --
> Erik Lindahl <erik.lindahl at dbb.su.se <mailto: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/20171030/9a19fdff/attachment.html>
More information about the gromacs.org_gmx-developers
mailing list