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
Received on Thu Feb 17 2022 - 09:30:02 PST
