New residue top

Continuing the discussion from Energy min grompp:

I have converted the itp of my gutathione to rtp format and copied the pdb to my strucuture along with using edit conf to renumber the file with a repacement residue. When I try to convert the pdb I recieve a error

Fatal error:
Residue 175 named GSH of a molecule in the input file was mapped
to an entry in the topology database, but the atom C used in
that entry is not found in the input file. Perhaps your atom
and/or residue naming needs to be fixed.

I checked and residue 175 is a TYR and I cant find any issue with my carbons besides the segmented structure.

The problem is with the GSH residue; the next residue in the sequence expects it to have a C atom to create bonds (e.g. -C in the TYR .rtp entry) and apparently it does not, so pdb2gmx cannot build the expected bonded structure.

OK. I went back and realized the source of my issue was the lack of labeling or missing labels when I convert my pdb to mol2 for cgenff. This issue would result in a str with atom names I dont recognize and cant configure into a proper rtp.

I have since started over confirming my atomnames
I added hydrogens to my residue amine group to satisfy cgenff

.attype warning: monovalent nitrogen not supported;skipped molecule

Then I removed the amine hydrogens in before using the python charmm2gmx script to create an itp (H02 and H12) which I am trying to turn to an rtp.


[ GSH ]
;Gutathionylated cysteine				 3

[ atoms ]
;    atomN      type	 charge     cgnr
        N      NG321 -1.015   1  
       CA      CG311  0.140   2   
        C      CG321 -0.256   3
        O      OG312 -0.923   4
       CB      CG321  0.062   5
       SG      SG301 -0.238   6
      C01      CG321  0.039   7
      N01      NG2S1 -0.395   8
      O01      OG2D1 -0.514   9
      S01      SG301 -0.189  10
      C02      CG311 -0.010  11
      N02      NG2S1 -0.457  12
      O02      OG2D1 -0.516  13
      N03      CG2O1 -1.077  14
      O03      OG2D1 -0.568  15
      O04      OG311 -0.598  16
      C05      CG2O1  0.536  17
      O05      OG2D1 -0.548  18
      C06      CG321 -0.078  19
      O06      OG311 -0.553  20
      C07      CG2O2  0.743  21
      C08      CG321 -0.199  22
      C09      CG321 -0.268  23
      C10      CG311  0.326  24
      C11      CG2O2  0.779  25
      H16      HGA1   0.090  26 
      H21      HGA2   0.090  27
      H01      HGA2   0.090  28
      H02      HGA2   0.090  29
      H03      HGA1   0.090  30
      H04      HGP1   0.333  31
      H05      HGA2   0.090  32  
      H06      HGA2   0.090  33
      H07      HGA2   0.090  34
      H08      HGA2   0.090  35
      H09      HGA2   0.090  36
      H10      HGP1   0.277  37
      H11      HGA2   0.090  38
      H12      HGA2   0.090  39
      H13      HGA2   0.090  40
      H14      HGA2   0.090  41
      H15      HGP1   0.430  42
      H17      HGP1   0.290  43
      H18      HGPAM2 0.396  44
      H19      HGA1   0.090  45
      H20      HGPAM2 0.396  46

do I need to also write the [ bonds] ? I started to convert them from the itp but they dont make sense in that there are repeats and connections not present in the pdb.

I dont fully understand how to read or write hdb so I based it off cys since the backbone matches

GSH        4
1       1       HN      N       CA      -C     
1       5       HA      CA      N       CB      C      
2       6       HB      CB      CA      SG     
1       2       HG1     SG      CB      HB1       

gmx pdb2gmx -f db_2dwp_gsh_.pdb -o 2dwp_gsh.gro

Program: gmx pdb2gmx, version 2022
Source file: src/gromacs/gmxpreprocess/resall.cpp (line 469)

Fatal error:
in .rtp file in residue moleculetype at line:

For more information and tips for troubleshooting, please check the GROMACS
website at Common Errors — GROMACS webpage documentation

This warning suggests something intrinsically wrong with your .mol2 file that should be rectified before even thinking about proceeding. It means the bonded structure is wrong.

I don’t understand this. You should never be making manual edits to coordinates or topologies at this point.

