[gmx-users] Force-field checking options

[Justin A. Lemkul]

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

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