GROMACS version:2020beta
GROMACS modification: Yes/No yes
Here post your question Dear Gromacss users
I want to add a new molecule (ligand) to force field directory. I am struggling to understand with the formation of rtp file. consulting Adding a Residue to a Force Field — GROMACS 2021.2 documentation and File formats — GROMACS 2021.2 documentation rtp file format mentions to add atoms information mandatory and other parameters as optional. My question is if I am adding atoms information only from where will it get information for bonds or angles? (I mean which atom is bonded to which in a particular angle or dihedral?
Looking forward to your kind reply, please.
Regards
Sadaf
Dear Sadaf,
as far as I know, bond information is not optional, you should add them by hand. Angles and proper dihedrals are typically auto-generated from the bond information, as well as exclusions. Improper dihedrals and cmap interactions (in the CHARMM force field) should also be specified explicitly in the rtp file.
Hope that helps.
Kind regards,
Andras
Dear @awacha thank you very much for your reply. I have kept only atoms and bond information in rtp file and doing grompp gives me warning for missing angles, dihedrals as below:-
No default Proper Dih. types.
However when I put this information to ffbonded.itp I did not get any such warning. Should not all of this information to be kept in rtp file?
Regards
Sadaf
Dear @sadaf,
As far as I understand you have a molecule with dihedrals and angles which do not have interaction parameters in the standard CHARMM force field distribution.
Typically you put only topology information in the rtp files and grompp will get the numeric values from the itp files, based on the atom types of the participating atoms. In this case you need only to write out the parameters, function type index and the atom types only once, and grompp will pick them up every time it encounters an angle/dihedral consisting of these atom types.
You can also include interaction parameters in the rtp files (and in the topology as well), if these angles/dihedrals occur only infrequently, but this is uncommon practice. The drawback is that you have to write the function type index and the values explicitly for each angle, and of course this does not work for autogenerated angles and dihedrals.
I’d suggest to put the parameters in the itp file. And since using wrong or unreliable parameters are dangerous, double check the source from where you got them, also take care to the correct units (CHARMM uses kcal, GROMACS kJ) and functional form (e.g. 1/2 factor in the harmonic potential in CHARMM or GROMACS).
Kind regards,
Andras
Dear @awacha thank you very much for your kind reply. I want to add my residue as a part of protein which needs a single topology using -merge command of gromacs in this case what do you suggest about setting parameters?
Regards
Sadaf
Dear @sadaf,
could you be a little more specific? I’m only guessing now as of your problem:
-
if you are unsure about the parameters for pdb2gmx:
if you want to add a ligand, you do not really need the -merge switch. All it does is that it puts all chains / molecules in a single [ moleculetype ] entry in the topology file. If the ligand is not covalently bonded to the protein, you do not need to -merge them. If you try it, pdb2gmx will most probably complain about the different residue types (Protein for the protein and Unknown/Other for the ligand, see residuetypes.dat in $GMXTOP).
You can also split the structure into the protein and ligand part, do “pdb2gmx” separately, and re-merge the topology and the coordinate files. -
if the problem lies with the interaction parameters: you need to find a parametrization of the ligand. E.g. if you use the CHARMM force field, CGenFF is your friend. You can convert the CGenFF output stream file to GROMACS format with the Python script found on Alexander MacKerell’s home page.
Hope that helps,
Andras
Dear @awacha
Thank you again for your kind and comprehensive reply. yes, that is exactly what I am interested in, basically, I am adding ligand and protein in a single molecule type to generate a single topology in which I could select few atom indexes belonging to ligand and protein for restraining later on during simulation. For this, I have added ligand in residuetype.dat file and using -merge function to get the topology as mentioned above.
Could you please explain a bit if I do pdb2gmx separately how can I emerge the topology and coordinate files?
Dear @sadaf,
GROMACS handles the topology in the same way as C or C++ files are preprocessed by the compiler, i.e. #include, #define #ifdef etc. directives can be used. #include is especially useful: you can move out any range from a topology file to a separate file, and grompp will see it as if it was still there (for your education, check the help of gmx grompp, and check the file written with the switch “-pp”).
If you do pdb2gmx separately (split the pdb file by hand and apply pdb2gmx on each one), you will end up with two topol.top files (preferably do each part in different directories). Copy out the [ moleculetype ] section of each topology file to a separate file (e.g. protein.itp and ligand.itp). Construct a new topology file. It should contain everything before the [ moleculetype ] sections of the two top files generated by pdb2gmx (lines “#include”-ing the forcefield.itp file from the FF directory and the water/ion topologies). Afterwards you should #include protein.itp and ligand.itp. Finally, make a [ system ] block with the appropriate [ molecules ] section where you give the count of molecules of each molecule type.
You should also save the two conf.gro files written by pdb2gmx: rename one to protein.gro, the other to ligand.gro. Their importance: in contrast to the pdb files you have made by splitting the original one, atoms in these files are in the same order as they appear in the topology (moleculetype section). Atom positions are not changed, so you can simply concatenate the atom information in a text editor, taking care to edit the header (2 lines) and the box information (last line) appropriately.
Read the following thoroughly:
https://manual.gromacs.org/current/reference-manual/topologies/topology-file-formats.html#topology-file
and
https://manual.gromacs.org/current/reference-manual/topologies/topology-file-formats.html#coordinate-file
If you want to restrain atoms, you can generate position restraints by molecule type (gmx genrestr), or you can add [ intermolecular_interactions ] in the [ system ] section of the topology in the same way as in the [ moleculetype ] / [ bond ], [ angle ] etc. blocks, but using global atom indices.
Kind regards,
Andras