[gmx-users] Electrostatic forces
dommert at icp.uni-stuttgart.de
Wed Jun 17 06:05:16 CEST 2009
* Berk Hess <gmx3 at hotmail.com> [2009-06-16 13:56:46 +0200]:
>I don't know what is going wrong in your calculations.
>I repeated your calculations with a script for many distances
>and the forces perfectly match the numerical derivative of the potential.
>Example: distance, pot, force:
>3.998 -40.569943757001 -5.14772e+00
>4.000 -40.559660943670 -5.13585e+00
>4.002 -40.549401856871 -5.12399e+00
>The numerical derivative of pot at d=4 (a=2) gives 5.135475,
>which gives a difference of 0.0004 kJ/mol.
>I expect this difference to go down if I take the points closer together.
At first many thanks for performing the testsi and the advice. This procedure allowed
me to track down, where difference arises from. It showed the at small
distances my code got inconsitent, when comparing energies and forces.
Unfortunately I was using the gmx_erfc and it seems that for small
values of r I somehow got wrong values. I tried the different #define
statements but nothing changed considerable. So finally it was clear
that my summation method in the reciprocal part was buggy. However now I
fixed this by changing the summation. The effect is, that rounding
errors cancel out for k and -k and I achieve the same results like
gromacs. Now have to print out a high precision example to check how
good the precision hold, when reaching machine precision.
However I am very confident and in case of success, that there will be
soon an error estimate for the Ewald Sum available, which will be the first
step to the an implementation a tuning routine for the SPME paramters to achieve optimal
balance between performance and accuracy ;)
>> Date: Sun, 14 Jun 2009 19:11:08 +0200
>> From: dommert at icp.uni-stuttgart.de
>> To: gmx-users at gromacs.org
>> Subject: [gmx-users] Electrostatic forces
>> Hello Gromacs users and developers,
>> I have a few questions concerning the calculation of the electrostatic
>> forces in Gromacs.
>> Unfortunately I am not able to reproduce the electrostatic forces
>> computed by the gromacs code ( VERSION 4.0.99_DEAD_CVS_BRANCH_20090604
>> (double precision), compiled on MacOS X 10.5 using its native gcc and
>> no fortran compiler ) by a self written code , that also uses double
>> precision, though the energies derived with the different codes match
>> perfectly. To get an idea where the deviation stems from Berk already
>> gave me the hint to setup following configuration:
>> A positive charge is located at (5.0+a,5.0,5.0) and a negative charge
>> at (5.0-a,5.0,5.0) in a cubic box of size 10x10x10 and this setup for
>> the electrostatic forces with an arbitrary parameter a, that allows me
>> to control the strength of the interaction. Furthermore following
>> mdp-parameter were applied to achieve a very high accuracy of the Ewald
>> sum ( delta < 10^-18 kJ / (mol nm) compared to an infinite large cutoff ):
>> ewald_rtol=1.348768e-21 ( -> 1/beta=0.6666667 )
>> Now I checked different cases:
>> 1. a = 2.5 -> only the reciprocal part of the Ewald sum is taken into
>> account and we have small terms in the sums:
>> In this case everything is perfect, apart from the forces in y and z
>> direction. The forces as well as the energies match up the last
>> 2. a=2.0 -> both contribtions from direct and recipr. space and the
>> amount of the terms is much larger.
>> The energies also match perfectly, however the forces differ
>> already by an amount of 10^-2 kJ / (mol nm).
>> 3. a=0.2 -> very large terms
>> Furthermore the energies fit perfectly, but the forces differ
>> As a first resumee I have assumed the implementation of the calculation
>> of the direct forces is wrong. Than I checked the code and compared the
>> analog part in gromacs. To stay as compatible as much is also use
>> gmx_erfc() in my code.
>> The energies are derived by summing over all pairs within the cutoff.
>> In this routine also the forces are calculated and added to the involved
>> Unfortunately there are so many sources where an an error can easily be
>> introduce. So at first I checked terms get lost during the summation of the
>> forces. But then I would not obtain exactly the same energies like
>> gromacs. To exclude the possibility of a wrong calculation of the
>> direct forces I already compared the code lines for the force
>> calculation in gromacs (src/mdlib/tables.c line 655) and the other
>> code, but they also match.
>> Can somebody give me an idea, where this deviations can arise from ?
>> Obviously an increase of the forces yields significant differences.
>> This can also be observed, when I decrease the box size and test case
>> 1.). For a 4x4x4 box the deviation is as large as the actual force.
>> Florian Dommert
>> Institute for Computational Physics
>> University Stuttgart
>> Pfaffenwaldring 27
>> 70569 Stuttgart
>> Tel: +49 - 711 / 6856-3613
>> Fax: +49 - 711 / 6856-3658
>> EMail: dommert at icp.uni-stuttgart.de
>> Home: http://www.icp.uni-stuttgart.de/~icp/Florian_Dommert
>> !! PGP-ENCODED emails preferred !!
>What can you do with the new Windows Live? Find out
>gmx-users mailing list gmx-users at gromacs.org
>Please search the archive at http://www.gromacs.org/search before posting!
>Please don't post (un)subscribe requests to the list. Use the
>www interface or send it to gmx-users-request at gromacs.org.
>Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Institute for Computational Physics
Tel: +49 - 711 / 6856-3613
Fax: +49 - 711 / 6856-3658
EMail: dommert at icp.uni-stuttgart.de
!! PGP-ENCODED emails preferred !!
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Size: 194 bytes
Desc: not available
More information about the gromacs.org_gmx-users