<html>
<head>
<style>
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 10pt;
font-family:Verdana
}
</style>
</head>
<body class='hmmessage'>
Hi,<br><br>I don't know what is going wrong in your calculations.<br><br>I repeated your calculations with a script for many distances<br>and the forces perfectly match the numerical derivative of the potential.<br><br>Example: distance, pot, force:<br>3.998 -40.569943757001 -5.14772e+00<br>4.000 -40.559660943670 -5.13585e+00<br>4.002 -40.549401856871 -5.12399e+00<br><br>The numerical derivative of pot at d=4 (a=2) gives 5.135475,<br>which gives a difference of 0.0004 kJ/mol.<br>I expect this difference to go down if I take the points closer together.<br><br>Berk<br><br>> Date: Sun, 14 Jun 2009 19:11:08 +0200<br>> From: dommert@icp.uni-stuttgart.de<br>> To: gmx-users@gromacs.org<br>> Subject: [gmx-users] Electrostatic forces<br>> <br>> Hello Gromacs users and developers,<br>> <br>> I have a few questions concerning the calculation of the electrostatic<br>> forces in Gromacs.<br>> <br>> Unfortunately I am not able to reproduce the electrostatic forces<br>> computed by the gromacs code ( VERSION 4.0.99_DEAD_CVS_BRANCH_20090604<br>> (double precision), compiled on MacOS X 10.5 using its native gcc and<br>> no fortran compiler ) by a self written code , that also uses double<br>> precision, though the energies derived with the different codes match<br>> perfectly. To get an idea where the deviation stems from Berk already<br>> gave me the hint to setup following configuration:<br>> <br>> A positive charge is located at (5.0+a,5.0,5.0) and a negative charge<br>> at (5.0-a,5.0,5.0) in a cubic box of size 10x10x10 and this setup for<br>> the electrostatic forces with an arbitrary parameter a, that allows me<br>> to control the strength of the interaction. Furthermore following<br>> mdp-parameter were applied to achieve a very high accuracy of the Ewald<br>> sum ( delta < 10^-18 kJ / (mol nm) compared to an infinite large cutoff ):<br>> <br>> rlist=4.5<br>> rcoulomb=4.5<br>> rvdw=4.5<br>> <br>> coulomb-type=ewald<br>> knx=kny=knz=70<br>> ewald_rtol=1.348768e-21 ( -> 1/beta=0.6666667 )<br>> <br>> Now I checked different cases:<br>> <br>> 1. a = 2.5 -> only the reciprocal part of the Ewald sum is taken into<br>> account and we have small terms in the sums:<br>> <br>> In this case everything is perfect, apart from the forces in y and z<br>> direction. The forces as well as the energies match up the last<br>> digit.<br>> <br>> 2. a=2.0 -> both contribtions from direct and recipr. space and the<br>> amount of the terms is much larger.<br>> <br>> The energies also match perfectly, however the forces differ<br>> already by an amount of 10^-2 kJ / (mol nm).<br>> <br>> 3. a=0.2 -> very large terms<br>> <br>> Furthermore the energies fit perfectly, but the forces differ<br>> significantly.<br>> <br>> As a first resumee I have assumed the implementation of the calculation<br>> of the direct forces is wrong. Than I checked the code and compared the<br>> analog part in gromacs. To stay as compatible as much is also use<br>> gmx_erfc() in my code.<br>> <br>> The energies are derived by summing over all pairs within the cutoff.<br>> In this routine also the forces are calculated and added to the involved <br>> particles.<br>> <br>> Unfortunately there are so many sources where an an error can easily be<br>> introduce. So at first I checked terms get lost during the summation of the<br>> forces. But then I would not obtain exactly the same energies like<br>> gromacs. To exclude the possibility of a wrong calculation of the<br>> direct forces I already compared the code lines for the force<br>> calculation in gromacs (src/mdlib/tables.c line 655) and the other<br>> code, but they also match. <br>> <br>> <br>> Can somebody give me an idea, where this deviations can arise from ?<br>> Obviously an increase of the forces yields significant differences.<br>> This can also be observed, when I decrease the box size and test case<br>> 1.). For a 4x4x4 box the deviation is as large as the actual force.<br>> <br>> Cheers,<br>> <br>> Flo<br>> <br>> Florian Dommert<br>> Dipl.-Phys.<br>> <br>> Institute for Computational Physics<br>> University Stuttgart<br>> <br>> Pfaffenwaldring 27<br>> 70569 Stuttgart<br>> <br>> Tel: +49 - 711 / 6856-3613<br>> Fax: +49 - 711 / 6856-3658<br>> <br>> EMail: dommert@icp.uni-stuttgart.de<br>> Home: http://www.icp.uni-stuttgart.de/~icp/Florian_Dommert<br>> <br>> !! PGP-ENCODED emails preferred !!<br><br /><hr />What can you do with the new Windows Live? <a href='http://www.microsoft.com/windows/windowslive/default.aspx' target='_new'>Find out</a></body>
</html>