[gmx-users] Force-field checking options

[mcgrath]

Hi Justin, thanks for the response.  For the big system, the dihederal
is still fine, while now the improper, U-B, and nonbonded (CP2K groups
all nonbonded, 1-4, and real space Coulomb into the same term) are off
by a factor of 2-3.  Is there any way to print the energy due to angles
in GROMACS?  The nonbonded interactions dominate (for GROMACS, the short
range Coulomb term is an order of magnitude bigger than everything
else).  I guess it's possible my Ewald parameters are off, though that
wouldn't affect the SR term.  My box is 85x85x105 and I use a grid of
88x88x108, along with a 4th order interpolation for both simulations,
along with an ewald_rtol of 1e-5.  Playing with these doesn't seem to
change the numbers much, nor does switching between PME, SPME, or
regular Ewald (for CP2K).  rlist=rvdw=rcoulomb=1.2.  Is there a
parameter that I'm missing?  I've tried playing around with vdwtype and
coulombtype (I believe CP2K truncates and shifts nonbond and Coulombic
realspace potentials), but nothing had a significant effect.

Thanks for the gmxdump tip.  I'll start looking at the parameters for
the smaller system, and hopefully I can hunt down where the differences
are popping up. CP2K doesn't have CMAP implemented, but that only a
small part according to GROMACS.

                          Cheers, Matt

> > -------- Forwarded Message --------
> > From: Justin A. Lemkul <jalemkul at vt.edu>
> > Reply-to: jalemkul at vt.edu, Discussion list for GROMACS users
> > <gmx-users at gromacs.org>
> > To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> > Subject: Re: [gmx-users] Force-field checking options
> > Date: Wed, 03 Aug 2011 00:12:35 -0400
> > 
> > 
> > mcgrath wrote:
> > > Hi.  I'm fairly new to GROMACS, and I've been using it to run some
> > > classical MD simulations.  In order to be sure that I'm using the
> > > correct force field (I had to add a molecule to CHARMM27), I'm comparing
> > > it to another simulation code that I know well, CP2K (using GROMACS
> > > 4.5.4 and a recent CVS version of CP2K).  Sadly, they are giving me
> > > energy differences of about a factor of 2 for a 75000 atom protein+water
> > > system.  As far as I can tell, I'm using the same PME parameters, and
> > > there's not a big change in energy when I change those, anyway.  I've
> > > been able to confirm that it's not a global problem (computing the
> > > energy of only the 25,000 TIP3P waters give a result to within 1%, which
> > > is not perfect, but better...I seem the same 1% difference if I only use
> > > 29 waters).  For a smaller system (using the molecule I added), the
> > > total energy is incorrect, but the torsion and improper torsions are
> > > good to within 1%.  So it looks like my topology is perhaps being
> > > incorrectly specified, or the parameters for them.
> > 
> > Which energy terms show the discrepancy in the case of the twofold difference?
> > 
> > > 
> > > What I would like to do is get GROMACS to print out all the charges (the
> > > electrostatics/nonbonded are also different) that it is actually using,
> > > as well as the force field parameters being used.  By using mdrun -v
> > > -debug I get some of that, but lines like
> > > 
> > 
> > You don't need mdrun -debug for this.  The topology is static, so run gmxdump on 
> > your .tpr file.  This will map all parameters to each atom.
> > 
> > -Justin
> > 
> > > c6= 1.48497790e-03, c12= 6.58807176e-06
> > > c6= 3.78914556e-04, c12= 4.08982345e-07
> > > c6= 2.02168059e-03, c12= 4.42785949e-06
> > > c6= 1.71635114e-03, c12= 4.54480232e-06
> > > c6= 2.60732602e-03, c12= 6.42257237e-06
> > > c6= 3.75109637e-04, c12= 2.77185705e-07
> > > c6= 1.98187144e-03, c12= 5.24788538e-06
> > > c6= 6.18919948e-05, c12= 7.54611396e-09
> > > c6= 1.73354731e-03, c12= 5.36550624e-06
> > > c6= 4.41942073e-04, c12= 4.93156790e-07
> > > c6= 6.79507386e-03, c12= 2.55060841e-05
> > > c6= 1.61714898e-03, c12= 3.18964840e-06
> > > c6= 2.48678902e-04, c12= 2.13336705e-07
> > > 
> > > are somewhat unclear.  They are obviously non-bonded terms, but
> > > corresponding to which atoms?  The same goes for the bond angles and
> > > torsions printed later.  The worst-case scenario is that I have to poke
> > > around in the source files, which I would to avoid so I'm hoping there
> > > is some documentation or more switches I can flip somewhere.  Thanks!
> > > 
> > >                    Cheers, Matt
> > > 
> > 
>