<< first time running with GROMACS >>

[I. Camps]

Hello,

I'm working with collaborators that had done a molecular dynamic with GROMACS 2019.3 and left for me to calculate the binding energy.

I had never worked with GROMACS and I am not familiar with the specificity of the different file types.

In order to accomplish my task, I successfully installed gmx_MMPBSA but now I am in trouble setting up the calculations. Obviously, I went through the tutorials, examples, etc.

My problem now is to obtain some files needed to run gmx_MMPBSA:

- A Reference Structure file

- The Antechamber output mol2 file of ligand parametrization

The calculations where done using CHARMM36 force fields.

Is there any place where I can read about how to get such files?

Thanks in advance.

Camps

[Mario Sergio Valdes]

Hi Ihosvany and welcome!

We are very happy to have you here.

Let's see...

First of all, since it is a parameterized system with CHARMM ff then the antechamber output is not necessary. This is because in this case the topology is converted from GROMACS to Amber. The ligand mol2 file is only needed when the topology is built from the structure extracted from the dynamics using AMBER ff. You can check this here (https://valdes-tresanco-ms.github.io/gmx_MMPBSA/howworks/) and a brief explication about the files types (https://valdes-tresanco-ms.github.io/gmx_MMPBSA/howworks/#required-input-files)

The reference structure is extracted from the dynamics itself or the tpr. This structure is not mandatory, it is only used to correctly assign the chain IDs. In some cases, it must be defined to avoid ambiguities (chains merged in GROMACS, structure extracted from tpr of old versions of GROMACS, the use of gro files which do not have string identifiers, etc). gmx_MMPBSA can assign it automatically, but we recommend that the user check the structure beforehand to make sure it is consistent. In most cases (as in the examples), it is unnecessary to use it, since they are systems without greater complexity or the extracted structure (from a tpr or PDB defined with the -cs option) in the process is sufficient.

Extract the pdb file (reference structure) from a tpr:

gmx editconf -f com.tpr -o com.pdb -n index.ndx

select: receptor + ligand group

from dynamics:

gmx trjconv -f md.xtc -s com.tpr -n index.ndx -b 0 -e 0 -o com.pdb

or

gmx trjconv -f md.xtc -s com.tpr -n index.ndx -dump 0 -o com.pdb

In this particular case, you can use the same command line as this example (https://valdes-tresanco-ms.github.io/gmx_MMPBSA/examples/Protein_ligand_CHARMMff/)

We have plans to extend gmx_MMPBSA for NAMD and we have already reviewed most of the requirements to do so. I think we will be able to implement it soon and start testing it. At the moment we are about to send the article ...

Any questions or difficulties we will gladly help you

HTH!

Mario S.

[I. Camps]

Hello Mario,

Thank you very much for your help!

I could run my test using your tips. But, the calculation stop with the error:

CalcError /home/icamps/miniconda3/envs/gmx_MMPBSA/bin/sander failed with prmtop COM.prmtop!
If you are using sander and PB calculation, check the *.mdout files to get the sander error and check this thread (http://archive.ambermd.org/201303/0548.html).
Check the gmx_MMPBSA.log file to report the problem.

I am attaching the full log file here.

Thank you in advance.

Camps

PS: Good to know you are porting gmx_MMPBSA to work with NAMD files :)

[Mario Sergio Valdes]

This error usually occurs because sander (the program that carries out the calculation) fails with some component of the topology (It cannot assign the PBRadii to an atom, there are inconsistencies in the topology, etc.)

Does the topology have any unusual atoms? (Metals, dummy or pseudo atom, etc. CHARMM has a Lone Pair (LP) atom that is not supported)

Did you define the option radiopt = 0? This option causes the PBRadii to be used from the topology (prmtop) instead of estimating it with sander

Please check the end of the * .mdout files, there sander shows the error that occurred and shows us what it is specifically about

HTH!

Mario S.

[I. Camps]

I took a look at the topology files and didn't find anything unusual. But, as I am not a GROMACS user, maybe my untrained eyes are missing something, so I decided to attach them here.

About the radiopt, yes, I'm used as recommended: radiopt = 0. The full mmpbsa.in file is:

&general
startframe=1, endframe=11, verbose=2,
/
&pb
# radiopt=0 is recommended which means using radii from the prmtop file
# for both the PB calculation and for the NP calculation

istrng=0.15, fillratio=4.0, radiopt=0
/

The only message in the *.mdout files is:

           Error opening unit    5: File "_GMXMMPBSA_pb.mdin" is missing or unreadable

but the file is there with the content:

 File generated by gmx_MMPBSA
 &cntrl
  ntb=0, nsnb=99999, cut=999.0, imin=5,
  igb=10, ipb=2, dec_verbose=0,
 /
 &pb
  istrng=150.0, radiopt=0, solvopt=2,
  maxitn=1000, fillratio=4.0,
 /

Thanks.

Camps

[Mario Sergio Valdes]

Hi Ihosvany

I think that it is a problem with the working directory location. To facilitate identification, if possible send us a sample of the files (topology + itps or parm folder, tpr, xtc [only 10 frames], index, groups)

In this way, we can carefully review each step and provide you with better advice

Best!

Mario S.

[marioe911116]

Hi Ihosvany! 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/

don't hesitate to contact us if you have any questions...

best!

Mario E.