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


_______________________________________________
gmx-users mailing list    gmx-users@gromacs.org
http://www.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