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

Ilja Khavrutskii ivkhavru at gmail.com
Tue Jul 27 17:06:50 CEST 2010


I think the bond terms differ simply due to the round off error. The gro
file is worse than pdb when it comes to round off in Cartesian coordinates.
You can confirm this is the case by taking your gro file and converting it
back to pdb using gromacs tools and then use the resulting pdb (with the
same round of as the gro file now) file from gromacs in NAMD. You should see
the same numbers for bonded terms then.

On Tue, Jul 27, 2010 at 10:56 AM, Berk Hess <gmx3 at hotmail.com> wrote:

>  Hi,
>
> The kinetic energy is non-zero because Gromacs uses a leap frog integrator
> and the kinetic energy
> at t=0 is the average of the kinetic energies at -dt/2 and +dt/2 (as
> described in the manual).
>
> Eric Sorin did comparisons between Amber and Gromacs (uses his Amber force
> field port in Gromacs)
> and found differences only of the size of roughly the machine precision. So
> I don't think there is an issue
> in Gromacs. But I have no experience with the route you took.
>
> Berk
>
> > Date: Tue, 27 Jul 2010 10:43:30 -0400
> > From: chris.neale at utoronto.ca
> > To: gmx-users at gromacs.org
> > Subject: [gmx-users] gromacs energies very different from other MD engine
> for the very same system and conditions
>
> >
> > 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
> >
> >
> >
> > --
> > gmx-users mailing list gmx-users at gromacs.org
> > http://lists.gromacs.org/mailman/listinfo/gmx-users
> > Please search the archive at http://www.gromacs.org/search before
> posting!
> > Please don't post (un)subscribe requests to the list. Use the
> > www interface or send it to gmx-users-request at gromacs.org.
> > Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
> --
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100727/924f29d4/attachment.html>


More information about the gromacs.org_gmx-users mailing list