Pdb2gmx does not make disulfide bond between peptide chains

Gromacs version: 2021.5
Modification: No
Force Field: Charmm C36M


I’m trying to generate the topology for a peptide with 2 identical chains linked by a disulfide bond at the n-terminus. But I can’t get pdb2gmx to create the disulfide bond. pdb2gmx creates the topology without the disulfide bond. I even tried to modify the specbond.dat file to allow disulfide bonds for longer distances but it doesn’t work.

I’ve tried the following commands:

gmx pdb2gmx -f dimer.pdb -o dimer_ss.gro -ignh -ss

gmx pdb2gmx -f dimer.pdb -o dimer_ss.gro -ignh -ss -merge all

gmx pdb2gmx -f dimer.pdb -o dimer_ss.gro -ignh -ss -merge inter

gmx pdb2gmx -f dimer.pdb -o dimer_ss.gro -ignh -ss -merge all -chainsep ter

None of them worked.

My files can be accessed in the following link: Google Drive

Any help will be appreciated.

This is the correct approach, but the file you have provided in the Google Drive link called “dimer.pdb” has two identical chains (with identical coordinates) separated by an END card, so it seems to me only one chain would ever be read. But you can’t have two chains with the same coordinates anyway, so something is simply wrong with the coordinate file.

It is essential when diagnosing pdb2gmx issues to provide the complete, copied and pasted, terminal output. There is a ton of useful information there that can be used for troubleshooting.

Good evening, thank you for your reply professor.

I’ve updated the files in the link to include the fulll pdb2gmx output pdb2gmx_output.txt. Also I’ve tried creating a new structure to see if it would work, but it didn’t. It’s odd that pdb2gmx has written the topology for both chains but no S-S bond. Also it’s weird that visualizing with pymol, showed 2 chains with different postions.

Except for the forcefield options that I couldn’t fit in here the output was:

Back Off! I just backed up topol.top to ./#topol.top.1#

Processing chain 1 (92 atoms, 14 residues)
Analysing hydrogen-bonding network for automated assignment of histidine
 protonation. 18 donors and 14 acceptors were found.
There are 14 hydrogen bonds
Will use HISE for residue 6
Will use HISE for residue 6

Identified residue CYS2 as a starting terminus.

Identified residue LEU8 as a ending terminus.

Identified residue CYS2 as a starting terminus.

Identified residue LEU8 as a ending terminus.
8 out of 8 lines of specbond.dat converted successfully
Special Atom Distance matrix:
                    CYS2    HIS6    CYS2
                     SG6   NE230    SG52
    HIS6   NE230   0.529
    CYS2    SG52   0.093   0.470
    HIS6   NE276   0.487   0.943   0.529
Start terminus CYS-2: NH3+
End terminus LEU-8: COO-
Start terminus CYS-2: NH3+
End terminus LEU-8: COO-
Opening force field file ./charmm36-jul2022-napht.ff/aminoacids.arn

Checking for duplicate atoms....

Generating any missing hydrogen atoms and/or adding termini.

Now there are 14 residues with 192 atoms

Making bonds...

Number of bonds was 192, now 192

Generating angles, dihedrals and pairs...
Before cleaning: 478 pairs
Before cleaning: 488 dihedrals
Keeping all generated dihedrals

Making cmap torsions...

There are   10 cmap torsion pairs

There are  488 dihedrals,   38 impropers,  346 angles
           478 pairs,      192 bonds and     0 virtual sites

Total mass 1367.658 a.m.u.

Total charge 0.000 e

Writing topology

Back Off! I just backed up posre.itp to ./#posre.itp.1#

Writing coordinate file...

Back Off! I just backed up dimer_ss.gro to ./#dimer_ss.gro.1#

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

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

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

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

GROMACS reminds you: "The Feeling of Power was Intoxicating, Magic" (Frida Hyvonen)

Although it has sucessfully written the topology it does not contain the ss bond

The two SG atoms are too close to one another. Their distance in the coordinate file is 0.093 nm, according to the screen output. The typical S-S length is on the order of 0.2 nm, which is what specbond.dat uses. The tolerance for this value is ±10%, so the distance needs to be within 0.18 - 0.22 nm for it to write the bond.

Thank you very much for your replay professor, it finally worked.

Hello @jalemkul

Is there any entry of S-S bonds in topology file?
If it is, under which section/header can I confirm?

I have created topology with and without -ss flag. Both gave me same results.