incorrect evaluation of system energy

What steps will reproduce the problem?

1. Edit or run the test job ambersim.py to add the following instructions after 
the script has completed
2. system.energy()
   system.energies()
   system.mustNowRecalculateFromScratch()
   system.energy()
   system.energies()

3. What is the expected output? What do you see instead?

The total energy recomputed from scratch should be equal to the total energy of 
the system recorded for the last computed configuration. 
Instead systematic differences are observed, the actual differences will depend 
on the trajectory, but to illustrate. 

In [2]: system.energy()
Out[2]: -17623.6 kcal mol-1
In [3]: system.energies()
Out[3]: 
{ E_{protein:solvent}^{CLJ} == -265.784, E_{protein:solvent}^{LJ} == -97.1268, 
E_{protein:solvent}^{coulomb} == -168.658, E_{protein_intraclj}^{CLJ} == 
-4350.77
  E_{protein_intraclj}^{LJ} == -227.492, E_{protein_intraclj}^{coulomb} == -4123.27, E_{protein_intraff}^{angle} == 4553.08, E_{protein_intraff}^{bend-bend} == 0
  E_{protein_intraff}^{bond} == 908.359, E_{protein_intraff}^{dihedral} == 2608.51, E_{protein_intraff}^{improper} == 281.156, E_{protein_intraff}^{internal} == 8351.1
  E_{protein_intraff}^{stretch-bend-torsion} == 0, E_{protein_intraff}^{stretch-bend} == 0, E_{protein_intraff}^{stretch-stretch} == 0, E_{protein_intraff}^{urey-bradley} == 0
  E_{solute:protein}^{CLJ} == -20.6086, E_{solute:protein}^{LJ} == -7.37734, E_{solute:protein}^{coulomb} == -13.2313, E_{solute:solvent}^{CLJ} == -4.75784
  E_{solute:solvent}^{LJ} == -3.81741, E_{solute:solvent}^{coulomb} == -0.940426, E_{solute_intraclj}^{CLJ} == -44.2914, E_{solute_intraclj}^{LJ} == 44.3799
  E_{solute_intraclj}^{coulomb} == -88.6713, E_{solute_intraff}^{angle} == 41.2044, E_{solute_intraff}^{bend-bend} == 0, E_{solute_intraff}^{bond} == 59.483
  E_{solute_intraff}^{dihedral} == 17.6433, E_{solute_intraff}^{improper} == 0.604145, E_{solute_intraff}^{internal} == 118.935, E_{solute_intraff}^{stretch-bend-torsion} == 0
  E_{solute_intraff}^{stretch-bend} == 0, E_{solute_intraff}^{stretch-stretch} == 0, E_{solute_intraff}^{urey-bradley} == 0, E_{solvent:solvent}^{CLJ} == -21407.4
  E_{solvent:solvent}^{LJ} == 39446.8, E_{solvent:solvent}^{coulomb} == -60854.3, E_{total} == -17623.6 }

In [4]: system.mustNowRecalculateFromScratch()

In [5]: system.energy()
Out[5]: -17466.7 kcal mol-1
In [8]: system.energies()
Out[8]: 
{ E_{protein:solvent}^{CLJ} == -264.851, E_{protein:solvent}^{LJ} == -97.4211, 
E_{protein:solvent}^{coulomb} == -167.429, E_{protein_intraclj}^{CLJ} == 
-4350.77
  E_{protein_intraclj}^{LJ} == -227.492, E_{protein_intraclj}^{coulomb} == -4123.27, E_{protein_intraff}^{angle} == 4553.08, E_{protein_intraff}^{bend-bend} == 0
  E_{protein_intraff}^{bond} == 908.359, E_{protein_intraff}^{dihedral} == 2608.51, E_{protein_intraff}^{improper} == 281.156, E_{protein_intraff}^{internal} == 8351.1
  E_{protein_intraff}^{stretch-bend-torsion} == 0, E_{protein_intraff}^{stretch-bend} == 0, E_{protein_intraff}^{stretch-stretch} == 0, E_{protein_intraff}^{urey-bradley} == 0
  E_{solute:protein}^{CLJ} == 135.372, E_{solute:protein}^{LJ} == 149.034, E_{solute:protein}^{coulomb} == -13.662, E_{solute:solvent}^{CLJ} == -4.75784
  E_{solute:solvent}^{LJ} == -3.81741, E_{solute:solvent}^{coulomb} == -0.940426, E_{solute_intraclj}^{CLJ} == -44.2914, E_{solute_intraclj}^{LJ} == 44.3799
  E_{solute_intraclj}^{coulomb} == -88.6713, E_{solute_intraff}^{angle} == 41.2044, E_{solute_intraff}^{bend-bend} == 0, E_{solute_intraff}^{bond} == 59.483
  E_{solute_intraff}^{dihedral} == 17.6433, E_{solute_intraff}^{improper} == 0.604145, E_{solute_intraff}^{internal} == 118.935, E_{solute_intraff}^{stretch-bend-torsion} == 0
  E_{solute_intraff}^{stretch-bend} == 0, E_{solute_intraff}^{stretch-stretch} == 0, E_{solute_intraff}^{urey-bradley} == 0, E_{solvent:solvent}^{CLJ} == -21407.4
  E_{solvent:solvent}^{LJ} == 39446.8, E_{solvent:solvent}^{coulomb} == -60854.3, E_{total} == -17466.7 }

The energy difference comes from protein:solvent, solute:protein & 
solute:solvent. 
The other terms are correct, so the bug is likely in the class InterGroupCLJFF. 

This also explains why the solute is able to drift too close to the protein 
atoms, this is most evident if visualizing snapshots of longer MC simulations.
Original issue reported on code.google.com by julienm...@googlemail.com on 12 Nov 2010 at 4:42
PiRSquared17 / sire

incorrect evaluation of system energy #3
