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

[Mark Abraham]

On 10/03/2010 7:56 AM, Qiong Zhang wrote:
> Hi dear Mark,
>
> Thank you very much for your reply!
>
> Yes, you are right that I should have stated the gromacs version in my first
> mail. I am sorry that I did not notice this issue. I will pay attention
> to this next time.
>
> As for the electrostatic interaction energy in the long range, I
> am afraid that I have some different opinion which I am not sure
> if it is correct or not. I think for some systems with strong
> electrostatic interaction, for example, the interaction between a
> Rutile (TiO2) surface and a protein, it seems that the electrostatic
> interaction energy in the long range plays a very important role in
> the total interaction energy as one of my colleagues shows. In such
> cases, I think the electrostatic interaction energy in the long
> range can not be neglected. What is your
>   opinion please?

Important, yes - you need long-range electrostatics to sample the right 
ensemble. Numerically meaningful when extracted from the whole 
condensed-phase ensemble, no. If it's low, then the total energy has 
sloshed into other degrees of freedom - so what? This is not gas-phase 
ab initio quantum chemistry at 0K, where internal energy correlates with 
something useful, because there are no other energetic degrees of 
freedom. The frequency of occurrence of a region of structure space in a 
converged trajectory can tell you something, i.e. relative free energies.

> And I think I understand now"the reciprocal-space
> calculation cannot be decomposed group-wise."  Maybe a better way
> to overcome this is using the formula:
>
> E_interact=E_tot(1-2)-E_tot(1)-E_tot(2)
>
> Do you agree with this?

No. The only term with long-range contributions is the reciprocal-space 
term and it cannot be decomposed. There is no way around this.

If you can find a published article explaining the usefulness of the 
analysis you're trying to do, they'll have used a forcefield and 
electrostatics model that are consistent with doing it. You should copy 
their method, in that case. I've given my advice three times, and am 
going to desist :-)

Mark

> I am highly appreciative for all your help!
>
> Qiong
>
>
>
> 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
>
>