[gmx-users] Electrostatic forces
gmx3 at hotmail.com
Tue Jun 16 13:56:46 CEST 2009
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.
> 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
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the gromacs.org_gmx-users