Issue with CHARMM-GUI generated system

[Christian Meyer]

Hello,

Thank you for your wonderful software!  I'm getting an error for a membrane-protein-ligand system I generated using CHARMM-GUI software. Attached is the log file.  

The command I am using is:

gmx_MMPBSA -O -i mmpbsa.in -cs step7_1.tpr -ci index.ndx -cg 1 22 -ct step7_1_pbc.xtc -cp topol.top -o FINAL_RESULTS_MMPBSA.dat -nogui

The index it says is missing is the first index of the ligand (85459). 

I tried this solution (link) but it didn't work for me.

The files for recreating the error will be sent separately.

Thanks again,

Christian

[Christian Meyer]

A google drive link to all the files to recreate error:

https://drive.google.com/file/d/127YwT0ecg1Gm6XVGTnzeNn0CjYi5wQYP/view?usp=drive_link

Thanks

C

[marioe911116]

It doesn't seem you are using PB optimization for membrane systems. Take a look at this tutorial and let us know if you encounter any issues. 

When working with membrane systems is always good to remove all information regarding lipids, water and ions from the topology and the input structure to ease the processing step. Then, you create a brand new index with this new "clean" structure (Protein+ ligand) and generate the tracjectory for only the complex. Then, proceed to use PB optimization for membranes.

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

[Christian Meyer]

Thanks!

I was able to get it to run and the results seemed reasonable.  
Based on this thread (link), I decided that I don't need to include the membrane since the binding site is >15A from the membranes which saved a lot of memory/computation.

For others reference, here is how I rebuilt the system without the lipids/solvents:

# Variables for index numbers
ligand_group=22  # Assuming 22 is the index for the ligand group, adjust if necessary
protein_group=1  # Assuming 1 is the index for the protein group

#Remove periodic boundary conditions. Fit protein (1), write whole system (0)

echo -e "1\n0" | gmx trjconv -s step7_1.tpr -f "$input_xtc" -o step7_1_pbc.xtc -n -fit rot+trans 

#Make a new index of only the protein ligand

echo -e "$protein_group | $ligand_group\nq" | gmx make_ndx -f  step7_1.tpr -o complex.ndx  

# Generate new tpr file for protein-ligand complex. 24= protein+ligand
echo -e "24" | gmx convert-tpr -s step7_1.tpr -o complex.tpr -n complex.ndx

# Generate new PDB file of protein-ligand complex
echo -e "24" | gmx trjconv -s step7_1.tpr -f step7_1.gro -o complex.pdb -n complex.ndx

# Modify the new PDB to add chain information for the ligand

sed -i 's/UNL  /UNL A/g' complex.pdb
       
# Generate a clean index file for the complex
echo -e "1 | 13\nq" | gmx make_ndx -f complex.tpr -o complex_clean.ndx
       
# Create a new trajectory file for the complex
echo -e "24" | gmx trjconv -s step7_1.tpr -f step7_1_pbc.xtc -o complex.xtc -n complex.ndx

# Run MMPBSA calculations
mpirun -np 10 gmx_MMPBSA -O -i mmpbsa.in -cs complex.tpr -ci complex_clean.ndx -ct complex.xtc -cg 1 13 -cp topol_complex.top -o FINAL_RESULTS_MMPBSA.dat -cr complex.pdb -nogui

A couple of notes:

1) For reasons I don't understand, the reference pdb had to have both the ligand and the protein.  As a result, I had to assign a chain to the ligand using the sed command above

2) I manually edited the CHARMM-GUI topol.top file to omit everything but the ligand and protein and saved it as topol_complex.top (see attached).

3) I used the defaults specified by gmx_MMPBSA --create_input pb

Thanks for the great software!

Blessings,

C

[Christian Meyer]

Mistake in first gmx trjconv command.  Fixed here:

echo -e "1\n0" | gmx trjconv -s step7_1.tpr -f  step7_1.xtc -o step7_1_pbc.xtc -n -fit rot+trans 

[Subha Yegnaswamy]

Thanks for the detailed steps.This is really helpful.

I am also working on Membrane-protein-ligand system created using charmm-gui.

I am not able to open the mmpbsa.in file in your link. Have you taken protein-membrane system input file (mmpbsa.in) or protein-ligand system example input file for mmpbsa calculation.

I think since we are calculating only protein and ligand and have deleted the other systems, would it be sufficient if we use the protein-ligand system?

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

May i am missing something.Please clarify.

Thank you

Subha.

[Subha Yegnaswamy]

Dear All,

I tried running the MMPBSA with the steps mentioned above creating a complex.pdb, tpr and clean index file. Still i am getting the following error. Pls find below and also have the attached input files taken for calculation except the xtc.

