[gmx-users] Electrostatic forces
Berk Hess
gmx3 at hotmail.com
Tue Jun 16 13:56:46 CEST 2009
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>
More information about the gromacs.org_gmx-users
mailing list