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

[Justin Lemkul]

On 8/7/16 2:21 AM, Christopher Neale wrote:
> 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.
>

You can send these questions directly to me; I maintain the port and did a lot
of the validation of it.

> 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.
>

If you send me all of your input files I will look into it.  Nothing obvious 
stands out, but I'll need to pull some more information, at least from the 
CHARMM side, to figure out what's going on.

-Justin


-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================