My apologies for the previous locked topic, the topic was created in the midst of drafting the question.
GROMACS version: 2020
GROMACS modification: No
My goal is to run a simulation on a protein with a ligand covalent bonded to cysteine using CHARMM36. I’ll be getting reliable parameters for the cys-ligand from a CHARMM expert collaborator so obtaining parameters is not an issue.
From reading the documentation and various forum discussions, my understanding is that I have to do the following:
My cys-ligand will be handed over as a .str which I can convert to .itp and .prm.
Add this as a new residue to my merged.rtp file in the charmm36.ff, this will take some manual editing but the .rtp is nearly identical to the .itp so this shouldn’t be too difficult. Does this have to be in alphabetical order or can I just stick it on the end?
Add any new atom types to atomtypes.dat (I don’t think this applies in my instance)
Copy any new bondtypes, angletypes, or dihedraltypes to lines to the corrosponding section in ffbonded.itp
Add a line for my residue to residuetypes.dat
Then I should set up my pdb with the new residue named with the corresponding new residue name and run pdb2gmx to generate the system topology.
Is this the right approach for this? Is there anything I am missing?
In my .top file, the bond information is listed as atom index numbers but all the bond information in the residue database are listed as atom names. Is there some way to convert this .top to an .rtp or .itp format to make it easier to add to this custom residue to merged.rtp? Or can index numbers be used in the .rtp instead of atom names?
The .rtp format requires atom names. The conversion should be straightforward - establish a hash (Perl), dictionary (Python), etc. of atom names and numbers, and print out the range of values in the bond list that correspond to your residue.
I wrote a python script for converting formats from top2rtp that worked out very well. Thank you for the suggestion.
I’ve added the atoms and bonds for my new ‘residue’ to merged.rtp and added it to residuetypes.dat. My custom residue is Cys covalently bound to a peptidomimetic inhibitor, thus had atom with identical atom names that required renaming (atom types did not change). This mostly involved adding a number to the end of the atom name, eg. the second N became N2 but the atom type remained NH1. No atom names for the Cys part were changed.
pdb2gmx creates a topology with no errors raised. However, upon closer inspection, I see that no peptide bond gets created between the backbone carbonyl C of my new residue and the amine N of the next residue. A dihedral and cmap gets created though. The distance between the C of my new residue and the N of the next residue is 1.3 angstroms. The peptide bond with the preceding residue gets created with no issue.
How does pdb2gmx know when to generate interesidue backbone bonds?
Why might such an issue occur and how to remedy it?
I clipped this section out of my protein to run the topology generating tests on a smaller peptide. Below is the section of the topology concerning the C (43) atom in question where atom 114 is the N that should be bonded to it.
EDIT:
I think I might have found my problem.
[ bonds ] requires a ‘+’ or ‘-’ to tell pdb2gmx to connect backbone atoms.
When I added
C +N
to my merged.rtp entry for my residue, pdb2gmx was able to create the necessary bond.
I don’t see this detail in the documentation, it might be good idea describe this maybe in the .rtp file formats section.
My follow up question is this: Why does CYS use C +N while the example in the documentation (GLY) use -C N?
The second number is a function type, which dictates the geometry of the constructed H atoms. The PDF manual has a full description of the functions and how the reference atoms are used.
I suppose when adding [ dihedraltypes ] to ffbonded.itp, the order of atoms does not make a difference. E.g. a dedral between atoms “c2 cc cc ha” is the same as reversed one “ha cc cc c2”, right?
I see there are two sections [ dihedraltypes ] in ffbonded.itp, one with func 4 and one with func 9. My dihedrals are func 3; should I make a new [ dihedraltypes ] or can I just append it to either of the above. Or should I keep them in the rtp file?
Having added my improper dihedrals to the rtp file, where do I add the associated values? Directly into the rtp (i.e. adding [ dihedral section ]?) or into the ffbonded.itp?
Ok, so I figured out my issue with the improper dihedrals.
However, I’m still not sure I’m understanding correctly how to deal with my proper dihedrals. Using antechamber I got dihedrals of type 3 (Ryckaert-Bellemans), but it seems that in ffbonded.itp or aminoacids.rtp only type 9 dihedrals (multiple proper dihedrals) are used. I suppose I have to translate my dihedrals to type 9 (and use the #define_torsion… to use multiples) or is there another way to deal with that?
You shouldn’t be parametrizing your species using antechamber if you’re using the CHARMM force field. The conventions used by AMBER and CHARMM are quite different, and indeed you can’t directly use any R-B dihedrals. You need to use proper dihedrals; the only somewhat unique thing about CHARMM is that it allows more than one term per torsion, unlike some other force fields, hence the difference in function type. I have no idea what #define_torsion is, sorry.
thanks for your response.
I should have mentioned that I am trying to do this for AMBER. Sorry for the missing info.
Ok, then I will try to add some lines to acpype.py to obtaining proper dihedrals - if I’m not mistaken, antechamber provides proper dihedrals which are then translated by acpype.
Thanks for the help, I think I’m almost there.
Best,
Jacek
BTW: #define_torsion_… lines can be found in amber99SBildn for specific amino acid residues, which (I think) overwrite the defaults for certain atomtypes if they don’t match those defaults.
OK, then ignore what I said above about CHARMM. You were posting on a thread about a CHARMM topology so I assumed it was a related issue. Better to start a new thread if it’s a new/unrelated problem.
I wanted to define a ligand covalently bonded to an Arginine, I defined a ligand attached Arg by CHARMM website and have ipt and prm files, but I don’t know how to convert itp file to the rtp format for adding to the merged.rtp file. It seems you wrote Python code to convert itp to rpt format. I am not good at scripting so I wanted to know if you could help me with this conversion.
I really appreciate any help you can provide.
I am attempting to add a covalent ligand to an existing force field in GROMACS, and I’m facing some difficulties in the process. I have already explored the available force fields, but none of them seem to have the specific parameters for my ligand. So, I believe the best approach would be to modify an existing force field to incorporate the necessary parameters for my covalent ligand. However, I’m uncertain about the exact steps to follow. Has anyone successfully accomplished a similar task before? Thanks. found information