At the bare minimum, yes. Otherwise pdb2gmx has no idea how the molecule is connected and you’ll just get a bunch of atoms in space that aren’t actually bonded in any way. You’ll also need [impropers] if there are any to assign.

I don’t see where this is coming from based on what you’ve provided, but there’s something syntactically wrong with the residue definition.

“attype warning: monovalent nitrogen not supported;skipped molecule”
this was a warning I resolved by adding hydrogens to my nitrogen attached to my CA.
I figured this was because of CGENFF’s intended use on ligands rather than individual members of a polypeptide.
I was trying to match the other rtp files while satisfying CGENFF

I will retry with a new model, the nitrogen in its correct trivalent form bonded to the Calpha and a hydrogen

It concerns me that the resulting mol2 file doesnt display the double bonded oxygens
also can I use atom types to write the bonds of the rtp?
here is the resulting rtp with bonds
`[ GSH ]
;Gutathionylated cysteine 3

[ atoms ]
; atomN type charge cgnr
N NG2D1 -0.498 1
CA CG311 0.104 2
C CG321 -0.279 3
O OG312 -0.923 4
CB CG321 0.012 5
SG SG301 -0.200 6
C01 CG321 0.041 7
N01 NG311 -0.673 8
O01 OG312 -0.873 9
S01 SG311 -0.220 10
C02 CG311 0.073 11
N02 NG311 -0.789 12
O02 OG312 -0.874 13
C03 CG311 -0.049 14
N03 NG321 -1.000 15
O03 OG312 -0.879 16
O04 OG311 -0.609 17
C05 CG311 0.030 18
O05 OG312 -0.878 19
C06 CG321 -0.016 20
O06 OG311 -0.606 21
C07 CG311 0.044 22
C08 CG321 -0.188 23
C09 CG321 -0.217 24
C10 CG311 0.172 25
C11 CG311 0.021 26
H12 HGP1 0.334 27
H17 HGA1 0.090 28
H23 HGA2 0.090 29
H24 HGA1 0.090 30
H25 HGA1 0.090 31
H26 HGA1 0.090 32
H27 HGA1 0.090 33
H01 HGA2 0.090 34
H03 HGA2 0.090 35
H04 HGA1 0.090 36
H05 HGPAM1 0.349 37
H06 HGA2 0.090 38
H07 HGA2 0.090 39
H08 HGA2 0.090 40
H09 HGA2 0.090 41
H10 HGA2 0.090 42
H11 HGPAM1 0.361 43
H13 HGA2 0.090 44
H14 HGA2 0.090 45
H15 HGA2 0.090 46
H16 HGA2 0.090 47
H18 HGP1 0.370 48
H19 HGP1 0.370 49
H20 HGPAM2 0.390 50
H21 HGA1 0.090 51
H22 HGPAM2 0.390 52

[ bonds ]
N CA ; NG321 CG311
N H12 ; CG311 CG321
CA C ; CG311 CG321
CA CB ; CG311 HGP1
CA H17 ; CG321 OG312
C O ; CG321 HGA2
C H23 ; CG321 HGA1
C H16 ; CG321 SG301
CB SG ; CG321 HGA2
CB H09 ; CG321 HGA2
CB H10 ; SG301 SG301
SG S01 ; CG321 SG301
C01 S01 ; CG321 CG311
C01 C02 ; CG321 HGA2
C01 H01 ; CG321 HGA2
C01 H03 ; NG2S1 CG311
N01 C02 ; NG2S1 CG2O1
N01 C03 ; NG2S1 HGP1
N01 H05 ; OG2D1 CG2O1
O01 C03 ; CG311 CG2O1
C02 C05 ; CG311 HGA1
C02 H04 ; NG2S1 CG2O1
N02 C05 ; NG2S1 CG321
N02 C06 ; NG2S1 HGP1
N02 H11 ; OG2D1 CG2O1
O02 C05 ; CG2O1 CG321
C03 C08 ; NG321 CG311
C03 H24 ; NG321 HGA1
N03 C10 ; NG321 HGPAM2
N03 H20 ; OG2D1 CG2O2
N03 H22 ; OG311 CG2O2
O03 C07 ; OG311 HGP1
O04 C07 ; OG2D1 CG2O2
O04 H18 ; CG321 CG2O2
C05 H25 ; CG321 HGA2
O05 C11 ; CG321 HGA2
C06 C07 ; OG311 CG2O2
C06 H13 ; OG311 HGPAM2
C06 H14 ; CG321 CG321
O06 C11 ; CG321 HGA2
O06 H19 ; CG321 HGA2
C07 H26 ; CG321 CG311
C08 C09 ; CG321 HGA2
C08 H07 ; CG321 HGA2
C08 H08 ; CG311 CG2O2
C09 C10 ; CG311 HGPAM2
C09 H06
C09 H15
C10 C11
C10 H21
C11 H27 `
I dont know if the charge#s are correct. I copied them from the str, but they dont match the patterns of the residues in the charmmff rtp.

