Re: [AMBER] Bug in mmpbsa.py internal energy when idecomp is 3 or 4

From: Jason Swails <jason.swails.gmail.com>
Date: Tue, 22 Feb 2022 12:23:36 -0500

On Fri, Feb 18, 2022 at 5:38 PM Tue Boesen <alyflex.gmail.com> wrote:

> 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.
>

There are many things you would find surprising about the pairwise
numbers. For instance, if you deleted all residues except residues 2 and
3, you would naively expect the pairwise interaction between the two of
them to remain unchanged. While this will be true for some of the terms
(like lennard-jones), it will not be true for the solvation terms which are
themselves many-body terms.

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.
>

But that *is* an option., meaning there's a viable workaround for those
instances where the decomposed internal energy contributions are either
qualitatively or quantitatively useful (I can't think of any situation
where those values would be useful outside of satisfying this sanity check).

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.


Consistently only in a narrow sense (per-res internal with pairwise
nonbonded does not sound consistent to me). Correctly, arguably not.
Residues are often bonded to each other, and residues connected to each
other shouldn't have non-zero pairwise internal energy contributions if you
claim to be computing such contributions. I don't actually this proposal
as an improvement over the current behavior (outside of the internal energy
omissions being undocumented right now).

Whether

(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.)


But angles and torsions are not, and those are more numerous (and
important) than bonds. If you only did bonds, the totals would still not
add up because angles and torsions would be omitted. If you did pairwise
bonds, but per-residue angles and torsions, what would those values even
mean?

I quibble here to point out that any proposal on how to perform pairwise
decomposition on a potential energy function that is not pairwise
decomposable is impossible to define in a way that most people would agree
is correct.

I definitely agree that the implemented behavior should be documented (and
once documented, it becomes a feature :-D). I will add these notes.

Thanks,
Jason

-- 
Jason M. Swails
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Feb 22 2022 - 09:30:03 PST
Custom Search