when I checked  the  atom number mentioned in the error is present in the pdb as well as in the index. but it is mentioned it is not in the .top file. i have provided all the inputs in the complex_topol.top.But still getting the error.Have also checked earlier threads.Not able to get it.

Could you please let me know how to solve it.

Thank you,

Subha.

gmxMMPBSA) bioproject@bioproject:~/anaconda3/envs$ gmx_MMPBSA -O -i gmxMMPBSA.in -cs complex.pdb -ci complex_clean.ndx -cg 1 13 -ct complex.xtc -cp complex_topol.top -o MERSM_MMPBSA.dat -eo MERSM_MMPBSA.csv
[INFO   ] Starting gmx_MMPBSA v1.6.2
[INFO   ] Command-line
  gmx_MMPBSA -O -i gmxMMPBSA.in -cs complex.pdb -ci complex_clean.ndx -cg 1 13 -ct complex.xtc -cp complex_topol.top -o MERSM_MMPBSA.dat -eo MERSM_MMPBSA.csv

tets
[INFO   ] Checking gmxMMPBSA.in input file...
[INFO   ] Checking gmxMMPBSA.in input file...Done.

[INFO   ] Checking external programs...
[INFO   ] cpptraj found! Using /home/bioproject/anaconda3/envs/gmxMMPBSA/bin/cpptraj
[INFO   ] tleap found! Using /home/bioproject/anaconda3/envs/gmxMMPBSA/bin/tleap
[INFO   ] parmchk2 found! Using /home/bioproject/anaconda3/envs/gmxMMPBSA/bin/parmchk2
[INFO   ] sander found! Using /home/bioproject/anaconda3/envs/gmxMMPBSA/bin/sander
[INFO   ] Using GROMACS version > 5.x.x!
[INFO   ] gmx found! Using /usr/local/gromacs/bin/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 Protein_LIG (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 Receptor: Saving group Protein (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 LIG (13) in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_LIG.pdb
[INFO   ] Checking the structures consistency...
[INFO   ]
[INFO   ] Using topology conversion. Setting radiopt = 0...
[INFO   ] Building Normal Complex Amber topology...
[ERROR  ] MMPBSA_Error

The atom 6979 in the complex index is not found in the topology file. Please check that the files are consistent.

Check the gmx_MMPBSA.log file to report the problem.
  File "/home/bioproject/anaconda3/envs/gmxMMPBSA/bin/gmx_MMPBSA", line 8, in <module>
    sys.exit(gmxmmpbsa())
  File "/home/bioproject/anaconda3/envs/gmxMMPBSA/lib/python3.8/site-packages/GMXMMPBSA/app.py", line 98, in gmxmmpbsa
    app.make_prmtops()
  File "/home/bioproject/anaconda3/envs/gmxMMPBSA/lib/python3.8/site-packages/GMXMMPBSA/main.py", line 682, in make_prmtops
    self.FILES.mutant_receptor_prmtop, self.FILES.mutant_ligand_prmtop) = maketop.buildTopology()
  File "/home/bioproject/anaconda3/envs/gmxMMPBSA/lib/python3.8/site-packages/GMXMMPBSA/make_top.py", line 125, in buildTopology
    tops = self.gmxtop2prmtop()
  File "/home/bioproject/anaconda3/envs/gmxMMPBSA/lib/python3.8/site-packages/GMXMMPBSA/make_top.py", line 555, in gmxtop2prmtop
    com_top = self.cleantop(self.FILES.complex_top, self.indexes['COM']['COM'])
  File "/home/bioproject/anaconda3/envs/gmxMMPBSA/lib/python3.8/site-packages/GMXMMPBSA/make_top.py", line 853, in cleantop
    GMXMMPBSA_ERROR(f'The atom {i} in the {id} index is not found in the topology file. Please check that '
  File "/home/bioproject/anaconda3/envs/gmxMMPBSA/lib/python3.8/site-packages/GMXMMPBSA/exceptions.py", line 171, in __init__
    raise exc('\n\n' + msg + '\n\nCheck the gmx_MMPBSA.log file to report the problem.')
MMPBSA_Error:

The atom 6979 in the complex index is not found in the topology file. Please check that the files are consistent.

Check the gmx_MMPBSA.log file to report the problem.
Exiting. All files have been retained.

[mariosergi...@gmail.com]

Please, attach the itp files 

[Subha Yegnaswamy]

Hello Mario,

I am attaching the itp files i have used for the calculation.Kindly let me know where i am going wrong.

Thank you,

Best Regards,

Subha.

[mariosergi...@gmail.com]

Your tpr file contains a complex dimer+ligand. The ligand had 52 atoms, but the itp files had only 44 atoms. So, probably the error is because you are using the wrong itp files. Please, make sure the ligand itp file is the one used in the complex structure

[Subha Yegnaswamy]

ok.thank you.Shall check it.

Best Regards,

Subha.