I ran through pdb2gmx

Program:     gmx pdb2gmx, version 2022
Source file: src/gromacs/gmxpreprocess/add_par.cpp (line 165)

Fatal error:
Atom HN not found in rtp database in residue GSH, it looks a bit like H12

For more information and tips for troubleshooting, please check the GROMACS
website at

I took this as gmx expecting an HN for my backbone amine with the name HN so I renamed my H12 HN
and re ran pdb2gmx

Opening force field file ./charmm36-jul2021.ff/aminoacids.arn

Checking for duplicate atoms....

Generating any missing hydrogen atoms and/or adding termini.
Segmentation fault (core dumped)

thanks for your support through this

going over the errors

Program:     gmx pdb2gmx, version 2022
Source file: src/gromacs/gmxpreprocess/add_par.cpp (line 165)

Fatal error:
Atom HN not found in rtp database in residue GSH, it looks a bit like H12

it seems like there is a process of confirming residues by locating the expected HN which my residue lacks.

I changed H12 to HN in my rtp and pdb and reran

Program:     gmx pdb2gmx, version 2022
Source file: src/gromacs/gmxpreprocess/pdb2gmx.cpp (line 873)

Fatal error:
Atom H in residue GSH 205 was not found in rtp entry GSH with 52 atoms
while sorting atoms.

For a hydrogen, this can be a different protonation state, or it
might have had a different number in the PDB file and was rebuilt
(it might for instance have been H3, and we only expected H1 & H2).
Note that hydrogens might have been added to the entry for the N-terminus.
Remove this hydrogen or choose a different protonation state to solve it.
Option -ignh will ignore all hydrogens in the input.

For more information and tips for troubleshooting, please check the GROMACS
website at

i added ingh option with the HN substitution and

Program:     gmx pdb2gmx, version 2022
Source file: src/gromacs/gmxpreprocess/add_par.cpp (line 165)

Fatal error:
Atom HA not found in rtp database in residue GSH, it looks a bit like HN

For more information and tips for troubleshooting, please check the GROMACS
website at

How can I resolve my HA HN and any other hydrogen bonding issue? are there specifed names for the backbone atoms of all residues that I need to adhere to when designing my GSH?

Your molecule is essentially a polypeptide with an unusual linkage. That’s all glutathione is. I would name everything like a standard amino acid, using all the same nomenclature. Honestly, you can assign all the atom types and charges directly from existing residues rather than going through all this mess.

How would you do this? is there an easier way than transcribing all the cys and gly and Glu bonds first then combining them in no specific order and pulling the s-s from my gsh.
Ive never written an rtp file and I dont know whats wrong with my current one.

Im currently running pdb2gmx and correcting every hydrogen name it spits out as “not found in rtp”

Program:     gmx pdb2gmx, version 2022
Source file: src/gromacs/gmxpreprocess/add_par.cpp (line 165)

Fatal error:
Atom HG1 not found in rtp database in residue GSH, it looks a bit like HN

For more information and tips for troubleshooting, please check the GROMACS
website at

why is gmx trying to build cys from GSH HG1 is the hydrogen that caps the sulfur of cys but doesnt exist in my GSH

