On 3/08/2011 9:20 PM, Justin A. Lemkul wrote:
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.
Different forcefields/programs use the label "Urey-Bradley" to refer to
either the sum of the harmonic angle potential and the harmonic 1-3
potential (e.g. CHARMM27 in GROMACS), or just the latter (e.g. CHARMM27
in CHARMM). So mileage does vary.
-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).
You didn't say how you were measuring the energies, but if it was MD
averages, then your averages are over different ensembles, even from the
same starting point with notionally the same algorithm.
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.
Probably. Getting the output from grompp -pp, and then choosing atom(s)
to have zero charge, and/or zero VDW parameters is a good way to
trouble-shoot, assuming your other software can be massaged to replicate
the effect. Before that, check that the energy of a single amino-acid
can be replicated.
Mark
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
--
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
