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

Qiong Zhang qiongzhang928 at yahoo.com
Tue Mar 9 11:32:24 CET 2010

```

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.

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?

Thank you very much!

Qiong
[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,

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100309/170b584d/attachment.html>
```