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

[Qiong Zhang]

			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,

  

  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



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