[gmx-users] Good method to turn off electrostatics

[Andrew DeYoung]

Hi, 

I am running NVE simulations of a system of 1000 water molecules.  I am
attempting to run analogous simulations on Gromacs 4.5.4 and also on another
MD program.  Once finished, I will compare the results of the two programs.
It was suggested to me that to compare the results of the two programs in a
step-by-step manner, as a first step I should run the simulations with the
electrostatics (both short-range and long-range interactions) turned off.
In other words, run the simulations considering only Lennard-Jones
interactions.  My question is, what is the best way to do this in Gromacs?

A (very much more experienced) colleague suggested that I "turn off
electrostatics" by simply setting all the charges to zero.  However, when I
look into the topology and structure files, I see some surprising (to me)
things -- the charges do not appear to be specified in the .top and .gro
files.  

To back up, I will first try to explain my system.  I have a PDB file
specifying 1000 water molecules that someone gave me.  This PDB file looks
like the following:

TITLE     water
REMARK    THIS IS A SIMULATION BOX
CRYST1   47.100   47.100   47.100  90.00  90.00  90.00 P 1           1
MODEL        1
ATOM      1  OW  SOL     1     -21.177 -21.195 -21.179  1.00  0.00
ATOM      2  HW1 SOL     1     -21.671 -21.195 -22.083  1.00  0.00
ATOM      3  HW2 SOL     1     -21.863 -21.195 -20.412  1.00  0.00
ATOM      4  OW  SOL     2     -16.467 -21.195 -21.179  1.00  0.00
ATOM      5  HW1 SOL     2     -16.961 -21.195 -22.083  1.00  0.00
ATOM      6  HW2 SOL     2     -17.153 -21.195 -20.412  1.00  0.00
...

And so on.  According to the PDB file format specification
(http://www.wwpdb.org/format/sect9.html), the last two columns in my file
correspond to the occupancies and temperature factors, respectively, of my
atom.  So, evidently, partial charges on the oxygen and hydrogens are not
specified in this PDB file.  

Next, I execute the command: 

pdb2gmx -f water.pdb -p water.top -i water.itp -o water_processed.gro

And, when prompted, I choose the OPLS-AA force field and the SPC/E water
model.  Apparently, Gromacs recognizes from my PDB file that my system
contains water (and only water).  Alternatively, I could have created a
topology file manually, without pdb2gmx, by simply creating a .top file
containing:
------
#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/spce.itp"

[ system ]
water

[ molecules ]
SOL   1000
------
and then simply feed grompp the topology file and the original PDB file when
preparing for a dynamics run (pdb2gmx, I understand, is meant for chains
rather than multiple molecules).  But, in any case, I took the long road and
used pdb2gmx to create the .top topology and a .gro coordinate file.  

What surprises me (a novice in Gromacs and MD in general), though, is that
it seems that neither the .top file nor the .gro file contains any reference
to the charges on the oxygens and hydrogens in this water-only system.  For
example, my .top file contains:
------
; Include forcefield parameters
#include "oplsaa.ff/forcefield.itp"

; Include water topology
#include "oplsaa.ff/spce.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

; Include topology for ions
#include "oplsaa.ff/ions.itp"

[ system ]
; Name
water

[ molecules ]
; Compound        #mols
SOL              1000
------

And, my .gro file contains, for example: 
------
water
 3000
    1HOH     OW    1  -2.118  -2.119  -2.118
    1HOH    HW1    2  -2.167  -2.119  -2.208
    1HOH    HW2    3  -2.186  -2.119  -2.041
    2HOH     OW    4  -1.647  -2.119  -2.118
    2HOH    HW1    5  -1.696  -2.119  -2.208
    2HOH    HW2    6  -1.715  -2.119  -2.041
...
   4.71000   4.71000   4.71000
------

So, as far as I can tell, neither my .top nor my .gro file specify the
charges on the atoms.  I guess that this means that the charges on the atoms
are only specified in the spce.itp file.  When I find and open the spce.itp
file, I see: 
------
[ moleculetype ]
; molname       nrexcl
SOL             2

[ atoms ]
;   nr   type  resnr residue  atom   cgnr     charge       mass
     1  opls_116   1    SOL     OW      1      -0.8476
     2  opls_117   1    SOL    HW1      1       0.4238
     3  opls_117   1    SOL    HW2      1       0.4238

#ifndef FLEXIBLE
[ settles ]
; OW    funct   doh     dhh
1       1       0.1     0.16330

[ exclusions ]
1       2       3
2       1       3
3       1       2
#else
[ bonds ]
; i     j       funct   length  force.c.
1       2       1       0.1     345000  0.1     345000
1       3       1       0.1     345000  0.1     345000

[ angles ]
; i     j       k       funct   angle   force.c.
2       1       3       1       109.47  383     109.47  383
#endif
------

So, there the charges are specified: -0.8476 on the oxygen and 0.4238 on
each hydrogen.  Coming back to my original goal--to "turn off
electrostatics" by zeroing the charges--I could in principle just set the
charges in this .itp file to 0.  However, on our cluster, this file appears
to be READ ONLY.  So to change the charges to 0, I need to see if I can get
permission from our administrator. 

Based on my story above, my questions are:

1.  Is there a better way to "turn off electrostatics"--without modifying
the spce.itp file?

2.  Why does it seem that charges are not specified in the .top and .gro
files?  Is this just because I am boringly simulating only a solvent, and
the topology file tells Gromacs to look in spce.itp when simulating the
water molecules?

Thank you very much for your time!

Andrew DeYoung
Carnegie Mellon University