Because your .hdb entry tells it to. But that entry is poorly constructed and shouldn’t be used because the atom names are not correct and it is incomplete. Pattern the Cys portion of GSH after the disulfide form of Cys in the force field. I would suggest leaving the other Cys moiety of the disulfide off of this molecule so it is treated as a normal part of the protein, and add on the GSH as a new entry in specbond.dat. It will make your life far easier.

I am following your instrucution.

I am renaming the atom names and types from another entry with a disulfide bond.
I deleted my hbd and am writing a new one

Why is it that gmx cant read the hydrogens in my rtp or pdb?

I rewrote my specbond.dat

CYS	SG	1	HEM 	CAC	1	0.18   CYS2	HEME
HIS	       NE2	1	HEM          FE       1	0.2	HIS1	HEME
MET	SD	1	HEM          FE	     1	0.24  MET	HEME
CO          C       1    HEME         FE      1       0.19    CO      HEME
CYM      SG      1     CYM          SG      1       0.2     CYS2    CYS2
GSH      SG      1	GSH          SD      1	      0.2	GSH	GSH

the formatting is correct in my file

I wrote a better hbd and edited my rtp
here is my .gro

I pushed it through pdb2gmx -ingh as to follow the hbd i wrote instead and it looks properly assembled.
First off thank you for your help!
my only question is how would you confirm the assembly?

grompp finishes with no warnings and an energy minimization finishes with a sensible geometry.

since completing my .gro i also added to the ffbonded.itp based on the gsh.prm and itp data.
I have run into errors /w trying to add ion’s with grompp

ERROR 1 [file, line 36837]:
  No default U-B types

ERROR 2 [file, line 36838]:
  No default U-B types

ERROR 3 [file, line 51480]:
  No default Proper Dih. types

ERROR 4 [file, line 51481]:
  No default Proper Dih. types

ERROR 5 [file, line 51482]:
  No default Proper Dih. types

ERROR 6 [file, line 51483]:
  No default Proper Dih. types

ERROR 7 [file, line 51484]:
  No default Proper Dih. types

ERROR 8 [file, line 51485]:
  No default Proper Dih. types

ERROR 9 [file, line 51486]:
  No default Proper Dih. types

ERROR 10 [file, line 63449]:
  Unknown cmap torsion between atoms 2772 2814 2816 2824 2826

There was 1 note

Program:     gmx grompp, version 2022
Source file: src/gromacs/gmxpreprocess/toppush.cpp (line 1715)

Fatal error:
There were 10 errors in input file(s)

For more information and tips for troubleshooting, please check the GROMACS
website at

I added my dihedrals from the created prm ( charmm2gmx) along with my bonds. How would you address these?
also is there a helpful guide for adding to the ffnonbonded.itp?

I looked over another discussion

and reading the info of the errors the top is listing contacts between my new residue and the backbone atomistic contacts that are missing description that should be covered in the dih types of the ffbonded.itp
do i use similar entries to describe these dihedrals?
for example

ERROR 4 [file, line 51483]:
  No default Proper Dih. types

I see in the topology file ln 51483
2748 2766 2768 2770 9
which refers to the solv.gro
204LYS N 2746 10.735 7.696 4.086
204LYS HZ2 2764 10.981 7.390 4.350
204LYS C 2766 10.837 7.764 3.868
205GSH CA 2770 10.763 7.793 3.652

so i would look up another instance where lys is next to cys and add it to the ffbonded.itp?
otherwise build a 3 peptide chain of the residues surrounding gsh and push it through cgen to copy the missing dihedrals?

Where did you add them? What are the lines you added? It is sufficient to #include a topology file that lists any new bonded parameters rather than necessarily hack ffbonded.itp.

Thank you. I thought that wouldnt work for residues
the ligand topology from charmm doesnt match the rtp I wrote. I dont think I can just copy it over to my rtp because the backbone doesnt match

would I be able to just add the itp and prm for the r-group since the bonded and nonbonded parameters of the backbone match cys (and should be already included) or do the parameters need to be specified.

should I rewrite the rtp (mainly the r-group, keeping the cys-based backbone) and rerun pdb2gmx with the renamed retyped rtf and hbd (based on the charmmgui topology)?
could I just make the addition to the R-group a separate residue in the rtf and the hbd, break up the construction into 2 parts?

