[gmx-users] Electrostatic forces

[Berk Hess]

Hi,

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.

Berk

> 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 ):
> 
>     rlist=4.5
>     rcoulomb=4.5
>     rvdw=4.5
> 
>     coulomb-type=ewald
>     knx=kny=knz=70
>     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
>       digit.
> 
>     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
>       significantly.
> 
>    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 
>    particles.
> 
>    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.
> 
> Cheers,
> 
> Flo
> 
> Florian Dommert
> Dipl.-Phys.
> 
> 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
http://www.microsoft.com/windows/windowslive/default.aspx
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090616/8faad5b9/attachment.html>