[gmx-users] Ewald summation: n=0

[Steve Fiedler]

Gromacs Users:

While porting field parameters to Gromacs, I ran across unexpected 
electrostatic energies returned by the ewald summation method.  To 
simplify the problem, I reduced the system to two partially charged 
particles.  For the simplest case of only one box, n = 0, I expected to 
see contributions from the direct sum (Vdir) and self-interaction 
correction (V0) terms.  I can exactly match the analytical values for 
this system using DL_POLY, but the Gromacs Coulombic values obtained 
with gmx_dump appear to differ.  Most notably the presence of a 
canceling term (Vrec) was surprising in this case.

To obtain n = 0, I set fourier_nx = fourier_ny = fourier_nz = 0, since m 
= 2*pi*n/L^2 = 0. A subsequent gmx_dump yields:
Coulomb (SR):  -10.4417 kJ/mol
Coul. recip.  -63.6614 kJ/mol
Coulomb + Coul. recip. = -74.2207 kJ/mol.

My back of the envelope calculations are as follows:
Beta = 3.47046 nm^-1 (from md.log)
Vdir = -10.4416 kJ/mol, V0 = -2*43.5255 kJ/mol
V = Vdir + V0 = -97.4927 kJ/mol.  Again, this matches the output from a 
DL_POLY calculation, but differs from the -74.2207 kJ/mol from above.

As an experiment, I found a close, but inexact match to the analytical 
value by setting fourier_nx = fourier_ny = fourier_nz = 1. 
The below input files were used with Gromacs 4.0.4 in double precision. 

conf.gro
-----------
Two_Atoms
  2
    1DMP     AA    1   0.100   0.000   0.000
    1DMP     BB    2   0.400   0.000   0.000
  10.00000  10.00000  10.00000

topol.top
-----------
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1     2     no      0.125  1.0

[ atomtypes ]
;name     mass     charge ptype  sig           eps
   AA    10.0000   -0.400  A     0.39057      0.500
   BB    10.0000    0.400  A     0.39057      0.700
The below input files were used with Gromacs 4.0.4 in double precision.
[ moleculetype ]
TWOA 1

[ atoms ]
;   nr   type  resnr residue  atom   cgnr     charge       mass
     1     AA    1      TWOA   AA      0       0.4000      10.0000
     2     BB    1      TWOA   BB      0      -0.4000      10.0000

[ system ]
; name
Two atoms

[ molecules ]
; name   number
TWOA      1

grompp.mdp
-----------------
integrator               = md
dt                       = 0.002
nsteps                   = 1
ns_type                  = grid
rlist                    = 0.9
rcoulomb                 = 0.9
rvdw                     = 0.9
tcoupl                   = no
pbc                      = xyz
coulombtype              = Ewald
;fourierspacing           = 0
ewald_rtol               = 1e-5
fourier_nx                = 0
fourier_ny                = 0
fourier_nz                = 0
pcoupl                   = no


Thank you,

Steve Fiedler