[gmx-users] Nucleic Acid in vacuo

[Matthias Ernst]

Dear Gromacs users,

I am currently working on porting forcefields from Gromacs to another
program with the aim of simulating DNA. Therefore, I wrote a little
script which reads in a PDB file, evaluates the Lennard-Jones and the
Coulomb energy in vacuo with the forcefield data from Gromacs.  I tried
to implement the search for the neighbours by comparing inter-atomic
distances to the van-der-Waals-radii and a flexible list-based factor
weighting of the various contributions (where -1 is a default value, 0
is the atom itsself, 1 a bond to a next neighbour and so on, so my
factor "3" should be equivalent to the "1-4-interaction" in Gromacs).

To verify what I have done, I tried to do a calculation of a single "DA"
nucleic acid with Gromacs to compare the LJ and Coulomb energy.
Unfortunately, the results do not match and I do not understand why.

Please find the PDB code I used here http://pastebin.com/zB3sQHdZ and
the python code I used for evaluating the potentials here
http://pastebin.com/2v4nVdg3

The results I got with my code are 
-> Total Lennard-Jones Energy [kJ/mol]:   0.712752493372
-> Total 1-4-Lennard-Jones Energy [kJ/mol]:   37.1121824156
-> Total Coulomb Energy [kJ/mol]:   -4354.55160305
-> Total Coulomb 1-4 Energy [kJ/mol]:   -207.712332451

Using Gromacs, I got
   Energies (kJ/mol)
Angle    Proper Dih.  Improper Dih.          LJ-14     Coulomb-14
2.00975e+01    6.05485e+01    4.28869e-04    3.69509e+01   -2.06010e+02
LJ (SR)   Coulomb (SR)      Potential Pressure (bar)   Constr. rmsd
-1.45022e+00   -5.67518e+01   -1.46615e+02    0.00000e+00    1.80690e-04

The option in my mdp are
title           =       INITIAL_EM0
cpp             =       /lib/cpp
constraints     =       all-bonds
integrator      =       steep
nsteps          =       50
nstcomm         =       1
ns_type         =       grid
rlist           =       0.0
rcoulomb        =       0.0
rvdw            =       0.0

;
; Energy minimizing stuff
;
emtol           =       1500.0
emstep          =       0.01

Tcoupl          =       no
Pcoupl          =       no
gen_vel         =       no
pbc = no

I have not yet done a lot of calculations with Gromacs so maybe I did
something wrong. If so, please correct me.

What I wonder is that the 1-4-energies seem quite reasonable (LJ: 37.11
vs. 36.95, Coulomb -207.7 vs. -206.0) but the rest is rather different,
especially the sign of the LJ term and the order of magnitude of the
Coulomb term. Why is that? What does the "SR" in LJ and Coulomb in the
Gromacs output mean?

Any help would be thankfully apreciated.

Best regards,
Matthias Ernst