[gmx-users] gromacs energies very different from other MD engine for the very same system and conditions

[chris.neale at utoronto.ca]

Dear Alan:

I don't have a solution for you, but I do have a simpler test case: a  
single spc water molecule. I have not addressed your noted energy  
difference between programs because I am not familiar with them  
(although you should try gromacs compiled in double precision).

I can reproduce the non-zero kinetic energy with a single spc water,
as long as I add define=-DFLEXIBLE to the .mdp file. The output .gro
contains all zeroes for velocities (as long as the input .gro does not  
contain any velocities!):

title
      3
      1SOL     OW    1   0.230   0.628   0.113  0.0000  0.0000  0.0000
      1SOL    HW1    2   0.137   0.626   0.150  0.0000  0.0000  0.0000
      1SOL    HW2    3   0.231   0.589   0.021  0.0000  0.0000  0.0000
    10.00000  10.00000  10.00000

but the run indicates a non-zero kinetic energy (log and .edr file, gromacs
4.0.5 and 3.3.3).

The reported kinetic energy is small though, 1.16470e-03, with a
reported temperature of 9.33871e-02.

I added unconstrainted_start=yes to your .mdp file but that did not help.

I tried in double precision, and I get the same reported kinetic  
energy (1.16473e-03)

Chris.

-- original message --

Hi there,

I have a very simple case of a tri alanine (AAA). You can use tleap from
ambertools to create such peptide and save the pdb and amber's parameters.

Then I do a single point calculation to get the potential energy. I am using
mdin (amber input file) like:

cat << EOF >| mdin
Single point
&cntrl
imin=0, maxcyc=0,
ntmin=2,
ntb=0,
igb=0,
cut=999,
/
EOF

Then I test the same input parameters with Namd2, using namd.conf (Amber is
not freely available, but AmberTools and Namd are):

cat << EOF >| AAAamb_namd.conf
outputEnergies 1  # Energy output frequency
DCDfreq        1  # Trajectory file frequency
timestep       2  # in unit of fs
temperature    300  # Initial temp for velocity assignment
cutoff         999
switching      off  # Turn off the switching functions
PME            off  # Use PME for electrostatic calculation
amber          on  # Specify this is AMBER force field
parmfile       AAAamb.prmtop  # Input PARM file
ambercoor      AAAamb.inpcrd  # Input coordinate file
outputname     AAAamb  # Prefix of output files
exclude        scaled1-4
1-4scaling     0.833333  # =1/1.2, default is 1.0
minimize       0
EOF

I have only a 0.001 kcal/mol absolute diff (or relative 0.0037%) for the
Elec term, everything else is absolutely the same.

Then comes gromacs. I convert my prmtop and inpcrd to gromacs top and gro
with either acpype or amb2gmx and then doing a single point with file:

cat << EOF >| SPE.mdp
integrator               = md
nsteps                   = 0
dt                       = 0.001
constraints              = none
emtol                    = 10.0
emstep                   = 0.01
nstcomm                  = 1
ns_type                  = simple
nstlist                  = 0
rlist                    = 0
rcoulomb                 = 0
rvdw                     = 0
Tcoupl                   = no
Pcoupl                   = no
gen_vel                  = no
nstxout                  = 1
pbc                      = no
nstlog = 1
nstenergy = 1
nstvout = 1
nstfout = 1
nstxtcout = 1
comm_mode = ANGULAR
EOF

I got something like:
     Energies (kJ/mol)
             Bond          Angle Ryckaert-Bell.          LJ-14     Coulomb-14
      1.77662e+02    3.19281e+01    4.55707e+01    2.48926e+01    9.40060e+02
          LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total Energy
     -9.21735e+00   -1.05005e+03    1.60848e+02    8.84934e+00    1.69697e+02
      Temperature Pressure (bar)
      2.28887e+01    0.00000e+00

First question, why I have temperature and kinetic energy? Remember, it's a
single point energy calculation, without any sort of coupling, 0 velocities
and and gen_vel = no.

Secondly, converting gromacs energies (kJ/mol) by dividing the terms by
4.184 and they don't match amber/namd results. See relative diffs
(abs(g-a)/max(abs(a),abs(g))):

Bond:    5.87%
Angle:   3.49%
Dihe:    0.14%
vdW:     0.52%
Elec:    0.17%
Pot.Tot: 5.90%

I am very surprised to find this difference, in special for Bond since it's
for this term essentially the same equation and parameters. Any ideas of
what could be missing here?

Many thanks,

Alan