Failure in MT protein-protein interaction

[wrmartin]

I am trying to do PBSA on a protein-protein complex using the CHARMM forcefield, but I'm getting numerous errors. I've been able to "fix" some of them, but I'm stuck here. The receptor for the system is also a multimer, and using a single-trajectory approach on the receptor only (and also the complex) work fine. A multi-trajectory approach on the receptor also fails with a similar error. The logfile is as follows:

[INFO   ] Starting
[INFO   ] Command-line
  gmx_MMPBSA -O -i mmpbsa.in -cs ../complex.pdb -ci ../complex.ndx -cg 32 34 -ct ../complex.xtc -cp complex.top -rs ../receptor.pdb -ri ../receptor.ndx -rg 32 -rt ../receptor.xtc -rp receptor.top -ls ../ligand.pdb -li ../ligand.ndx -lg 26 -lt ../ligand.xtc -rp ligand.top -nogui

[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 /cm/shared/apps/amber/21/bin/cpptraj
[INFO   ] tleap found! Using /cm/shared/apps/amber/21/bin/tleap
[INFO   ] parmchk2 found! Using /cm/shared/apps/amber/21/bin/parmchk2
[INFO   ] sander found! Using /cm/shared/apps/amber/21/bin/sander
[INFO   ] Using GROMACS version > 5.x.x!
[INFO   ] gmx found! Using /cm/shared/apps/gromacs/2021.1/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 32_34 in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_COM.pdb
[INFO   ] A receptor structure file was defined. Using MT approach...
[INFO   ] Normal receptor: Saving group 32 in ../receptor.ndx file as _GMXMMPBSA_REC.pdb
[INFO   ] A ligand structure file was defined. Using MT approach...
[INFO   ] Normal Ligand: Saving group 26 in ../ligand.ndx file as _GMXMMPBSA_LIG.pdb
[INFO   ] Building Normal Complex Amber Topology...
[INFO   ] Writing Normal Complex Amber Topology...
[INFO   ] A Receptor topology file was defined. Using MT approach...
[INFO   ] Building AMBER Receptor Topology from GROMACS Receptor Topology...
  File "/cm/shared/apps/amber/21/miniconda/bin/gmx_MMPBSA", line 8, in <module>
    sys.exit(gmxmmpbsa())
  File "/cm/shared/apps/amber/21/miniconda/lib/python3.9/site-packages/GMXMMPBSA/app.py", line 95, in gmxmmpbsa
    app.make_prmtops()
  File "/cm/shared/apps/amber/21/miniconda/lib/python3.9/site-packages/GMXMMPBSA/main.py", line 576, in make_prmtops
    self.FILES.mutant_receptor_prmtop, self.FILES.mutant_ligand_prmtop) = maketop.buildTopology()
  File "/cm/shared/apps/amber/21/miniconda/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 117, in buildTopology
    tops = self.gmxtop2prmtop()
  File "/cm/shared/apps/amber/21/miniconda/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 349, in gmxtop2prmtop
    rec_top.coordinates = self.receptor_str.coordinates
  File "/cm/shared/apps/amber/21/lib/python3.9/site-packages/parmed/structure.py", line 1699, in coordinates
    coords = coords.reshape((-1, len(self.atoms), 3))
ValueError: cannot reshape array of size 83385 into shape (27899,3)
Error occured on rank 0.
Exiting. All files have been retained.
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
[warn] Epoll MOD(1) on fd 48 failed.  Old events were 6; read change was 0 (none); write change was 2 (del): Bad file descriptor
[warn] Epoll MOD(4) on fd 48 failed.  Old events were 6; read change was 2 (del); write change was 0 (none): Bad file descriptor

For what it's worth, the receptor is 27795 atoms; I'm not sure where the 27899 is coming from.

[marioernestovaldes]

Hello there! Usually the fastest way for us to deal with these errors is to send around the files you are using, so we can check closely... Could you pleae send a compressed file with all the files you are using? Make sure to include all the itp files associated with the topologies... Send only few frames in the trajectories...

Cheers!

[wrmartin]

Best I can do is a link to some files; it's too big to attach here (the trajectories aren't the problem).

https://1drv.ms/u/s!AhqzCZUniK1ShYY7h7gSP-V_la8LOw?e=c7TN7G

[Valdes-Tresanco-MS]

Hi Martin The error is related to the number of atoms in the topology and structure file. The structure file only contains the atoms of the selected group (27795) while the topology also contains the atoms that are not removed, for example, CAL, ZN2, SOD, POPS Currently, only the following groups of atoms are automatically removed:

sol_ion = [
     # standard gmx form
     'NA', 'CL', 'SOL',
     # charmm-GUI form ??
     'SOD', 'CLA', 'TIP3P', 'TIP4P', 'TIPS3P', 'TIP5P', 'SPC', 'SPC/E', 'SPCE', 'TIP3o', 'WAT']

The way to fix this is by eliminating those molecules and leaving only the ones of interest, for example:

#include "topparr/forcefield.itp"
#include "topparr/PROA.itp"
#include "topparr/PROB.itp"
#include "topparr/PROC.itp"

[ system ]
; Name
Title

[ molecules ]
; Compound  #mols
PROA               1
PROB               1
PROC               1

We are working to improve the handling of topologies, given that in some cases eliminating the molecules is not a viable option and focused on the implementation of the consideration of explicit water molecules. However, there is bad news and that is that the current stable version (v1.4.3) has a bug for the MT approach. Currently, we fixed it in part but the development version (v1.5.0) is not ready for production yet. With your system, we managed to identify another problem related to this, so we hope that the new version will be able to do the calculations without difficulty. We are making a huge effort to release the v1.5.0 as soon as possible. We regret not being able to give you the solution immediately, but we assure you that as soon as the new version is released, you will be able to do your calculations without problems. Sincerely.

[wrmartin]

I don't believe there is an issue regarding the first part. As I noted, the complex works just fine in the ST approach with the files as they are, and only fails when using the MT approach. Sounds like I will just have to be patient for the update, thank you!

[Valdes-Tresanco-MS]

Ok, I have identified the error. In version 1.4.3 it is still impossible to carry out the calculation since an error related to the way the topology is generated from the index persists.

You are defining the topology of the ligand with the -rp option. Currently, gmx_MMPBSA does not check if the option is duplicated, which means that ligand.top is taken as an argument of -rp, which will always cause a mismatch between the structure and the topology of the receptor. In version 1.4.3, although the modification of the topology is not necessary, its manipulation is quite limited. In the current version, no modification is necessary and we have much more robust management of the topology, so there should be no future related errors. We have corrected these problems in the development version. I'll let you know when we're in the final stages of testing for you to install it.

[Valdes-Tresanco-MS]

Fixed via f6a3ad0, 6ebdeab

[Valdes-Tresanco-MS]

We are in the final stages of testing the new version. You can test if it runs fine for your system. Our tests did not give errors, however, your verification is important to us. You can test the development version as follows: conda create --name AmberTools21 python=3.9 conda activate AmberTools21 conda install -c conda-forge mpi4py ambertools=21 compilers python -m pip install PyQt5 python -m pip install git+https://github.com/Valdes-Tresanco-MS/gmx_MMPBSA

I will close this issue, but if you have any doubts or errors, feel free to open a new one or reopen this one.