That is very true Thomas.
I wasn't aware that MMPBSA was open source, such that I even had the option
of looking under the hood. If it has been written in Fortran 90 I might
take a look, since I did know fortran very well half a decade ago at
least...
I assume that the same is true for the ambertools manual then?
Because if nothing else, I would add something like Jasons explanation to
the manual under idecomp 3-4, since there really is nothing in there saying
that it won't calculate internal energy or why, nor even a warning that the
total energy is actually not the full total energy.
In any case, thank you both for your feedback
Tue
On Fri, Feb 18, 2022 at 8:48 PM Thomas Cheatham <tec3.utah.edu> wrote:
>
> Poor designs or bugs can be fixed ... we always welcome volunteers and/or
> you can write your own analysis tools. My first pass with MM-PBSA was Perl
> scripts. Later scripts to validate other's scripts. Research requires
> constant validation and assessment - trust, but verify. If you've got
> better ideas, please implement them and convince the community. --tec3
>
> ________________________________________
> From: Tue Boesen <alyflex.gmail.com>
> Sent: Friday, February 18, 2022 3:38 PM
> To: AMBER Mailing List
> Subject: Re: [AMBER] Bug in mmpbsa.py internal energy when idecomp is 3 or
> 4
>
> 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
>
> _______________________________________________
> 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 Sat Feb 19 2022 - 10:00:02 PST
