[gmx-users] Re: Problem of Kinetic Energy
Size Zheng
alexzheng42 at gmail.com
Wed Aug 3 19:33:22 CEST 2011
Hi, Mark,
Thank you for your response. I see, you meant g_traj with -ekt and -ekr options just operates on group which is investigated. Then is there any other tool to analyze the rotational and translational kinetic energy of the system or some groups on atoms or molecules?
Many thanks
Size Zheng
On Aug 3, 2011, at 3:00 AM, gmx-users-request at gromacs.org wrote:
> Send gmx-users mailing list submissions to
> gmx-users at gromacs.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> or, via email, send a message with subject or body 'help' to
> gmx-users-request at gromacs.org
>
> You can reach the person managing the list at
> gmx-users-owner at gromacs.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of gmx-users digest..."
>
>
> Today's Topics:
>
> 1. Re: Force-field checking options (Mark Abraham)
> 2. Re: Problem of Kinetic Energy (Mark Abraham)
> 3. Re: Force-field checking options (Mark Abraham)
> 4. Re: Force-field checking options (mcgrath)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Wed, 03 Aug 2011 15:45:13 +1000
> From: Mark Abraham <Mark.Abraham at anu.edu.au>
> Subject: Re: [gmx-users] Force-field checking options
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID: <4E38E069.8010809 at anu.edu.au>
> Content-Type: text/plain; charset=UTF-8; format=flowed
>
> 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
>>>>>
>
>
>
> ------------------------------
>
> Message: 2
> Date: Wed, 03 Aug 2011 15:54:17 +1000
> From: Mark Abraham <Mark.Abraham at anu.edu.au>
> Subject: Re: [gmx-users] Problem of Kinetic Energy
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID: <4E38E289.8030405 at anu.edu.au>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
> On 3/08/2011 3:25 PM, Size Zheng wrote:
>> Hi, there,
>>
>> I simulate a simple system of pure water with rigid TIP4P model in NVT at 273K.
>>
>> Why the total kinetic energy, reported from g_energy, is not equal to
>> the sum of the rotational kinetic energy and the translational kinetic
>> energy, which are obtained from g_traj with options of -ekt and -ekr.
>
> Read g_traj -h. The -ekt and -ekr options do not operate on atoms.
>
>>
>> Even though I normalize the total kinetic energy (reported from g_energy)
>> with molecular number, its value is still not equal to the sum of the
>> rotational and translational kinetic energy ,obtained from g_traj
>> with options of -ekt and -ekr.
>
> A stationary swarm of bees has high KE, and yet the swarm has low
> rotational and translational KE.
>
> Mark
>
>
> ------------------------------
>
> Message: 3
> Date: Wed, 03 Aug 2011 15:55:29 +1000
> From: Mark Abraham <Mark.Abraham at anu.edu.au>
> Subject: Re: [gmx-users] Force-field checking options
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID: <4E38E2D1.6010107 at anu.edu.au>
> Content-Type: text/plain; charset=UTF-8; format=flowed
>
> 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
>>>>>>
>>
>
>
>
> ------------------------------
>
> Message: 4
> Date: Wed, 03 Aug 2011 18:32:15 +0900
> From: mcgrath <mcgrath at theory.biophys.kyoto-u.ac.jp>
> Subject: Re: [gmx-users] Force-field checking options
> To: gmx-users at gromacs.org
> Message-ID: <1312363935.22887.53.camel at mcgrath-TP55>
> Content-Type: text/plain; charset="UTF-8"
>
> 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.
>
> 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
>>>>>>
>>
>
>
>
> ------------------------------
>
> --
> gmx-users mailing list
> gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>
> End of gmx-users Digest, Vol 88, Issue 23
> *****************************************
More information about the gromacs.org_gmx-users
mailing list