[gmx-users] charmm36 gromacs/charmm simulation package gives different electrostatic energies with same charges (probably I am using charmm incorrectly)

[Christopher Neale]

Dear Gromacs users:

This is a hybrid gromacs/charmm question, but there is not such a mailing list and I am hoping this is a suitable place for this question.

I am trying to convert some charmm36 parameters (from a new paper) to gromacs format and for this I am doing single-point energy evaluations in both gromacs and charmm to ensure accuracy. After finding some differences, I reduced my test system to a 2-residue protein sequence with zwitterionic termini (not using any new parameters at all). I find that although coordinates (both via .pdb for equal precision) and charges (via gmx .top and charmm .psf) are the same, there is some difference in the electrostatic energy computed. For this 2-residue test system, the difference is only 0.4 kcal/mol, but it adds up and the difference gets up to 28 kcal/mol for a few hundred protein residues (again, not anything other than standard charmm36 for this test), while all of the other energy terms are within the noise of floating point precision (maybe 0.001 kcal/mol for the ~200 residue system).

For the gromacs evaluation:
gmx_d grompp -f md.mdp -p topol.top -c conf.pdb -o MD.tpr -maxwarn 
gmx_d mdrun -s MD.tpr -rerun conf.pdb -nt 1
## the grompp maxwarn -1 is to avoid a warning based on oscillation period, which I don't care about for this evaluation.

For the charmm evaluation:
(1) I had to modify the .pdb from gmx editconf to remove the weird gromacs wrapping of atom names:
cat this.pdb |awk '{v=substr($0,13,1); $0=substr($0,1,12)substr($0,14,3)""v""substr($0,17,length($0)) ; print $0}' > new.pdb
(2) I had to modify H1, H2, and H3, to HT1, HT2, and HT3 in the .pdb file
(3) Then I ran charmm like this:

* single point energy evaluation
*

bomblev -1

OPEN UNIT 1 CARD READ NAME "toppar/top_all36_prot.rtf"
READ RTF CARD UNIT 1 
CLOSE UNIT 1

OPEN UNIT 1 CARD READ NAME "toppar/par_all36_prot.prm"
READ PARA CARD flex UNIT 1
CLOSE UNIT 1

open read card unit 1 name "new.pdb"
read sequence pdb unit 1

generate PROT setup warn first NTER last CTER

rewind unit 1
read coor pdb unit 1 
close unit 1

IC PARAmeters

PRINt COOR SELEct .NOT. INITialized END
 
OPEN UNIT 1 WRITe CARD NAME new.psf
WRITe PSF CARD UNIT 1
* PSF File
*
close unit 1

OPEN UNIT 1 WRITe CARD NAME new.pdb
WRITe COOR PDB UNIT 1
*PDB with Coordinates
*
close unit 1
 
ENERgy CTONNB 980 CTOFNB 990 CUTNB 1000 ATOM 

stop

########################################################

gromacs version is 5.1.2
charmm version is c40b1 compiled as gnu_xlarge

############## gromacs .mdp file (for generating .tpr for mdrun-rerun):

integrator = sd
dt = 0.002
tinit = 0

constraints = none 
cutoff-scheme = Verlet
vdwtype = cutoff
rlist = 10
rvdw = 10
coulombtype = cutoff
rcoulomb = 10
DispCorr = no

ns_type = grid
nstcomm = 100

nstxout = 0
nstvout = 0
nstfout = 0
nstxtcout = 10000
nstenergy = 10000
nstlist = 10
nstlog=0 ; reduce log file size

ewald-rtol = 1e-5
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4


######### gromacs energies from the .log file are:

   Energies (kJ/mol)
           Bond            U-B    Proper Dih.  Improper Dih.          LJ-14
    2.38966e+02    1.14273e+02    5.62574e+01    1.00406e+02    2.87319e+01
     Coulomb-14        LJ (SR)   Coulomb (SR)      Potential    Kinetic En.
    6.15434e+02   -6.87961e+00   -6.89232e+02    4.57957e+02    6.78967e+01
   Total Energy    Temperature Pressure (bar)
    5.25853e+02    1.64972e+02    2.80686e-02


######### charmm energies are:

ENER ENR:  Eval#     ENERgy      Delta-E         GRMS
ENER INTERN:          BONDs       ANGLes       UREY-b    DIHEdrals    IMPRopers
ENER EXTERN:        VDWaals         ELEC       HBONds          ASP         USER
 ----------       ---------    ---------    ---------    ---------    ---------
ENER>        0    109.89151      0.00000     44.82315
ENER INTERN>       57.11421     19.53553      7.77627     13.44585     23.99764
ENER EXTERN>        5.22280    -17.20080      0.00000      0.00000      0.00000
 ----------       ---------    ---------    ---------    ---------    ---------


####################

evaluation of electrostatic difference is:

gromacs:
6.15434e+02 -6.89232e+02 = -73.798 kJ/mol

charmm:
-17.20080 kcal.mol = -71.96814720 kJ/mol




Thank you very much for your help,
Chris.