[gmx-users] Force-field checking options

[Mark Abraham]

On 3/08/2011 3:34 PM, mcgrath wrote:
> 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?

Uh, they're available wherever you were getting energy break-downs from 
(e.g. in the .log file, or via g_energy on the .edr file). You can't get 
a break-down for each interaction, however. Simplify your system to 
probe things here. Do zero-step MD (not EM) without constraints, to 
evaluate the energy of a single conformation, and compare that with your 
other software. Complex things are complex to compare. :-) Reduce the 
complexity.

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

...and the magnitude of its energy contribution is available.

Mark

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