Ligand only pairwise per-residue energy decomposition syntax

[Harold Grosjean]

Hello,

I am trying to calculate the interaction energy of a ligated (only) against all protein residues to identify favorable and/ or unfavorable binding regions.  To do so, I am trying to modify the provide "print_res" section so that it performs the aforementioned calculation. The calculation works well when I use print_res="within 4". I have tried multiple combination (see infile for more details) of print_res= using this command:

mpirun -np 5 gmx_MMPBSA MPI -O -i protein_ligand_ST.in -cs F11-PHIPA-x20023.tpr -ci protein_ligand_sep.ndx -cg 24 25 -ct md_fit.xtc -cp ../complex/setup/complex.top -lm ../model_building/LIG_antechamber.mol2

but it always returns the following error:

[INFO   ] Starting
[INFO   ] Command-line
  gmx_MMPBSA -O -i protein_ligand_ST.in -cs F11-PHIPA-x20023.tpr -ci protein_ligand_sep.ndx -cg 24 25 -ct md_fit.xtc -cp ../complex/setup/complex.top -lm ../model_building/LIG_antechamber.mol2

[WARNING] protein_forcefield and ligand_forcefield variables are deprecate since version 1.4.1 and will be remove in the next version. Please, use forcefield instead.
[WARNING] entropy_seg variable is deprecate since version 1.4.2 and will be remove in v1.5.0. Please, use ie_segment instead.
[INFO   ] Checking external programs...
[INFO   ] cpptraj found! Using /biggin/b193/lina3397/anaconda3/envs/followups_docking/bin/cpptraj
[INFO   ] tleap found! Using /biggin/b193/lina3397/anaconda3/envs/followups_docking/bin/tleap
[INFO   ] parmchk2 found! Using /biggin/b193/lina3397/anaconda3/envs/followups_docking/bin/parmchk2
[INFO   ] sander found! Using /biggin/b193/lina3397/anaconda3/envs/followups_docking/bin/sander
[INFO   ] Using GROMACS version > 5.x.x!
[INFO   ] gmx found! Using /sbcb/packages/opt/Linux_x86_64/gromacs/2020.3_GCC6.2_CUDA10.1.AVX_512//bin/gmx
[INFO   ] Checking external programs...Done.

[INFO   ] Building AMBER Topologies from GROMACS files...
[INFO   ] Checking gmxMMPBSA data folder exists in Amber data...
[INFO   ] Get PDB files from GROMACS structures files...
[INFO   ] Making gmx_MMPBSA index for complex...
[INFO   ] Normal Complex: Saving group 24_25 in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_COM.pdb
[INFO   ] Generating ligand parameters from ../model_building/LIG_antechamber.mol2 file...
[INFO   ] No receptor structure file was defined. Using ST approach...
[INFO   ] Using receptor structure from complex to generate AMBER topology
[INFO   ] Normal Complex: Saving group 24 in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_REC.pdb
[INFO   ] No ligand structure file was defined. Using ST approach...
[INFO   ] Using ligand structure from complex to generate AMBER topology
[INFO   ] Normal ligand: Saving group 25 in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_LIG.pdb
[INFO   ] Building Normal Complex Amber Topology...
[INFO   ] Writing Normal Complex Amber Topology...
[INFO   ] No Receptor topology files was defined. Using ST approach...
[INFO   ] Building AMBER Receptor Topology from Complex...
[INFO   ] No Ligand Topology files was defined. Using ST approach...
[INFO   ] Building AMBER Ligand Topology from Complex...
[ERROR  ] MMPBSA_Error We expected something like this: A/2-10,35,41 but we get B/1435;A/1316-1434 instead.

I would be grateful if you could indicate how to get the correct synthax so that I can run this calculation. Please find attached relevant files files.

ideally, I would like to reproduce figure 4 of this webpage:

https://valdes-tresanco-ms.github.io/gmx_MMPBSA/analyzer/#per-residue-decomposition

but with only 1 line which should speed up the calculation too.

Additional question:

I see that two options concerning the idecomp for pairwise analysis: idecomp=3 or idecomp=4. Which of the 2 would you pick in order to quantify residue energetic contributions to ligand binding.

Please, let me know if you need clarification.

I thank you very much in advance for your helps,

Regards,

Harold Grosjean

[Harold Grosjean]

Update:

Sometimes the error changes to:

  File "/biggin/b193/lina3397/anaconda3/envs/followups_docking/bin/gmx_MMPBSA", line 8, in <module>
    sys.exit(gmxmmpbsa())
  File "/biggin/b193/lina3397/anaconda3/envs/followups_docking/lib/python3.7/site-packages/GMXMMPBSA/app.py", line 95, in gmxmmpbsa
    app.make_prmtops()
  File "/biggin/b193/lina3397/anaconda3/envs/followups_docking/lib/python3.7/site-packages/GMXMMPBSA/main.py", line 576, in make_prmtops
    self.FILES.mutant_receptor_prmtop, self.FILES.mutant_ligand_prmtop) = maketop.buildTopology()
  File "/biggin/b193/lina3397/anaconda3/envs/followups_docking/lib/python3.7/site-packages/GMXMMPBSA/make_top.py", line 122, in buildTopology
    self.reswithin()
  File "/biggin/b193/lina3397/anaconda3/envs/followups_docking/lib/python3.7/site-packages/GMXMMPBSA/make_top.py", line 528, in reswithin
    dist, exclude, res_selection = selector(self.INPUT['print_res'])
  File "/biggin/b193/lina3397/anaconda3/envs/followups_docking/lib/python3.7/site-packages/GMXMMPBSA/utils.py", line 147, in selector
    ri = [chain, int(ci[0]), ''] if len(ci) == 1 else [chain, int(ci[0]), ci[1]]
ValueError: invalid literal for int() with base 10: '1435;'
Error occured on rank 0.
Exiting. All files have been retained.

When I use: print_res="B/1435; A/1316-1434 in the infile

Thanks for the help,

Regards,

Harold

[marioe911116]

Hello Harold! try with this:

print_res="A/1316-1434 B/1435"

The variable idecomp determines what type of decomposition is going to be performed... if you want per-residue, then idecomp=2, if you want per-wise (computationally more expensive, and more specific), then idecomp=4. Take into account that the residues beyond 6~7 angstroms from the ligand generally don't have a significant contribution and then the per-residue energy will be close to zero. That's why we implement (and prefer) going with the "within 6" option. Last but not least, when performing per-wise decomposition, the number of calculations scales really quick and can take quite some time and generate heavy files.

I strongly recommend to check the decomposition session in the documentation (https://valdes-tresanco-ms.github.io/gmx_MMPBSA/input_file/#decomp-namelist-variables) to get better clarity regarding notations, variables, etc...

cheers!

Mario E.