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.
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?
@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.
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:
Gaussian (or alternative DFT/QM program - you have to do a structure optimization)
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
Python to translate Amber parameters into Gromacs format.
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 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
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.
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 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