Could I just change my cys to a disulfide cys2 and add my gutathione as a ligand?
how would I describe this in my rtp and hbd? could I just write gutathione with a incomplete sulfur?

ultimately i dont understand how to add my new residue topology parameters and coordninates to the protein.
Also should I create a new residue to add in place of cys to insert into the protein or just a new R group to add as a covalent attachment ?
currently I understand my errors are due to my bonded and nonbonded information missing from my top file in the form of prm and itp files

I went forward and wrote a new pdb with cys2 in place of ptm site
i followed your protein ligand instructions in designing my gutathione as a ligand and using cgenff to parameterize it.
after parameterization I removed the backbone atoms and converted it to a .gro
i pasted it in to make my complex.gro
Im having trouble customizing the itp to only represent the R group of my original structure.
also how do i write this as covalently attached to cys2?
I rewrote my specbonds.dat but how do i call for that to be applied in this .top

ERROR 2 [file gsh3.itp, line 52]:
  Atom index (37) in bonds out of bounds (1-34).
  This probably means that you have inserted topology section "bonds"
  in a part belonging to a different molecule than you intended to.
  In that case move the "bonds" section to the right molecule.

currently I am correcting the bonds to only connect the remaining r group atoms

[ bonds ]
;	ai	  aj funct			  c0			c1			  c2			c3
    13    34	 1 ;       CG321       HGA2  **from atom 13 to 34 single bond?**
    5    28	 1 ;       CG321       HGA2
    7     9	 1 ;       CG311      NG2S1
    7    18	 1 ;       CG311      CG2O1
    7    26	 1 ;       CG311      CG321
    7    37	 1 ;       CG311       HGA1
    8    26	 1 ;       SG301      CG321
    9    10	 1 ;       NG2S1      CG2O1
    9    38	 1 ;       NG2S1       HGP1
   10    11	 1 ;       CG2O1      CG321
   10    12	 1 ;       CG2O1      OG2D1
   11    13	 1 ;       CG321      CG321
   11    41	 1 ;       CG321       HGA2
   11    42	 1 ;       CG321       HGA2
   13    14	 1 ;       CG321      CG311
   13    43	 1 ;       CG321       HGA2
   13    44	 1 ;       CG321       HGA2
   14    15	 1 ;       CG311      NG321

If my understanding is correct (above) how do I connect the SG301 of gsh to the SM of cys2?

Correct me if I’m wrong. I currently see 2 ways to build my structure. The recommended route of :

building my alternative residue pdb,

Formatting to determine the topology (cgenff),

Then formatting the structure based on the topologies already provided in the gmx rtp (identify the names and types then replace the names types and bond information of the itp to match a traditional residue, for me I could see this working by replacing the backbone of gsh with atoms matching cys)
paste the formatted itp in the rtp
write up an hbd for new residue
rewrite the .gro adding the coordinates of the gsh r group to the cys and renaming it GSH

the other way which I am currently working is to
change cys to cys2 with an available sulfur
treat the R-group as a ligand
convert the pdb to gro add it to the protein .gro
build the pdb, format it to push through cgenff then remove the backbone
change the sulfur to something indicating its attachment to the cys2 sulfur atom type sm in the ligand itp
edit the $gmx top specbonds.dat to allow the bond.

also could you help me understand the [pairs] section of the itp

[ pairs ]
;	ai	  aj funct			  c0			c1			  c2			c3
    1     4	 1
    1    30	 1
    1    31	 1
    1     6	 1
    1    35	 1
    1    36	 1
    2    32	 1
    2     8	 1

Is there a better way than reading through the itp and looking up each line of atoms bonds and dihedrals to rewrite it and then rewrite it again when the numbers refering to the codes change?

The simplest thing to do here is to simply create an all-in-one residue. Start with the disulfide Cys and build out from the backbone by adding the γ-Glu-Cys(disulfide form)-Gly. Pretty much all the bonded connectivity and associated .hdb information comes directly from the existing residues. The only difference is the γ-peptide bond, but you can take all the charges and atom types from the backbone of existing residues. There may be a missing dihedral or two, but those are simple to introduce by analogy.