Issue with Non-Sequential Residue Selection in gmx_MMPBSA (IndexError: List Index Out of Range)
James Davies
I hope you're doing well. I am encountering an issue while running gmx_MMPBSA on my system, and I would appreciate your assistance in resolving it.
The gmx_MMPBSA command works fine when I use a custom index file with sequential residues (e.g., residues 33 and 34), but I encounter an error when I modify the index file to include non-sequential residues that make up my binding site. The error I get is:
gmx_MMPBSA -O -i mmpbsa.in -cs new.tpr -ct 10.xtc -ci binding.ndx -cg 33 34 -cp CuO_solv.top -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv -do FINAL_DECOMP_MMPBSA.dat -deo FINAL_DECOMP_MMPBSA.csv
tets
[INFO ] Checking mmpbsa.in input file...
[INFO ] Checking mmpbsa.in input file...Done.
[INFO ] Checking external programs...
[INFO ] cpptraj found! Using /home/c.c1842421/.conda/envs/gmxMMPBSA/bin/cpptraj
[INFO ] tleap found! Using /home/c.c1842421/.conda/envs/gmxMMPBSA/bin/tleap
[INFO ] parmchk2 found! Using /home/c.c1842421/.conda/envs/gmxMMPBSA/bin/parmchk2
[INFO ] sander found! Using /home/c.c1842421/.conda/envs/gmxMMPBSA/bin/sander
[INFO ] Using GROMACS version > 5.x.x!
[INFO ] gmx found! Using /home/c.c1842421/.conda/envs/gmxMMPBSA/bin.AVX2_256/gmx
[INFO ] Checking external programs...Done.
[INFO ] Building AMBER topologies from GROMACS files...
[INFO ] Get PDB files from GROMACS structures files...
[INFO ] Making gmx_MMPBSA index for complex...
[INFO ] Normal Complex: Saving group coord_CU1_CU2 (33_34) 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 Receptor: Saving group coord (33) 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 CU1_CU2 (34) in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_LIG.pdb
[INFO ] Checking the structures consistency...
[WARNING] No reference structure was found and the complex structure not contain any chain ID. Assigning chains ID automatically...
File "/home/c.c1842421/.conda/envs/gmxMMPBSA/bin/gmx_MMPBSA", line 8, in <module>
sys.exit(gmxmmpbsa())
File "/home/c.c1842421/.conda/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/app.py", line 98, in gmxmmpbsa
app.make_prmtops()
File "/home/c.c1842421/.conda/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/main.py", line 682, in make_prmtops
self.FILES.mutant_receptor_prmtop, self.FILES.mutant_ligand_prmtop) = maketop.buildTopology()
File "/home/c.c1842421/.conda/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 123, in buildTopology
self.gmx2pdb()
File "/home/c.c1842421/.conda/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 514, in gmx2pdb
self.check_structures(self.complex_str, self.receptor_str, self.ligand_str)
File "/home/c.c1842421/.conda/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 1414, in check_structures
self._assign_chains_IDs(com_str, rec_str, lig_str)
File "/home/c.c1842421/.conda/envs/gmxMMPBSA/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 1473, in _assign_chains_IDs
lig_str.residues[i].chain = res.chain
IndexError: list index out of range
Exiting. All files have been retained.
This isn’t an issue with chain IDs, as my previous calculations, which include the entire protein and my ligand as separate index groups, worked fine. When I use the residues from my binding.ndx index file (which contains non-sequential residues 49, 51, 71....), I encounter this error. However, the command works perfectly fine if I use residues 33 and 34 in sequential order.
I’ve tested this thoroughly, and the issue only arises when the residues are non-sequential. Could you kindly advise on how to overcome this error when working with non-sequential residues in the .ndx file?
If further information is required, please let me know, and I will be happy to provide it.
Thank you for your time and assistance.
Best regards,
James
mariosergi...@gmail.com
This error is related to automatic chain assignment. If you look at the log, there is a warning that your structure has no chains so it will be assigned automatically. Since these are discontinuous residues, the mapping (annotation of the index of each atom) and the ligand index do not match. When the user does not provide a reference structure, we use a continuity method for chain assignment. If the ligand contains continuous residues, then the same chain is assigned to it, if they are discontinuous, then two different chains are attempted to be assigned. This should be fixed by adding the -cr reference_structure.pdb option. This file should have all the residues mapped to their respective chains.
James Davies
Thank you for your advice! It worked – gmx_mmpbsa is now correctly recognizing the residues I selected for the calculation.
To give you more context, I’m investigating the energetics of an active site in a copper-containing metalloprotein, specifically comparing the oxygenated and deoxygenated states. My focus is on the copper coordination site and the surrounding residues to explore how the energetics differ between these two states.
For the calculations, I’m aiming to determine the optimal dielectric constant for the system. I plan to run calculations with different dielectric constants (starting with values like 1, 3, and 5) and analyze the results to identify the value that offers the best balance of accuracy and consistency. This will be done both for the copper site and the coordinating residues. So far, I’ve observed significant variation in total energy, with higher dielectric constants yielding lower complex total energy (positive) and larger delta total energy (negative). I’d love to hear your interpretation of these results, and whether you think I should focus purely on the binding site or consider the whole protein, with decomposition performed on the binding site.
Do you think this approach is suitable for this kind of system, or would you suggest any modifications? Specifically, I’d appreciate your thoughts on the increments I should use between dielectric constants (e.g., 1 or 0.5) and what the upper boundary should be, especially given the influence of the copper coordination site.
Looking forward to your insights!
Best regards,
James