[gmx-users] Force-field checking options

[Justin A. Lemkul]

mcgrath wrote:
> Hi Mark.
> 
>> 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.
> 
> I must be missing something really obvious, then, because this is where
> I'm getting the energies from in the log file.
> 
>    Energies (kJ/mol)
>             U-B    Proper Dih.  Improper Dih.          LJ-14
> Coulomb-14
>     2.42228e+02    1.23963e+02    6.15742e-02    1.43130e+02   -5.01705e
> +03
>         LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total
> Energy
>     8.30424e+01   -4.63865e+02   -4.88849e+03    1.05098e+02   -4.78339e
> +03
>     Temperature Pressure (bar)   Constr. rmsd
>     3.00959e+02   -7.08408e+02    6.28488e-06
> 
> I don't see any bend energies there.  No bonds, either, but that's fine,
> because I was constraining them (constraints all-bonds).  Removing that
> constraint gives a bond energy in that section, but still no bend
> energy.
> 

All angles in CHARMM27 use Urey-Bradley type angles, so the angle bending 
energies are the U-B term.

-Justin

> Agree completely about simplification.  The first system I quoted for
> Justin was for two molecules (ATP + Mg), which is the smallest system
> that I see the problem for (as I mentioned, using 29 or 25000 waters
> gave a difference of only 1%, which is explainable by slight rounding of
> the coordinates...at least, possibly).  He then asked for the breakdown
> for the big system, so I gave him that, too.
> 
>> ...and the magnitude of its energy contribution is available.
> 
> Also agreed; that's how I concluded that its contribution was only a
> small part, and using only ATP there appears to be no CMAP (see the
> above energies).  But, even on this small system, I haven't been able to
> hunt down the differences.  The parameter files have the same values for
> the nonbonded parameters (converting between the GROMACS and CHARMM
> format) and charges, so I'm not seeing any simple solution.  Unless
> someone has an idea to save me, I guess I'm stuck with looking at
> individual contributions for this system.
> 
> Sigh.  Some days this field isn't so much fun.  :)  Especially since the
> bug will probably be something stupid I did in setting it all up.
> 
> Thanks!
> 
>                                      Cheers, Matt
> 
> 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
>>>>>>
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================