Protein-ligand simulation

GROMACS version:
GROMACS modification: Yes/No
Here post your question

I want to perform an MD simulation of a protein-ligand complex. I followed the Justin A. Lemkul tutorial for the protein-ligand complex. The tutorial went smoothly. Next, I tried with my complex that is, I have docked a ligand with the protein of interest using autodock 1.5.6

Upon running the pdb2gmx command, it shows an error. : Fatal error:

Atom HN1 in residue PHE 3 was not found in rtp entry PHE with 22 atoms
while sorting atoms.

I understand that the gromacs is not able to recognize the file due to atom numbering. But how to solve this issue, as I don’t know what forcefield will be compatible with this file.

Kindly help me in this regard.

Hi Sanah,

I never used autodock, so I can’t comment on any incompatibilities with Gromacs. However, did it change the naming of atoms in your output PDB file? I’d go into your autodock-output PDB file and see if the atom naming is consistent, or at least consistent with your input PDB file (the PDB obtained directly from the PDB homepage should compatible with pdb2gmx, at least from experience). If it’s not consistent, try to change the “NH1” for “H” (I think that’s the default name for backbone H of the N-H).
What force field did you plan to use initially?

Best,
Jacek

And one more thing:
Try to use the the -ignh option in pdb2gmx (see http://manual.gromacs.org/documentation/2018/onlinehelp/gmx-pdb2gmx.html). Since this causes to avoid all hydrogens, it may help potentially.

Best,
Jacek

Yes the atom types have changed from the original pdb.

-ignh worked. Thanks a lot!

@Jacek I got a new error.
Since i have the protein-ligand docked complex. i followed with the comment. as you said -ignh, next it gives this error.
Fatal error:
Residue ‘LIG’ not found in residue topology database

So, gromacs has a collection of parameters for all aminoacids, nucleotides and other typical biomolecules. However, your ligand may not be implemented in Gromacs’ parameter library or the naming is off (i.e. 3-letter-code and atom names must match something in the gromacs parameter files).

If gromacs does not contain parameters for your LIG, you have to do the following here (very similar to Justin Lemkul’s tutorial): Take you protein-ligand-PDB-file and separate protein and ligand into two separate PDB-files (keep both!). The protein PDB-file can be used in pdb2gmx, now. For the ligand you have to get parameters, either from literature, if available, or you use something like Antechamber within Ambertools (have you done that before?).
After you have your ligand parameters you just follow Justin’s tutorial. In order to use the extracted ligand PDB structure, you can use editconf to turn the PDB-file into a gro-file. Then you just combine “protein.gro” with “ligand.gro” and go ahead as described by Justin.

Let me know if you need more detailed info on certain steps, such as getting parameters using antechamber or so.

1 Like

@Jacek I have separated the protein and ligand as mentioned in the tutorial .
I haven’t used antechamber earlier. please help me with it.

Ok, I think I have somewhere a little manual written. I’ll see if I can find it.

In the meantime, to obtain parameters for your molecules you need:

  1. Gaussian (or alternative DFT/QM program - you have to do a structure optimization)
  2. You have to build AmberTools on your system (Maybe you have it already on the system you are using? You should check the modules that are installed on the cluster/system you use). AmberTools is free: http://ambermd.org/AmberTools.php
  3. Python to translate Amber parameters into Gromacs format.

Also: what force field are you planning to use?

The last point is crucial for force field compatibility. Ligand preparation should be done in a manner consistent with the protein force field. That may require simple use of a general force field but may also require some QM work.

I have amber tools and python in my system.
Currently i am trying to use CHARMM36 all-atom force field (July 2017) as mentioned in the tutorial .

I tried using CHARMM-GUI for the preparation of the input file. I thought it may help in removing the error.
But still the error is same :

Writing topology
Processing chain 2 ‘H’ (33 atoms, 1 residues)
Warning: Starting residue LIG0 in chain not identified as Protein/RNA/DNA.
Problem with chain definition, or missing terminal residues.
This chain does not appear to contain a recognized chain molecule.
If this is incorrect, you can edit residuetypes.dat to modify the behavior.
8 out of 8 lines of specbond.dat converted successfully


Program gmx_mpi, VERSION 5.0.4
Source code file: /root/locuz/gromacs-5.0.4/src/gromacs/gmxpreprocess/resall.c, line: 647

Fatal error:
Residue ‘LIG’ not found in residue topology database
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors

I have gaussian also.

Regarding force field:
Then you should generate CHARMM compatible parameter for your force field using CGenFF https://cgenff.umaryland.edu/. I haven’t used it myself, thought. I’ve used mainly Amber force fields, so probably it’s better if somebody else assists you with the details on that. However, maybe it’s straightforward, you should check the link.

Regarding your error:
Is this an error from pdb2gmx? Then it seems that it still reads your ligand and simply doesn’t know it.
I guess at that point it would be helpful to see your pdb, gro, top, itp files to help you more efficiently.

If you are using CHARMM, follow my tutorial directly. Antechamber and AmberTools would be inappropriate in this case.

I am following your tutorial exactly. after separating the ligand and protein file. Upon using the command pdb2gmx, this error comes up.

Kindly help me with this.

So, I guess there’s something wrong with your protein.pdb file, after you removed the ligand. As said, if you upload your pdb file, it would be easier to help. Otherwise, and just speculation, but just try a different force field (amber99 or so) and see if you have the same error. If you get the same error it should be the pdb file that causes issues

I am not able to upload my .pdb file as it does not support the format

The unknown residue error comes from the ligand, which means you did not separate its coordinates before passing the protein through pdb2gmx.

I understood. I used this command grep LIG clean.pdb > lig.pdb
It created a new file which containes the ligand atoms in lig.pdb
but the clean.pdb is still containing the ligand atoms.
Is it because my protein-complex.pdb (clean.pdb) is having no TER between my protein coordinates and ligand coordinates?
Attaching a few atoms of my pdb

ATOM 4623 OT1 VAL P 303 68.066 -18.170 36.188 1.00 0.00 PROA
ATOM 4624 OT2 VAL P 303 69.530 -16.788 37.135 1.00 0.00 PROA
ATOM 4625 C1 LIG H 0 70.691 5.133 5.607 1.00 0.00 HETA
ATOM 4626 C2 LIG H 0 69.584 4.519 4.822 1.00 0.00 HETA