Error running gmx_MMPBSA

[Nathaniel STILLSON]

Im using CHARMM36m and relatively new to this. I am trying to get results for a protein ligand system. The ligand was prepared prior to MD with CGenFF. There was a boundary condition with the results from my MD which I tried to resolve using trjconv -pbc no jump. When I use the .tpr file as the structure input, the receptor still appears broken in the _GMXMMPBSA outputs; however, when I use the .gro file this is no longer an issue. In either case I still get the following error: 

Thanks in advance

gmx_MMPBSA -O -i mmpbsa.in -cs gsk_1ns.gro -ci index.ndx -cg 1 13 -ct gsk_corrected.xtc -cp topol.top
[INFO   ] Starting
[INFO   ] Checking external programs...
[INFO   ] cpptraj found! Using /home/njs00/anaconda3/envs/amber/bin/cpptraj
[INFO   ] tleap found! Using /home/njs00/anaconda3/envs/amber/bin/tleap
[INFO   ] parmchk2 found! Using /home/njs00/anaconda3/envs/amber/bin/parmchk2
[INFO   ] sander found! Using /home/njs00/anaconda3/envs/amber/bin/sander
[INFO   ] Using GROMACS version > 5.x.x!
[INFO   ] gmx found! Using /usr/bin/gmx
[INFO   ] Checking external programs...Done.

