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: domm...@icp.uni-stuttgart.de
Home: http://www.icp.uni-stuttgart.de/~icp/Florian_Dommert

!! PGP-ENCODED emails preferred !!

Attachment: pgpNeWfhD3tEY.pgp
Description: PGP signature

_______________________________________________
gmx-users mailing list    gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
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-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php

Reply via email to