Pdb2gmx not able to form bonds in an unusual ligand

Ligand1 · January 27, 2025, 7:30pm

GROMACS version: 2022.5
GROMACS modification: No

Hi community members,

I am sorry but this is going to be a long post since I am going to be very detailed (I’ll try my best). I am working on a building the correct topology for an unusual ligand which is two glutathione molecules with a cadmium in the center. Here is how it looks:

To begin with, I tried to make the ligand topology by the CGenFF server. Since CGenFF does not take the Cadmium, I manually removed the cadmium from the mol2 file that I supplied. I also removed the -H from the sulfurs of the Glutathiones. Next, CGenFF made the structure for me without the cadmium and I downloaded the Gromacs ready files from the server. The downloaded package included a Charmm36 forcefield directory, one pdb file and one top file. Name of the ligand is “ful”. Then based on the coordinates of Cadmium (that I already know), I added the CAD manually in the pdb file. I am pasting the appropriate part of the pdb file below. I have used CAD for Cadmium since the atomtypes.atp file in the Charmm36 forcefield directory provided by CGenFF has it like this:
" CAD 112.41100 ; cadmium (II) cation".

I added an additional line in the pdb file as follows:
ATOM 89 CAD ful 1 116.210 116.207 104.438 1.00324.75 XXXX
TER
END

Next, I also added the CAD atom to the topology file provided by the CGenFF server. The file is pasted below:

[ atoms ]
; nr ; residue 1 ful rtp ful q -2.0
1 NG321 2 CG311 3 CG311 4 OG311 5 OG311 6 CG321 7 CG321 8 CG311 9 OG311 10 NG311 11 CG311 12 CG311 13 OG311 14 CG323 15 SG302 16 NG311 17 CG321 18 CG311 19 OG311 20 OG311 21 NG321 22 CG311 23 CG311 24 OG311 25 OG311 26 CG321 27 CG321 28 CG311 29 OG311 30 NG311 31 CG311 32 CG311 33 OG311 34 CG323 35 SG302 36 NG311 37 CG321 38 CG311 39 OG311 40 OG311 41 HGPAM2 42 HGPAM2 43 HGA1 44 HGA1 45 HGP1 46 HGP1 47 HGA2 48 HGA2 49 HGA2 50 HGA2 51 HGA1 52 HGP1 53 HGPAM1 54 HGA1 55 HGA1 56 HGP1 57 HGA2 58 HGA2 59 HGP1 60 HGPAM1 61 HGA2 62 HGA2 63 HGA1 64 HGP1 65 HGP1 66 HGPAM2 67 HGPAM2 68 HGA1 69 HGA1 70 HGP1 71 HGP1 72 HGA2 73 HGA2 74 HGA2 75 HGA2 76 HGA1 77 HGP1 78 HGPAM1 79 HGA1 80 HGA1 81 HGP1 82 HGA2 83 HGA2 84 HGP1 85 HGPAM1 86 HGA2 87 HGA2 88 HGA1 89 type resnr residue atom cgnr charge mass typeB chargeB massB
1 ful N1 1 -1 14.007 ; qtot -1
1 ful C2 2 0.172 12.011 ; qtot 0.172
1 ful C3 3 0.351 12.011 ; qtot 0.351
1 ful O4 4 -0.636 15.9994 ; qtot -0.636
1 ful O5 5 -0.636 15.9994 ; qtot -0.636
1 ful C6 6 -0.217 12.011 ; qtot -0.217
1 ful C7 7 -0.19 12.011 ; qtot -0.19
1 ful C8 8 0.266 12.011 ; qtot 0.266
1 ful O9 9 -0.554 15.9994 ; qtot -0.554
1 ful N10 10 -0.707 14.007 ; qtot -0.707
1 ful C11 11 0.099 12.011 ; qtot 0.099
1 ful C12 12 0.343 12.011 ; qtot 0.343
1 ful O13 13 -0.553 15.9994 ; qtot -0.553
1 ful C14 14 -0.408 12.011 ; qtot -0.408
1 ful S15 15 -0.8 32.06 ; qtot -0.8
1 ful N16 16 -0.824 14.007 ; qtot -0.824
1 ful C17 17 -0.018 12.011 ; qtot -0.018
1 ful C18 18 0.378 12.011 ; qtot 0.378
1 ful O19 19 -0.639 15.9994 ; qtot -0.639
1 ful O20 20 -0.639 15.9994 ; qtot -0.639
1 ful N21 21 -1 14.007 ; qtot -1
1 ful C22 22 0.172 12.011 ; qtot 0.172
1 ful C23 23 0.351 12.011 ; qtot 0.351
1 ful O24 24 -0.636 15.9994 ; qtot -0.636
1 ful O25 25 -0.636 15.9994 ; qtot -0.636
1 ful C26 26 -0.217 12.011 ; qtot -0.217
1 ful C27 27 -0.19 12.011 ; qtot -0.19
1 ful C28 28 0.266 12.011 ; qtot 0.266
1 ful O29 29 -0.554 15.9994 ; qtot -0.554
1 ful N30 30 -0.707 14.007 ; qtot -0.707
1 ful C31 31 0.099 12.011 ; qtot 0.099
1 ful C32 32 0.343 12.011 ; qtot 0.343
1 ful O33 33 -0.553 15.9994 ; qtot -0.553
1 ful C34 34 -0.408 12.011 ; qtot -0.408
1 ful S35 35 -0.8 32.06 ; qtot -0.8
1 ful N36 36 -0.824 14.007 ; qtot -0.824
1 ful C37 37 -0.018 12.011 ; qtot -0.018
1 ful C38 38 0.378 12.011 ; qtot 0.378
1 ful O39 39 -0.639 15.9994 ; qtot -0.639
1 ful O40 40 -0.639 15.9994 ; qtot -0.639
1 ful H41 41 0.39 1.008 ; qtot 0.39
1 ful H42 42 0.39 1.008 ; qtot 0.39
1 ful H43 43 0.09 1.008 ; qtot 0.09
1 ful H44 44 0.09 1.008 ; qtot 0.09
1 ful H45 45 0.414 1.008 ; qtot 0.414
1 ful H46 46 0.414 1.008 ; qtot 0.414
1 ful H47 47 0.09 1.008 ; qtot 0.09
1 ful H48 48 0.09 1.008 ; qtot 0.09
1 ful H49 49 0.09 1.008 ; qtot 0.09
1 ful H50 50 0.09 1.008 ; qtot 0.09
1 ful H51 51 0.09 1.008 ; qtot 0.09
1 ful H52 52 0.403 1.008 ; qtot 0.403
1 ful H53 53 0.349 1.008 ; qtot 0.349
1 ful H54 54 0.09 1.008 ; qtot 0.09
1 ful H55 55 0.09 1.008 ; qtot 0.09
1 ful H56 56 0.403 1.008 ; qtot 0.403
1 ful H57 57 0.09 1.008 ; qtot 0.09
1 ful H58 58 0.09 1.008 ; qtot 0.09
1 ful H59 59 0.414 1.008 ; qtot 0.414
1 ful H60 60 0.361 1.008 ; qtot 0.361
1 ful H61 61 0.09 1.008 ; qtot 0.09
1 ful H62 62 0.09 1.008 ; qtot 0.09
1 ful H63 63 0.09 1.008 ; qtot 0.09
1 ful H64 64 0.414 1.008 ; qtot 0.414
1 ful H65 65 0.414 1.008 ; qtot 0.414
1 ful H66 66 0.39 1.008 ; qtot 0.39
1 ful H67 67 0.39 1.008 ; qtot 0.39
1 ful H68 68 0.09 1.008 ; qtot 0.09
1 ful H69 69 0.09 1.008 ; qtot 0.09
1 ful H70 70 0.414 1.008 ; qtot 0.414
1 ful H71 71 0.414 1.008 ; qtot 0.414
1 ful H72 72 0.09 1.008 ; qtot 0.09
1 ful H73 73 0.09 1.008 ; qtot 0.09
1 ful H74 74 0.09 1.008 ; qtot 0.09
1 ful H75 75 0.09 1.008 ; qtot 0.09
1 ful H76 76 0.09 1.008 ; qtot 0.09
1 ful H77 77 0.403 1.008 ; qtot 0.403
1 ful H78 78 0.349 1.008 ; qtot 0.349
1 ful H79 79 0.09 1.008 ; qtot 0.09
1 ful H80 80 0.09 1.008 ; qtot 0.09
1 ful H81 81 0.403 1.008 ; qtot 0.403
1 ful H82 82 0.09 1.008 ; qtot 0.09
1 ful H83 83 0.09 1.008 ; qtot 0.09
1 ful H84 84 0.414 1.008 ; qtot 0.414
1 ful H85 85 0.361 1.008 ; qtot 0.361
1 ful H86 86 0.09 1.008 ; qtot 0.09
1 ful H87 87 0.09 1.008 ; qtot 0.09
1 ful H88 88 0.09 1.008 ; qtot 0.09
CAD 1 ful CAD 89 2.00 112.4100