[INFO   ] Building AMBER Topologies from GROMACS files...
[INFO   ] Checking if supported force fields exists in Amber data...
[INFO   ] Get PDB files from GROMACS structures files...
[INFO   ] Making gmx_MMPBSA index for complex...
[INFO   ] Normal Complex: Saving group 1_13 in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_COM.pdb
[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 1 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 13 in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_LIG.pdb
[WARNING] No reference structure was found and a gro file was used for the complex structure. Assigning chains ID...
[INFO   ] Cleaning GROMACS topologies...
[INFO   ] Building Normal Complex Amber Topology...
/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/parmed/gromacs/gromacstop.py:1027: GromacsWarning: The [ pairs ] section contains 3 exceptions that aren't 1-4 pairs; make sure you know what you're doing!
  warnings.warn('The [ pairs ] section contains %i exceptions that '
[INFO   ] Writing Normal Complex Amber Topology...
  File "/home/njs00/anaconda3/envs/amber/bin/gmx_MMPBSA", line 8, in <module>
    sys.exit(gmxmmpbsa())
  File "/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/GMXMMPBSA/app.py", line 97, in gmxmmpbsa
    app.make_prmtops()
  File "/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/GMXMMPBSA/main.py", line 577, in make_prmtops
    self.FILES.mutant_receptor_prmtop, self.FILES.mutant_ligand_prmtop) = maketop.buildTopology()
  File "/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 118, in buildTopology
    tops = self.gmxtop2prmtop()
  File "/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 313, in gmxtop2prmtop
    action = ChRad(com_amb_prm, PBRadii[self.INPUT['PBRadii']])
  File "/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/parmed/tools/changeradii.py", line 295, in ChRad
    _call_method[radii_set](parm)
  File "/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/parmed/tools/changeradii.py", line 145, in mbondi2
    if atom.bond_partners[0].atomic_number == 7:
IndexError: list index out of range
Exiting. All files have been retained.

[Mario Sergio Valdes]

Hi and welcome!

Im using CHARMM36m and relatively new to this. I am trying to get results for a protein ligand system. The ligand was prepared prior to MD with CGenFF. There was a boundary condition with the results from my MD which I tried to resolve using trjconv -pbc no jump. When I use the .tpr file as the structure input, the receptor still appears broken in the _GMXMMPBSA outputs; however, when I use the .gro file this is no longer an issue. In either case I still get the following error: 

Let's see, the first problem you mention I opened as an issue some time ago https://github.com/Valdes-Tresanco-MS/gmx_MMPBSA/issues/18. 

If the tpr file has PBC, it is most likely that the solvation box was initially too small, so there will be multiple inconsistencies. In the issue that I mention you can see the reasons why I closed it.

...

 

[INFO   ] Building Normal Complex Amber Topology...
/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/parmed/gromacs/gromacstop.py:1027: GromacsWarning: The [ pairs ] section contains 3 exceptions that aren't 1-4 pairs; make sure you know what you're doing!
  warnings.warn('The [ pairs ] section contains %i exceptions that '

 

This warning is related to a strange element in the topology. I don't know exactly what it is, because there is no way to define it, except you who built the topology.

 

[INFO   ] Writing Normal Complex Amber Topology...
  File "/home/njs00/anaconda3/envs/amber/bin/gmx_MMPBSA", line 8, in <module>
    sys.exit(gmxmmpbsa())
  File "/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/GMXMMPBSA/app.py", line 97, in gmxmmpbsa
    app.make_prmtops()
  File "/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/GMXMMPBSA/main.py", line 577, in make_prmtops
    self.FILES.mutant_receptor_prmtop, self.FILES.mutant_ligand_prmtop) = maketop.buildTopology()
  File "/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 118, in buildTopology
    tops = self.gmxtop2prmtop()
  File "/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 313, in gmxtop2prmtop
    action = ChRad(com_amb_prm, PBRadii[self.INPUT['PBRadii']])
  File "/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/parmed/tools/changeradii.py", line 295, in ChRad
    _call_method[radii_set](parm)
  File "/home/njs00/anaconda3/envs/amber/lib/python3.9/site-packages/parmed/tools/changeradii.py", line 145, in mbondi2
    if atom.bond_partners[0].atomic_number == 7:
IndexError: list index out of range
Exiting. All files have been retained.

Here gmx_MMPBSA cannot assign the PBRadii to a certain atom in the topology. This is probably related to the previous warning.

gmx_MMPBSA may have such limitations for rare atoms, however, this should not be the case. Make sure the files you use are consistent. If your system contains rare atoms, it probably requires other methods that have it parameterized.

If your system does not contain rare atoms, you have no problems with the PBC and it still does not work properly, you can send the files to me for review.

Best!
Mario S.

 

[Nathaniel STILLSON]

Thanks for your help!

(1) From inspection it seems my solvation box is not too small to contain my protein, it just is not centered on it. For future MD's should I center the crystal structure beforehand? The gromacs tutorial made it seem that an off-center box is not a critical error.

(2) This warning only appears when I used the .gro file. From my understanding errors in 1-4 pairs concern atoms separated by (multiple?) bonds that are too far apart to still be realistically bonded. The system I am working with is a heterotetramer and it some files the ends of the chains are connected by a dotted line. When preparing the MD, however, 4 separate .itp files are generated for the four chains of the complex.

(3) The ligand I am working with has a chlorobenzyl group. When put through CGenFF it generated an LP atom bonded to the chlorine (which I am presuming is a lone pair). This LP atom is new to me and I was not clear on its involvement since a chlorine group would be expected to have 3 lone pairs in a tetrahedral configuration. Is this the likely source of error?

I will try to perform a short MD of just the relevant monomer with a different inhibitor to see if that resolves the issue.

Best,

Jake

[Mario Sergio Valdes]

Since we get the PDB files with editconf from the tpr, that the system is centered is almost a requirement

This error is not new to us, however, we have never seen it in this context. Please can you send me a copy of the files you use? (10 frames xtc, tpr, gro, index, and groups). So we can review them

[Nathaniel STILLSON]

I ran a different small molecule on the monomer and everything seems to have worked correctly (although I am still having some issues centering the monomer in the box).The files all together were too large to attach so they are shared here. The protein and ligand groups are 1, 13 respectively.

Thanks

[Mario Sergio Valdes]

Could you please, send me the topology file (and itps)?

[Nathaniel STILLSON]

Ive added those files to the folder.

Thanks

[Mario Sergio Valdes]

We are reviewing the files. We will try to get you an answer as soon as possible.

Best!

Mario S.

[Mario Sergio Valdes]

Hi again.

Please, Include the charmm36-feb2021.ff folder and the gsk.prm file

Best!

[Nathaniel STILLSON]

done!

[Mario Sergio Valdes]

Hello, unfortunately, I have bad news.

Since LP is a pseudo atom and parmed does not support lone pairs for Amber, the topology cannot be built, so gmx_MMPBSA is unable to process your system.

The error that you reported at the beginning, is due to the fact that the topology does not register any bond for that atom so that when gmx_MMPBSA tries to assign the PBRadii for this atom fail since it is non-bonded

Regarding your system and the PBC.

I really don't know if it's related, but in my experience and the tests I did to close the issue that I mentioned earlier, for long systems in one direction the box should be larger. That is, your system has 4 contiguous protein units which make it very long in one direction than the other two. In this same way, the box (if it is cubic) defined will have these characteristics. When pre-production procedures are performed, the molecule can rotate and translate. In case the longest side of the molecule comes out through the narrower side of the box, then it is very difficult to remove the PBC, since Gromacs will try to do it on the structure, which is longer than that side of the box, obtaining a structure still with PBC. The only way that worked to remove the PBC was to center on a residue.

Ways to fix the problems:

- The only possible way to use gmx_MMPBSA is by removing the LP atom. You must run the dynamic without including the LP atom and do different analyzes if necessary.

- Unless you are studying some effect of the ligand on the binding of protein units, I recommend simulating the system as simple as possible. The calculation time increases considerably, both for the simulation of molecular dynamics and for the calculation with the PB model, as well as the possibility of inconsistencies.

- Remove the PBC through centering in res. In this case, we select the residue A:249 (ri 249 in the make_ndx) and it worked perfectly

Best wishes!

Mario S.

[marioe911116]

Hi there! We came up with a solution for those systems that include LPH atoms to get a better representation of the halogen bond... check our new tutorial at:

https://valdes-tresanco-ms.github.io/gmx_MMPBSA/examples/Protein_ligand_LPH_atoms_CHARMMff/

let us know if you have any questions...

best!

Mario E.