[gmx-users] Force-field checking options

[Mark Abraham]

On 3/08/2011 3:45 PM, Mark Abraham wrote:
> 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.

...and trivial to turn off in pdb2gmx.

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