Hi Jason,
I appreciate the insight!
However, I would still classify this as a bug or at least as a somewhat
poor design choice for several reasons:
1)
The fact that the total energies no longer add up to the actual total
energy seems wrong.
2)
When you run the code with idecomp 3,4 you get no internal energy at all,
so if you still want internal energy (even if you can only get it per
residue, the only option you have for this is to run mmpbsa with both
idecomp 1-2 and 3-4, and then read in results from several different output
files and merge these.
If the internal energy does not make much sense decomposing for pairwise
residues, then I would be fine with it just being displayed as a single per
residue term in the pairwise decomposition (so the internal energy would be
zero for all residue pairs (i,j), except when i=j in which case it would
show the per residue internal energy for residue i. In that case the
internal energy would still be available, and the total energy would be
computed consistently. (Though as you also mentioned the bonded term should
be a pairwise term, so I would expect that one to be somewhat easily and
sensible to decompose for pairwise residues.)
Kind regards
Tue
On Fri, Feb 18, 2022 at 6:46 AM Jason Swails <jason.swails.gmail.com> wrote:
> 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
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Feb 18 2022 - 15:00:02 PST
