Decomposition Analysis of MD Simulation data (in CHARMM Forcefield)

[Sarthak Trivedi SVNIT]

Respected Fellows/ Seniors,

As I have successfully completed MMPBSA calculations for the protein-ligand system, now I want to perform stability calculations and a free binding energy decomposition analysis for my protein-ligand system. As my MD simulation is based on the CHARMM36 force field, I'm perplexed about the changes I need to make to the input file (MMPBSA.in). I would request that you put your valuable guidance to use and let me proceed further.

Thank you in advance. 

Regards,

Sarthak Trivedi. 

[Mario Sergio Valdes]

The stability option is rarely used, since it is implicitly already contained in the binding free energy calculation. It only requires specifying it on the command line (-s or --stability).

For energy decomposition it only requires specifying the namelist &decomp in the input file. Please see this example (https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/examples/Decomposition_analysis/) and the documentation (https://valdes-tresanco-ms.github.io/gmx_MMPBSA/ dev/input_file/#decomp-namelist-variables)

Mario S.

[Sarthak Trivedi SVNIT]

Hello,

I just applied the solution and got the error for the decomposition analysis. I am attaching the error message. I have checked the topol.top file for the extra spaces, everything seems ok in the topol.top file. Kindly guide me.

(base) xyzw:~/Complex_02_50ns_Final_Decomp$ mpirun -np 16 gmx_MMPBSA MPI -O -i mmpbsa.in -cs md_0_10.tpr -ci index.ndx -cg 20 21 -ct md_0_10_center.xtc -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv -do FINAL_DECOMP_MMPBSA.dat -deo FINAL_DECOMP_MMPBSA.csv
[INFO   ] Starting gmx_MMPBSA v1.5.7
[INFO   ] Command-line
  mpirun -np 16 gmx_MMPBSA MPI -O -i mmpbsa.in -cs md_0_10.tpr -ci index.ndx -cg 20 21 -ct md_0_10_center.xtc -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv -do FINAL_DECOMP_MMPBSA.dat -deo FINAL_DECOMP_MMPBSA.csv

[INFO   ] Checking mmpbsa.in input file...
[INFO   ] Checking mmpbsa.in input file...Done.

[INFO   ] Checking external programs...
[INFO   ] cpptraj found! Using /home/drr-18/miniconda3/bin/cpptraj
[INFO   ] tleap found! Using /home/drr-18/miniconda3/bin/tleap
[INFO   ] parmchk2 found! Using /home/drr-18/miniconda3/bin/parmchk2
[INFO   ] sander found! Using /home/drr-18/miniconda3/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   ] Get PDB files from GROMACS structures files...
[INFO   ] Making gmx_MMPBSA index for complex...
[INFO   ] Normal Complex: Saving group 20_21 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 20 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 21 in _GMXMMPBSA_COM_index.ndx file as _GMXMMPBSA_LIG.pdb
[ERROR  ] MMPBSA_Error gmx_MMPBSA does not support water/ions molecules in any structure, but we found 39023 molecules in the complex..
           Check the gmx_MMPBSA.log file to report the problem.
  File "/home/drr-18/miniconda3/bin/gmx_MMPBSA", line 8, in <module>
    sys.exit(gmxmmpbsa())
  File "/home/drr-18/miniconda3/lib/python3.9/site-packages/GMXMMPBSA/app.py", line 98, in gmxmmpbsa
    app.make_prmtops()
  File "/home/drr-18/miniconda3/lib/python3.9/site-packages/GMXMMPBSA/main.py", line 548, in make_prmtops
    self.FILES.mutant_receptor_prmtop, self.FILES.mutant_ligand_prmtop) = maketop.buildTopology()
  File "/home/drr-18/miniconda3/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 116, in buildTopology
    self.gmx2pdb()
  File "/home/drr-18/miniconda3/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 432, in gmx2pdb
    self.check4water()
  File "/home/drr-18/miniconda3/lib/python3.9/site-packages/GMXMMPBSA/make_top.py", line 454, in check4water
    GMXMMPBSA_ERROR(f'gmx_MMPBSA does not support water/ions molecules in any structure, but we found'
  File "/home/drr-18/miniconda3/lib/python3.9/site-packages/GMXMMPBSA/exceptions.py", line 169, in __init__
    raise exc(msg + '\nCheck the gmx_MMPBSA.log file to report the problem.')
MMPBSA_Error: gmx_MMPBSA does not support water/ions molecules in any structure, but we found 39023 molecules in the complex.
Check the gmx_MMPBSA.log file to report the problem.
Error occurred on rank 0.
Exiting. All files have been retained.
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0

STRUCTURE OF INPUT FILE (mmpbsa.in):

Sample input file for decomposition analysis
This input file is meant to show only that gmx_MMPBSA works. Althought,
we tried to used the input files as recommended in the Amber manual,
some parameters have been changed to perform more expensive calculations
in a reasonable amount of time. Feel free to change the parameters
according to what is better for your system.

&general
sys_name="Decomposition",
startframe=1,
endframe=5000,
forcefields="leaprc.protein.ff14SB"
/
&gb
igb=5, saltcon=0.150,
/
#make sure to include at least one residue from both the receptor
#and ligand in the print_res mask of the &decomp section.
#this requirement is automatically fulfilled when using the within keyword.
#http://archive.ambermd.org/201308/0075.html
&decomp
idecomp=2, dec_verbose=3,
print_res="within 5"
/

Thank you in advance.

[Mario Sergio Valdes]

This problem is related to the index groups... Check that you have correctly selected the molecules in the groups. Currently, see that you are including water or ions in one of the groups
mario s

[Sarthak Trivedi SVNIT]

Hello Sir,

Thank you so much. You are right, but my receptor (protein) is in group 1, and the ligand is in group 13. How can I change the default groups (1 and 13). I tried to run the code by changing it from 20 21 to 1 13. But still, the error persists. So how can I proceed further. 

Thank you.

Regards,

Sarthak 

[Mario Sergio Valdes]

Please, send me the files (*.tpr and *.ndx) to check them

[Sarthak Trivedi SVNIT]

Here, i am attaching a google drive link. Please try to open it and let me know if it's opening properly and working fine. The message was too long so I had to upload it on a drive.

Decomposition Analysis

[Mario Sergio Valdes]

Sorry for the huge delay, I really forgot all about it.
According to the files, you should have no problems if you define the groups as -cg 1 13. The groups only contain the molecules of interest so it should not give errors.
If the error still persists, please attach the files to verify (*.tpr, *.ndx, 10 frames *.xtc, *.top + *.itp + *.prm + charmm36-jul2022.ff folder)

Mario S.

[Sarthak Trivedi SVNIT]

Hello Sir,

I appreciate your help, i have updated the google drive folder for you, which i shared with you previously. Please, access it and feel free to ask for further files you require. 

Regards,

Mr. Sarthak Trivedi
Research Fellow,
Department of Physics,
Sardar Vallabhbhai National Institute of Technology, Surat.

--
You received this message because you are subscribed to the Google Groups "gmx_MMPBSA" group.
To unsubscribe from this group and stop receiving emails from it, send an email to gmx_mmpbsa+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/gmx_mmpbsa/462e6a30-fe5e-4ddd-a203-79facc24e3a2n%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Mario Sergio Valdes]

