[gmx-users] Ewald summation: n=0
Steve Fiedler
fiedler at umich.edu
Tue Apr 21 20:08:31 CEST 2009
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
More information about the gromacs.org_gmx-users
mailing list