[gmx-users] Re:problem with interaction energy calculated by g_energy

[Mark Abraham]

On 9/03/2010 9:32 PM, Qiong Zhang wrote:
>
>
>
>   Hi gmx users,
>
> I found the big discrepancy between the interaction energy I got from my
> first approach and send approach should be ascribed to a bug reported here:
>
> http://www.mail-archive.com/gmx-users@gromacs.org/msg20963.html
>
> The gromacs I am using now is exactly gmx4.0.4. I also reran with a
> parallel version and the energies never changed during the rerun stage.

Well that's why we tell people to report their GROMACS version. :-) 
Using the latest version, and announcing what you are using can help you 
avoid wasting people's time :-)

> Still, the discrepancy in the energies between the second approach and
> the third approach is still puzzled to me. Which one is the correct way
> of calculating interaction energy?

Like I said last time, you can't do this with PME. The reciprocal-space 
calculation cannot be decomposed group-wise. Go read up on PME if you 
don't understand this.

Also last time I pointed out this was a non-problem, for such an 
interaction energy doesn't mean much of anything anyway, even if you 
calculate it with some other electrostatics model.

Mark

>   [gmx-users] Re:problem with interaction energy calculated by g_energy
>
> Qiong Zhang
> Tue, 09 Mar 2010 01:17:02 -0800
>
>
>
> Hi dear Mark,
>
>
>
> Please ignor my last mail replied to you. I made some mistake there.
>
>
>
> Yes, you are right that I am using PME. The cutoff for the real space and
> reciprocal space is 1.2nm.
>
>
>
> The molecules I am simulating are carbohydrates. And I am using Glycam06 Force
> Field.
>
>
>
> I tried there
> different ways to calculate the interaction energy:
>
>
>
> The first approach is analyzed by directly using g_energy, summing up Coul_SR
> and LJ_SR of two groups, since in the .mdp file I have defined in energygrps 1
> 2.
>
> The interaction energy between 1 and 2 (E 1_2) = E
> Coul_SR + E LJ_SR =-170.048+(-232.719)=-402.767 kJ/mol
>
>
>
> The second approach is
> using "mdrun -rerun" option with the exactly the same energygrps 1 2 defined
> in .mdp, the same traj.xtc and the same index. Weird enough, this time, I got
> interaction
> energy between 1 and 2
>   (E 1_2) = E Coul_SR + E LJ_SR
> = -91.5234 + (-238.712) = -330.235 kJ/mol, which is quite far from the
> previously -402.767 kJ/mol!!!! But this -330.235 kJ/mol is the exact sum of the
> contributions of subunits. The contributions of subunits are also calculated in
> this approach with rerun. So the discrepancy I reported in my first mail is
> solved.
>
>
>
> But what is the reason for the huge discrepancy between
> the interaction energy from the original run and the “rerun”?? I think they
> should be exactly the same.
>
>
>
> The third approach, in order to include the long range interaction, I've also
> tried "mdrun -rerun" option with three
> "reruns" carried out for molecule 1(1st), molecules 2 (2nd) and
> molecule 1 and 2 (3rd). The interaction energy for molecule 1 and 2 is now
> calculated by:
>
>
>
> [Coul(SR+recip)+LJ(SR+Disper. corr.)]_3rd - [Coul(SR+recip)+LJ(SR+Disper.
> corr.)]_2nd -
>   [Coul(SR+recip)+LJ(SR+Disper. corr.)]_1st
>
> =Delta(Coul_SR)+Delta(Coul_recip)+Delta(LJ_SR)+Delta(LJ_Disper.corr.)
>
> =(-128.73) + (-30.33) +( -252.021) + (-39.9) = -450.217 kJ/mol
>
>
>
> If we neglect the long-range interactions, namely, Delta(Coul_recip) and
> Delta(LJ_Disper.corr.),
> we got the interaction energy -128.73
> -252.021= -380.751 kJ/mol. We see here the long-range
> contribution is not negligible. However, this short range energy -380.751
> kJ/mol is neither close to the -330.235 kJ/mol nor -402.767 kJ/mol.
>
>
>
> So Now I am confused. Which approach should be really
> adopted in the calculation of interaction energy? And what approach do you use
> in such interaction energy calculations?
>
>
>
> Thank you very much!
>
>
>
> Qiong
>
>
>
>
>
> --- On Tue, 3/9/10, Qiong Zhang<qiongzhang... at yahoo.com>
> wrote:
>
>
>
> From: Qiong Zhang
>   <qiongzhang... at yahoo.com>
>
> Subject: Re:problem with interaction energy calculated by g_energy
>
> To: gmx-users at gromacs.org
>
> Date: Tuesday, March 9, 2010, 4:27 PM
>
>
>
>
>    Hi dear Mark,
>
>
>
>    Thanks very much for your reply.
>
>
>
>    Yes, you are right that I am using PME.
>
>
>
>    The molecules I am simulating are carbohydrates. And I am using Glycam06
>    Force Field.
>
>
>
>    The interaction energy I got previously is analyzed by directly using
>    g_energy, summing up Coul_SR and LJ_SR of two groups.
>
>
>
>    In order to include the long range interaction, I've also tried "mdrun
>    -rerun" option.    So three "reruns" were carried out for
>    molecule 1(1st), molecules 2 (2nd) and molecule 1 and 2 (3rd). This time, I
>    found the long range Coul_recip between molecule 1 and 2 is a quite positive
>    value. So when only Coul_SR is included, the
>   electrostatic interaction
>    between molecule 1 and molecules 2 is much more negative (>  100 kj/mol)
>    than that when both Coul_SR and Coul_recip are included. I guess, for such
>    carbohydrate molecules, long range Coul_recip can not be excluded.
>
>    Am I right here?
>
>
>
>    For the second summing up problem, I am still checking all the input file,
>    especially the index file.
>
>
>
>    Thank you very much!
>
>
>
>    Qiong
>
>
>
>    ----- Original Message -----
>
>    From: Qiong Zhang<qiongzhang... at yahoo.com>
>
>    Date: Monday, March 8, 2010 20:35
>
>    Subject: [gmx-users] problem with interaction energy calculated by g_energy
>
>    To: gmx-users at gromacs.org
>
>
>
>    -----------------------------------------------------------
>
>    |>  Dear gmx users,
>
>    >
>
>    >  I am studying the adsorption behavior of a molecule ( molecule 1) on a
>
> surface
>    (molecules 2). Based on the production run, I calculated the interaction
>    energy between molecule 1 and molecules 2 by g_energy.
>
>    >  Here comes the first question: Why only short range interactions between
>    1 and 2 are displayed, namely, Coul_SR and LJ_SR? So the interaction energy E
>    1_2 I calculated is just the sum of Coul_SR+LJ_SR. Will this bring about huge
>    errors?
>
>
>
>    Guessing wildly (since you've not told us the nature of your simulation
>    protocol) you're using PME, and so the long-range contributions cannot be
>    decomposed group-wise. This is probably a good thing - I'm not aware of any
>    force field that has been parameterized so that small chunks of atoms
>    interaction energies correlate to anything useful.
>
>
>
>    >  After this, I'd like to know the individual contributions of the
>    components of molecule 1    to the interaction energy between 1 and 2.
>   For
>    example, molecule 1 is composed of A, B, C and D resdues. So again, by
>    g_energy, I got interaction energy between A, B, C and D with 2,
>    respectively, denoted by E A_2, E B_2, E c_2 and E D_2.
>    Still, these interaction energies are the sum of
>
>    Coul_SR+LJ_SR.
>
>    >  Then comes the second question: Why the sum of E A_2, E B_2, E c_2 and E
>    D_2 does not equal to E 1_2? I found there was big difference between them,
>    sometimes as large as 50 kJ/mol.
>
>    >
>
>    >  Could anybody give me some hints or suggestions please?
>
>
>
>    They should add up. Check your index group definitions and use in the .mdp
>    file.
>
>
>
>    Mark
>
>