It works great for me, actually. You have to include the -cp topol.top option so that it doesn't end in error since the FF is CHARMM.

The command-line

mpirun -np 8 gmx_MMPBSA -O -i mmpbsa.in -cs md_0_10.tpr -ci index.ndx -cg 1 13 -ct md_fit.xtc -cp topol.top

Run a few short tests to make sure everything is running smoothly. For example, setting endframe = 10 would perform the calculation only for the first 10 frames.

Please, make sure that the trajectory is consistent (check this section https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/gmx_MMPBSA_running/#before-running-gmx_mmpbsa)

Additionally, I think that the ligand has problems. There should be no bond between C2 and C20 due to the (sp2) configuration of both carbons. I think you should check it carefully.

Mario S.

[Sarthak Trivedi SVNIT]

Hello Sir,

Thank you so much for your response and accept my apology for delayed response. I will surely try it and check it and will update you in this regard.

Thank you for your time and consideration.

Regards,

Sarthak

[Sarthak Trivedi SVNIT]

Hello Sir,

As per your suggestion about adding flag -cp to define topology file, i was able to successfully run the command. But i am confused about the results, because my system is prepared with CHARMM Forcefield, The sample input file for Decomposition analysis on the section (Decomposition Analysis https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/examples/Decomposition_analysis/ ) has the structure as follows: 

Sample input file for decomposition analysis
This input file is meant to show only that gmx_MMPBSA works. Althought,
we tried to used the input files as recommended in the Amber manual,
some parameters have been changed to perform more expensive calculations
in a reasonable amount of time. Feel free to change the parameters
according to what is better for your system.

&general
sys_name="Decomposition",
startframe=1,

endframe=10,

forcefields="leaprc.protein.ff14SB"
/
&gb
igb=5, saltcon=0.150,
/
#make sure to include at least one residue from both the receptor
#and ligand in the print_res mask of the &decomp section.
#this requirement is automatically fulfilled when using the within keyword.
#http://archive.ambermd.org/201308/0075.html
&decomp
idecomp=2, dec_verbose=3,

print_res="within 4"
/

My question is, will these results ("leaprc.protein.ff14SB" forcefield) conflict with my CHARMM prepared system ? Because both forcefields are different. Do i need to mention PBRaddi with namelist &pb ? 

As per my understanding, decomposition input file structure is not compatible with CHARMM forcefield, is it correct ?

Please guide me in this regard, i would be pleased and thankful to you for clarification.

Thank you for your time and consideration.

Regards,

Sarthak.

[Mario Sergio Valdes]

My question is, will these results ("leaprc.protein.ff14SB" forcefield) conflict with my CHARMM prepared system ? Because both forcefields are different. Do i need to mention PBRaddi with namelist &pb ? 

No, because when CHARMM is used this parameter is ignored. Please, check this section (https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/howworks/#topology-preparation).

PBRadii is a global parameter, so it is only defined in the &general namelist. Don't matter what method you select, PBRadii is the same for all.

 

As per my understanding, decomposition input file structure is not compatible with CHARMM forcefield, is it correct ?

You can carry out the decomposition for each ff supported in gmx_MMPBSA. Although both PB and GB support decomposition, GB has limitations with less common atoms like Au, Mg, etc.

 

[Sarthak Trivedi SVNIT]

Hello Sir,

Hope you are doing fine. As per your suggestion to check into the ligand structure regarding the statement: 

Additionally, I think that the ligand has problems. There should be no bond between C2 and C20 due to the (sp2) configuration of both carbons. I think you should check it carefully. 

I looked into it, but i want some of your remarks about the ligand structure relaxation, the ligand was optimised using Gaussian 09 ( Using DFT, Basis set: 6-31 + G [d,p] ).  Because during molecular dynamics simulation, when i was preparing ligand parameter file using CGenFF server, it showed bad penalty score for ligand. By literature review and discussions on several forums, it's found that people don't do this type of quantum mechanical optimisation (Using DFT, HF, etc.). Most of the people perform ligand optimisation using molecular mechanics models (Using Avogadro Software etc.) and the get good penalty score during ligand topology preparation (On CGenFF server). Theoretically, Offcourse Density Functional Theory based ligand structure relaxation is considered to be best (As per the discussion with my PhD supervisor, as per theoretical point of view)  So can you comment that these OPTIMISATION is somehow conflicting this ligand problem ? Can you throw the light on the ligand structure optimisation, i mean is it conflicting something with the optimisation?

Thank you in advance. 

Regards,

Sarthak

[Mario Sergio Valdes]

One of the most frequent problems in this type of calculation is an erroneous parameterization of the ligand. It is common to start from defective or incomplete structures. For example, the PDB format is uninformative regarding chemical bonds. This translates into a serious problem since the bond order is defined by the distance between atoms, which is very ambiguous. For these cases, you must define the PDB CONECTs that explicitly indicate the bond. However, the bond order problem persists. In these cases, there are more suitable formats, such as mol2 or sdf. From my experience, optimizing geometry with DFT is not worth it since it consumes a considerable amount of time without the results being alarmingly good. As you saw in your review, Avogadro is widely used because it allows you to edit and relax structure. Keep in mind, that the most important aspect, in this case, is that the bonds and their orders are well-defined.
If you look at the figure, there are several atoms with worrying configurations. For example, as I mentioned before, C2 and C20,  visually and in the topology appear bound. This is completely wrong. The C1 should be an Sp2 carbon, similar to the C6. The C3 seems to have an Sp3 geometric, but with Sp2 valences. My recommendation is that you look for the 2D representation and follow it carefully to assign the bonds correctly. In case you don't have it available, I recommend using other programs to generate the structure, such as ChemAxon, Maestro, OpenBabel, RDKit, etc., and comparing the results.

Note that my opinion is supported only by the ligand visualization. I do not have more information, so there may be elements in which I may be partially wrong.