Inserting a modified residue into a protein/DNA chain using pdb2gmx

If you want to generate hydrogens automatically with gmx pdb2gmx with rtp/hdb files, the naming has to follow a strict scheme: for example, if Gromacs auto-adds 2 hydrogens with a basename H1, they will be named H11 and H12. This becomes problematic though because if you add 1 hydrogen with a basename H11, it will be H11 too, creating a conflict.

For this reason, I had to (a) introduce a renaming scheme and (b) make the basename non-numeric to avoid confusion.

If your original structure contains generically named hydrogens and it causes a conflict with the auto-generated .rtp entry, you can pass -ignh to gmx pdb2gmx to ignore all existing hydrogens and generate them from scratch. In the long run, I believe using the Gromacs convention is better than the typical H1, H2, …, H30 placed at the end of the structure, as the positions and names of the atoms are more intuitive.

Thus, the simple way to use to_rtp procedure for pdb and itp files where initially all the hydrogen atoms have different non-numerical? In my situation, when the PDB files are created with ā€˜wrong’ hydrogen names now I should to change the hydrogen atom names (I have only 40 names) in the PDB to files to new that to_rtp procedure generated by hand?

I think this tutorial should answer your questions:

If you want to preserve your hydrogen numbering, and your .pdb has all the hydrogens attached/named correctly, go for pathway 1.

In turn, if you want to have a replicable setup for fully automated creation of your systems together with small molecules, pathway 2 will do what you need.

In either case, you shouldn’t have to rename your atoms manually.

Hello, @milosz.wieczor,

Thank you for your help. I’m writing to ask for assistance with an issue I’m having while using the gromologist to build an rtp file from itp files generated by the LigParGen server.
The problem lies in the [dihedrals] that are written by function 3 in the itp files, which are related to the Ryckaert-Bellemans dihedral potential. The gromologist doesn’t convert this type of potential to rtp and leaves only the caption [dihedrals].
Could you please help me? I was wondering if there is a way to convert the dihedrals of function 3 from itp files to rtp files using the to_rtp procedure?

Hi @Yel21,

You’re right, the current version explicitly supports proper dihedral types 1 and 9 and impropers 2 and 4. I’d be happy to add support for type 3 dihedrals though - do you mind attaching an example .itp to your post, or sending one to milosz.wieczor@irbbarcelona.org?

Best,
Miłosz

Hi, @milosz.wieczor,

I’ve attached an example of an. itp file.
UNK_16C3F9.top (3.4 KB)

Since there are limitations on the file extensions that can be uploaded to the forum, I’ve changed the file extension from .itp to .top.

Sorry, that doesn’t look like an .itp file to me, rather an .rtf.

I’m sorry, @milosz.wieczor. It was my mistake — I sent you the wrong file. I’ve attached the correct .itp file below.

UNK_16C3F9.top (17.7 KB)

Sorry for the long delay, I was out on holiday.

Thanks! I pushed the modified version to the library’s gitlab, you can download it and install and check if it works.

Adding the dihedral type was actually trivial but I had to change the way in which hydrogen atoms are recognized in OPLS (normally the criterion was ā€œatom type starting with h/Hā€), so now also the .hdb file should be produced correctly; before it would come out empty.

Thank you, @milosz.wieczor! It’s working perfectly for me.

I was wondering if you could help me with gml.calc_gmx_energy. I’m trying to calculate the energy interactions between two groups for a specific part of the trajectory, for example, between frames 100 and 200 or for all 40,000 frames. However, the results are calculated only for the initial part of the trajectory.

Is it possible to modify the rerun.mdp file to achieve this?

Good to hear it worked!

But there should be no problem with longer trajectories, the default behavior is to rerun and read everything. Can you explain what exactly happens?

Sorry, this was a mistake on my part. The gml.calc_gmx_energy works.

1 Like

hi @milosz.wieczor , I wanted to ask that ,while incorporating the modified residue in the main chain, where do we derive the angle bond type and dihedral parameters from?

You can assign by analogy to get them from the standard protein/DNA force field if you think that makes sense (will depend on your case - you should use chemical intuitions for that). Whichever script/server you’re using will give you some guessed parameter values, but there can be compatibility issues at the boundaries - I assume that’s what you’re referring to.

In that case, find_missing_ff_params has the fix_by_analogy option that requires a dictionary of equivalents, as illustrated here:

Thanks!