Now theoretically CAD has to interact with the two sulfur atoms located on each glutathione molecule, S15 and S35 in this case. There is no exact tutorial or protocol mentioned that I could find to make this work. Instead there are bits and pieces of information on the forum or old Gromacs discussion pages that I have used. Firstly I tried to update my specbond.dat file as follows:
8
CYS SG 1 CYS SG 1 0.2 CYS2 CYS2
CYS SG 1 HEM FE 2 0.25 CYS2 HEME
CYS SG 1 HEM CAB 1 0.18 CYS2 HEME
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
fu**l S15 1 ful CAD 2 0.245 S15 CAD **
ful S35 1 ful CAD 2 0.250 S35 CAD

The last two lines were the only lines added by me other stuff was already present. I measured the distances between S15 and CAD and S35 and CAD using Chimera to get those values.

Next, I updated the ful_ffbonded.itp file located in the Charmm36 directory provided by CGenFF. I did the following modifications:

[ bondtypes ]
; i j func b0 kb
CG311 NG311 1 0.14740001 220078.41
CG311 CG323 1 0.15380000 186188.00
S15 CAD 1 0.2450000 885440.00
S35 CAD 1 0.2500000 961540.00

The values of b0 and kb were found in the literature. Next, I also added the angles in the same file as follows:

[ angletypes ]
; i j k func theta0 ktheta ub0 kub
S15 CAD S35 1 109.500000 384.960000 0.00000000 0.00

These values are also from the literature. Next, in the ful.rtp file I added the following lines to the already existing section:
[ ful ]
[ atoms ]
CAD CAD 2.000 89
and
[ bonds ]
S15 CAD
S35 CAD
Next, in the residuetype.dat file I added the following line:
CAD ful
Next when I tried to give Pdb2gmx command it says the following:
Command line:
gmx pdb2gmx -f full_gmx.pdb -o ligand.pdb

Select the Force Field:

From current directory:

1: CHARMM36 all-atom force field (November 2018)

From ‘/opt/ohpc/pub/apps/gromacs/2022.5/share/gromacs/top’:

2: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003)

3: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)

4: AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996)

5: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000)

6: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65, 712-725, 2006)

7: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010)

8: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)

9: CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins)

10: GROMOS96 43a1 force field

11: GROMOS96 43a2 force field (improved alkane dihedrals)
12: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)

13: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)

14: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)

15: GROMOS96 54a7 force field (Eur. Biophys. J. (2011), 40, 843-856, DOI: 10.1007/s00249-011-0700-9)

16: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
1

Using the Charmm36 force field in directory ./charmm36.ff
Opening force field file ./charmm36.ff/watermodels.dat

Select the Water Model:

1: TIP3P TIP 3-point, recommended, by default uses CHARMM TIP3 with LJ on H

2: TIP4P TIP 4-point

3: TIP5P TIP 5-point

4: SPC simple point charge

5: SPC/E extended simple point charge
6: None
1

going to rename ./charmm36.ff/merged.r2b
Opening force field file ./charmm36.ff/merged.r2b
Reading full_gmx.pdb…
Read ‘’, 89 atoms

Analyzing pdb file
Splitting chemical chains based on TER records or chain id changing.

There are 1 chains and 0 blocks of water and 1 residues with 89 atoms

chain #res #atoms

1 ’ ’ 1 89

All occupancies are one
All occupancies are one
Opening force field file ./charmm36.ff/atomtypes.atp

