Correct me if I'm wrong, but your main concern seems to be that the
internal energy contributions for the pairwise decomposition are all 0 and
so do not add up. The rest of the terms *do* add up the way you expect.
Is that correct?
If so, the explanation is pretty straightforward: Amber doesn't compute the
pairwise decomposition for internal terms. You can check inside the
$AMBERHOME/AmberTools/src/sander/ene.F90 source code file and you will see
that `decpair()` (the subroutine used to trigger the aggregation of
decomposed energy contributions) is only called when idecomp == 1 or 2. I
think this choice makes sense.
Most of the internal potential functions are not pairwise decomposable.
Only bonds are pairwise potentials -- angles, torsions, and coupled-torsion
correction map potentials are 3-, 4-, and 5-body potentials, respectively.
For per-atom (or per-residue) decomposition, you can simply divide the
contribution of each term by the number of atoms in the term. For example,
each atom in a valence angle is assigned 1/3 of the total potential energy
in that angle, each atom in a torsion term is assigned 1/4 of the total
potential energy in that term, and so on. There are obvious shortcomings:
imagine distorting just the middle atom in a valence angle from its
equilibrium position -- a strong argument could be made that that middle
atom deserves all of the "credit" for the strain energy. But simply
dividing the contributions equally is conceptually simple, well-defined,
and trivial to implement.
The situation is much more complex in a pairwise decomposition scheme. In
a valence angle, there are 3 distinct pairs of interactions (atoms 1&2,
atoms 1&3, atoms 2&3). In a torsion potential, there are 6, and in coupled
torsion correction maps there are 10. To make things worse, the most
common application of this energy decomposition is for MM/PBSA analyses
where internal potential energy terms have even less meaning (they often
cancel completely in these applications, and even for 3-trajectory
protocols where bound/unbound configurations are not identical, it's
usually the torsions that show the most distortion). It's certainly
possible to craft a scheme in which the internal potential energy
contributions are divided up between pairwise interactions and *do*
properly add up to per-residue and total contributions, but I don't think
there's an obvious "best" scheme, nor do I think any such scheme would even
be useful for anything.
The nonbonded potentials, on the other hand, *are* more pairwise
decomposable, since the interactions are computed over pairs of atoms.
I'll issue my usual caveat here that the solvation free energy contribution
are *not* strictly pairwise decomposable, because they rely on quantities
that are computed over all atoms (effective GB radius for GB and the
dielectric boundary for PB). Nevertheless, there are qualitatively useful
applications for the resulting pairwise decomposition of the nonbonded
potentials.
-----
tl;dr - this isn't a bug, it's pairwise-decomposed bond/angle/torsion/cmap
terms are not computed and so are always exactly 0.
HTH,
Jason
On Thu, Feb 17, 2022 at 12:24 PM Tue Boesen <alyflex.gmail.com> wrote:
> Hi Jason,
>
> Thank you for the reply!
>
> I have attached example.zip to this folder, this contains a pdb and prmtop
> file, and the command that I have run these files with in order to get the
> output. Furthermore there are two folders in the zip file, pairwise and
> per_residue, which contains the input configuration file for both runs
> (they are identical except for idecomp = 2 or 4), and the output files.
>
> The FINAL_RESULTS_MMPBSA.dat for both pairwise and per_residue contains
> identical output results:
>
> Energy Component Average Std. Dev. Std. Err. of
> Mean
>
> -------------------------------------------------------------------------------
> BOND 4.9069 0.0000
> 0.0000
> ANGLE 39.7646 0.0000
> 0.0000
> DIHED 363.0338 0.0000
> 0.0000
> VDWAALS 638.2793 0.0000
> 0.0000
> EEL -2049.8394 0.0000
> 0.0000
> 1-4 VDW 149.2980 0.0000
> 0.0000
> 1-4 EEL 1873.3509 0.0000
> 0.0000
> EGB -426.9272 0.0000
> 0.0000
> ESURF 15.0242 0.0000
> 0.0000
>
> G gas 1018.7941 0.0000
> 0.0000
> G solv -411.9030 0.0000
> 0.0000
>
> TOTAL 606.8911 0.0000
> 0.0000
>
> *For the per_residue case:*
> energy.csv contains the energy decomposition per residue, with lines
> looking like this:
> _______________________________________
> Generalized Born Decomposition Energies
>
> Total Energy Decomposition:
> Frame #,Residue,Internal,van der Waals,Electrostatic,Polar
> Solvation,Non-Polar Solv.,TOTAL
> 1,1,10.936,1.564,-25.44,-84.262,1.2995063999999998,-95.9024936
> 1,2,13.261,-0.122,7.857,-5.835,0.46726559999999995,15.628265599999997
> 1,3,12.866,2.025,13.57,-3.245,0.114948,25.330947999999996
> ________________________________________
> In this case we have various energy components: Internal, van der Waals,
> Electrostatic, Polar solvation, Non-Polar Solv, and finally TOTAL, which is
> just the sum of all the others.
> If I sum up TOTAL for all residues in this case, I get 606.8911 which
> matches the total in FINAL_RESULTS_MMPBSA.dat
> Similarly I can do this with all the other components:
>
> - internal energy = BOND+ANGLE+DIHED
> - van der Waals = VDWAALS+1-4 VDW
> - Electrostatic = EEL + 1-4 EEL
> - Polar solvation = EGB
> - Non-Polar Solv = ESURF
>
> *For the pairwise case:*
> energy.csv contains the energy decomposition per pairwise residues, with
> lines looking like this:
> ___________________________________________
> Generalized Born Decomposition Energies
>
> Total Energy Decomposition:
> Frame #,Resid 1,Resid 2,Internal,van der Waals,Electrostatic,Polar
> Solvation,Non-Polar Solv.,TOTAL
> 1,1,1,0.0,1.938336,-29.877415,-72.295107,2.7776784095999996,-97.4565075904
> 1,1,2,0.0,0.320052,-6.12456,-1.184676,-1.2055162104,-8.1947002104
> 1,1,3,0.0,-0.251325,-1.732596,1.432641,-0.119385396,-0.670665396
> 1,1,4,0.0,-0.066491,0.3531,-0.595142,0.0,-0.30853299999999995
> ___________________________________________
> In this case we have the same energy components, but the internal energy is
> now 0 at all points.
> If I take the van der Waals energy component in this file. I can sum up all
> the contributions where Resid 1 = 1. If I do this, I get exactly the van
> der Waals energy for residue 1 in the per_residue case, and similarly this
> means that if I sum up all van der Waals components for all residue pairs,
> I get the VDWAALS energy in the FINAL_RESULTS_MMPBSA.dat
> Similarly for electrostatic, Polar solvation, and non-polar solvation.
> However for the internal energy I of course do not get this, since all
> components are zero in the pairwise case.
> And similar problems arise for the total energy, since this is just a sum
> of all the previous components.
>
> Hence when I sum the total energy in the energy.csv file for pairwise
> residues I get a different energy than the total energy in
> FINAL_RESULTS_MMPBSA.dat, since I am missing the internal energy.
>
> I hope that explained the problem better for completeness sake I include
> the following python script that can read in the energy components from the
> csv files to make comparisons easier:
>
> import numpy as np
> import matplotlib.pyplot as plt
> def read_mmpbsa_results(file_energy):
> """
> This script reads the enregydeo.csv files produced by mmpbsa in
> ambertools.
> """
> with open(file_energy, 'r') as f:
> lines = f.readlines()[4:]
> internal = []
> vdw = []
> eel = []
> egb = []
> esurf = []
> total = []
> for line in lines:
> dat = line.strip("\n").split(',')
> if len(dat) < 8:
> break
> internal.append(float(dat[2]))
> vdw.append(float(dat[3]))
> eel.append(float(dat[4]))
> egb.append(float(dat[5]))
> esurf.append(float(dat[6]))
> total.append(float(dat[7]))
> internal = np.asarray(internal)
> vdw = np.asarray(vdw)
> eel = np.asarray(eel)
> egb = np.asarray(egb)
> esurf = np.asarray(esurf)
> total = np.asarray(total)
> return internal, vdw, eel, egb, esurf, total
>
> def read_mmpbsa_results_pair(file_energy):
> """
> This script reads the enregydeo.csv files produced by mmpbsa in
> ambertools.
> """
> with open(file_energy, 'r') as f:
> lines = f.readlines()[4:]
> internal = []
> vdw = []
> eel = []
> egb = []
> esurf = []
> total = []
> nres = 0
> for line in lines:
> dat = line.strip("\n").split(',')
> if len(dat) < 8:
> break
> nres = max(nres, int(dat[1]))
> internal.append(float(dat[3]))
> vdw.append(float(dat[4]))
> eel.append(float(dat[5]))
> egb.append(float(dat[6]))
> esurf.append(float(dat[7]))
> total.append(float(dat[8]))
> internal = np.asarray(internal).reshape(nres,nres)
> vdw = np.asarray(vdw).reshape(nres,nres)
> eel = np.asarray(eel).reshape(nres,nres)
> egb = np.asarray(egb).reshape(nres,nres)
> esurf = np.asarray(esurf).reshape(nres,nres)
> total = np.asarray(total).reshape(nres,nres)
> return internal, vdw, eel, egb, esurf, total
>
>
>
>
>
>
>
>
> file_per_residue = '/home/tue/test/energy.csv'
> file_pairwise = '/home/tue/test/1/energy.csv'
> # Eall_pair,Eside_pair,Eback_pair = read_mmpbsa_results_pair(file_pairwise)
> internal1, vdw1, eel1, egb1, esurf1, total1 =
> read_mmpbsa_results_pair(file_pairwise)
> internal, vdw, eel, egb, esurf, total =
> read_mmpbsa_results(file_per_residue)
>
> # plt.imshow(np.log(np.abs(total1)+1e-9))
> # plt.colorbar()
> # plt.show()
>
>
> plt.imshow(vdw1)
> plt.colorbar()
> plt.show()
>
> internal1_sum = np.sum(internal1,axis=1)
> vdw1_sum = np.sum(vdw1,axis=1)
> eel1_sum = np.sum(eel1,axis=1)
> egb1_sum = np.sum(egb1,axis=1)
> esurf1_sum = np.sum(esurf1,axis=1)
> total1_sum = np.sum(total1,axis=1)
>
> np.abs(internal - internal1_sum) < 0.001
> np.abs(vdw - vdw1_sum) < 0.001
> np.abs(eel - eel1_sum) < 0.001
> np.abs(egb - egb1_sum) < 0.001
> np.abs(esurf - esurf1_sum) < 0.001
> np.abs(total - total1_sum) < 0.001
>
> Kind regards
> Tue
>
> On Thu, Feb 17, 2022 at 7:19 AM Jason Swails <jason.swails.gmail.com>
> wrote:
>
> > Hi Tue,
> >
> > The decomposition code has not changed in years.
> >
> > It's unclear to me exactly what you're observing (which is why sending
> > actual output/files is so helpful). Examples of the questions I have are
> > posted in-line
> >
> > On Sun, Feb 13, 2022 at 9:47 PM Tue Boesen <alyflex.gmail.com> wrote:
> >
> > > I believe there is a bug when running mmpbsa with idecomp setup for
> > > pairwise energies (idecomp=3 or idecomp=4).
> > >
> > > In this mode the internal energies are reported at zero, while all
> other
> > > energies are given with reasonable values that add up to their total
> > value
> > > and matches the results you would get if you ran it with idecomp=1-2
> (per
> > > residue).
> > >
> >
> > What values are you looking at? The complex/receptor/ligand, or the
> > binding totals (that is, the difference between the bound and free
> > states)? All of them?
> >
> > For instance, the internal terms *should* always be zero for total
> binding
> > terms in a single-trajectory protocol, since the bonded geometries are
> > identical in the bound and unbound states. They are also likely to be
> very
> > close to zero in the residue-pairwise decomposition, since only those
> terms
> > that bridge two residues would actually contribute to this term, and if
> > they are close to equilibrium values the total energies will not be very
> > large.
> >
> > I also don't understand what you mean by "matches the results you would
> get
> > if you ran it with idecomp=1,2". These compute fundamentally different
> > quantities, so the numbers should be quite different. The Amber force
> > field potential is not pairwise decomposable, so idecomp=3,4 are only
> > approximations, anyway.
> >
> > The error appears in both the csv, and final decomposition output and
> means
> > > that the energies do not add up as they are supposed to. The energies
> in
> > > FINAL_RESULTS_MMPBSA, still include the summed total of the internal
> > > energies, and give identical results as when ran with idecomp=1-2.
> > >
> >
> > How do you expect the terms to add up? This is where actual output files
> > and precise descriptions help.
> >
> > Good luck,
> > Jason
> >
> > --
> > Jason M. Swails
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber
> >
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
--
Jason M. Swails
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Feb 18 2022 - 07:00:02 PST
