[gmx-developers] Bug in Ewald Code?

[Dave Rogers]

Developers,

Problem:
  The ewald potential of the input configuration changes from
sr=0.000000   recip=-0.094927 to sr=-460.603673 recip=-330.965952 for
a single water molecule (all system files included below) when the
bond type ([bonds] funct) is changed from 1 (default) to 6 (harmonic)
for both of the bonds in Voth's water model!  For now I am unable to
reproduce his numbers for that water model as he used a harmonic bond
potential.

Thanks for looking into this,
~ David Rogers

  Version: gromacs 3.3.1 compiled with options (which uses fftw3 by default):
./configure --prefix=/usr/local --enable-mpi --disable-float
--with-external-blas --with-external-lapack --program-suffix=
## --------- ##
## Platform. ##
## --------- ##

hostname = yogen.cros.uc.edu
uname -m = x86_64
uname -r = 2.6.17-prepYOGEN
uname -s = Linux
uname -v = #2 SMP Wed Aug 23 13:04:18 EDT 2006

Files:

<FILE: flex_voth.itp>
[ moleculetype ]
; molname       nrexcl
SOL             2

[ atoms ]
;   nr   type  resnr residue  atom   cgnr     charge       mass
     1     OW      1    SOL     OW      1      -0.82
     2     HW      1    SOL    HW1      1       0.41
     3     HW      1    SOL    HW2      1       0.41

[ bonds ]
; i     j       funct   length  force.c.
1       2       1       0.1012  443153.3808
1       3       1       0.1012  443153.3808

[ angles ]
; i     j       k       funct   angle   force.c.
2       1       3       1       113.24  317.5656
</FILE>

<FILE: topol.top>
#include "ffgmx.itp"
#include "flex_voth.itp"

[ system ]
Pure Water

[ molecules ]
SOL     1
</FILE>

<FILE: gro.mdp>
; MDP for computing Voth's flex-spc model properties with Gromacs.

; VARIOUS PREPROCESSING OPTIONS
title                    = Yo
cpp                      = /usr/bin/cpp
include                  =
define                   =

; RUN CONTROL PARAMETERS
integrator               = md
dt                       = 0.001
nsteps                   = 1
tinit                    = 0

comm-mode                = Linear
nstcomm                  = 1

nstxout                  = 1
nstvout                  = 1
nstfout                  = 1

nstlog                   = 1
nstenergy                = 1

nstlist                  = 1
ns_type                  = simple
pbc                      = xyz
domain-decomposition     = no

coulombtype              = PME
rcoulomb                 = 0.9
vdw-type                 = Cut-off
rvdw                     = 0.9
rlist                    = 0.9

pme_order                = 6
ewald_rtol               = 3.95403e-8
optimize_fft             = no
; FFT grid size, when a value is 0 fourierspacing will be used
fourier_nx               = 32
fourier_ny               = 32
fourier_nz               = 32
</FILE>

<FILE: wat.pdb>
TITLE     Single, lonely water molecule
REMARK    THIS IS A SIMULATION BOX
CRYST1   19.000   19.000   19.000  90.00  90.00  90.00 P 1           1
MODEL        1
ATOM      1  OW  SOL     1       2.300   6.280   1.130  1.00  0.00
ATOM      2  HW1 SOL     1       1.370   6.260   1.500  1.00  0.00
ATOM      3  HW2 SOL     1       2.310   5.890   0.210  1.00  0.00
ENDMDL
</FILE>

<FILE: run.sh>
#!/bin/sh

grompp -f gro.mdp -c wat.pdb -o gro.tpr
mdrun -s gro.tpr -g gro -e gro -o gro -c gro
cat <<END | g_energy -f gro.edr -o gro.en.xvg
Bond
Angle
LJ-(SR)
Coulomb-(SR)
Coul.-recip.
Potential
Kinetic-En.
Pressure-(bar)
Vir-XX
Vir-YY
Vir-ZZ

END

cat gro.en.xvg
</FILE>