Reading residue database… (Charmm36)
Opening force field file ./charmm36.ff/ful.rtp
Opening force field file ./charmm36.ff/merged.rtp
Opening force field file ./charmm36.ff/merged.hdb
Opening force field file ./charmm36.ff/merged.n.tdb
Opening force field file ./charmm36.ff/merged.c.tdb
Processing chain 1 (89 atoms, 1 residues)

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
Opening force field file ./charmm36.ff/merged.arn

Checking for duplicate atoms…

Generating any missing hydrogen atoms and/or adding termini.

Now there are 1 residues with 89 atoms

Making bonds…

Long Bond (35-89 = 0.250851 nm)

Number of bonds was 88, now 88

Generating angles, dihedrals and pairs…
Before cleaning: 230 pairs
Before cleaning: 230 dihedrals

Making cmap torsions…

There are 230 dihedrals, 0 impropers, 153 angles
230 pairs, 88 bonds and 0 virtual sites
Total mass 741.170 a.m.u.

Total charge 0.000 e

Writing topology

Writing coordinate file…

            --------- PLEASE NOTE ------------

You have successfully generated a topology from: full_gmx.pdb.

The Charmm36 force field and the tip3p water model are used.

            --------- ETON ESAELP ------------

GROMACS reminds you: “Don’t Follow Me Home” (Throwing Muses)

The issue is when I try to open the output file Ligand.pdb there are no bonds formed between S15-CAD and S35-CAD. I can see some lines about long bonds 35-89 but there are no bonds also why no bond between 15-89.

I have been trying to get help by posting similar questions in the last couple of months but I could not get any response. I am really looking forward for suggestions to navigate this problem. Kindly help me with this please.

Best regards

Ligand1 · January 28, 2025, 3:17am

Hi,

Any help in this regard will be highly appreciated.

Thanks in advance

Ligand1 · January 29, 2025, 2:29am

I am sorry but I think my post did not get any attention from experts. I would really need help from experts @jalemkul @MagnusL Please help to resolve this issue.

Thanks

MagnusL · January 29, 2025, 6:09am

If I had known I would have helped you already. Personally I’d avoid using pdb2gmx for anything that are not standard amino acid/RNA/DNA residues, even if you can add/modify rtp files. I would suggest modifying the itp files you get from CGenFF instead.

Your ligand is very unusual. How will you verify that your simulation results are reliable? I don’t think you can just add Cd2+ to glutathione and hope it will work well (even if you might get it to run). I think you will have to do the whole parameterization manually from the start. That would require an expert and quite a lot of effort.

Ligand1 · January 29, 2025, 8:06pm

Hi Magnus,

Many thanks for your comments and your suggestion. I agree that it would be a good idea to have an expert work on this problem since I do not have QM/MM experience. However, I do have a follow up question regarding my strategy. I agree that starting with Pdb2gmx might not be the best approach. Modifying the .itp file that I got from the CGenFF server looks a good option. Can you have a look at the modifications that I did for the other files such as specbond.dat, ful_ffbonded.itp, ful.rtp, and residuetype.dat and recommend if that looks OK to you. Is there anything else that needs modifications. With you inputs I can use the same strategy in the future projects if I have to include a new bond or a residue in the system.

Many thanks for helping me

Best regards

MagnusL · January 30, 2025, 8:08am

I’m afraid I won’t have time (tomorrow will be my last day at my current position) to look closely at your modifications (and I can’t find the files, perhaps the forum has not accepted them).

But in general, if you use (include it in the topology file) the generated itp file for a separate ligand, you usually don’t have to make any modifications to rtp and dat files. If you need to modify/add a residue in a protein/DNA/RNA etc then you’ll need to modify files used by pdb2gmx.

Topic		Replies	Views
Generating coordinates for ligand containing Cadmium using Charmm forcefield in Gromacs User discussions	1	30	January 17, 2025
Building an unusual ligand from Scratch User discussions	0	19	January 21, 2025
Add covalent ligand topology in CHARMM36 User discussions pdb2gmx	1	567	April 28, 2022
Alternative to CGenFF for generating large ligand topology User discussions forcefield	4	1051	April 6, 2024
How to create topology for a ligand User discussions	1	1550	June 15, 2020

Pdb2gmx not able to form bonds in an unusual ligand

Related topics