[gmx-users] Using specbond.dat in pdb2gmx

Matthew Stancea mstancea at ggc.edu
Thu Jun 12 22:53:52 CEST 2014

>On 6/11/14, 2:33 PM, Matthew Stancea wrote:
>> ? Hello,
>> I have been having a bit of issues generating an accurate gromacs topology
>> file ("topol.top") utilizing a pdb of a cyclic peptide in pdb2gmx. I have
>> been able to generate topologies that are almost identical to the original
>> pdb using the option "-ignh" at the end of my command, but doing so deletes
>> the bond between the first nitrogen and the last carbon (which should be
>> connected for this peptide to be cyclic) and adds two additional hydrogens
>> and a positive charge to the first nitrogen and an additional oxygen and a
>> negative charge to the last carbon, rendering this peptide as non-cyclic.
>> After searching around for quite a while, I found out that many others on
>> this mailing list were having the same issues as myself, and some replies to
>> their messages including the usage of a file known as "specbond.dat" which
>> may be helpful in retaining that bond.

>The -ignh flag is not relevant to those observations.  

In my working directory, I have a file named specbond.dat and it contains the following information:

CYS     N       1       VAL     C       1       0.25    CYS2    VAL2
CYS     SG      1       CYS     SG      1       0.25    CYS2    CYS2
CYM     SG      1       CYM     SG      1       0.25    CYS2    CYS2

When I input the command "pdb2gmx_mpi -ff amber99sb -water tip3p -f 1NB1.pdb -o conf.gro -p topol.top -i kalata.itp" while specbond.dat is in the working directory (Here is a link for the structure of peptide 1NB1 for reference if needed: http://www.ncbi.nlm.nih.gov/protein/1NB1_A ), I get the following message: 

Program pdb2gmx_mpi, VERSION 4.6.2
Source code file: /home/msaum/apps/gromacs-4.6.2/src/kernel/pdb2gmx.c, line: 727

Fatal error:
Atom HB3 in residue CYS 1 was not found in rtp entry NCYS with 13 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 http://www.gromacs.org/Documentation/Errors

(HB3 would be the hydrogen of a cysteine residue that has not formed a disulfide bond with another amino acid residue, however my pdb does not contain an HB3).

I then tried adding -ignh as per the suggestion in the twelfth line in that message. The following is a part of the message I received:

"Checking for duplicate atoms....
Generating any missing hydrogen atoms and/or adding termini.
Now there are 29 residues with 385 atoms
Making bonds...
Number of bonds was 390, now 389
Generating angles, dihedrals and pairs...
Before cleaning: 1010 pairs
Before cleaning: 1030 dihedrals
Keeping all generated dihedrals
Making cmap torsions...There are 1030 dihedrals,   73 impropers,  699 angles
          1007 pairs,      389 bonds and     0 virtual sites
Total mass 2916.376 a.m.u.
Total charge -0.000 e
Writing topology

Writing coordinate file...
                --------- PLEASE NOTE ------------
You have successfully generated a topology from: 1NB1.pdb.
The Amber99sb force field and the tip3p water model are used.
                --------- ETON ESAELP ------------

On the fifth line of that excerpt, I read "Number of bonds was 390, now 389" and questioned what that could mean. Because of that, I wanted to verify that the topology that was "successfully generated" contains all the necessary bonds, so I looked through the topol.top file and saw 3 total hydrogens on the terminal nitrogen on the beginning cysteine and two total oxygens on the terminal carbon on the ending valine. Also, I did not see a bond between that nitrogen and that carbon (atoms 1 and 383); these data lead me to believe that this peptide is no longer cyclic.

>Construction of termini and assignment of ionization 
>state is done with -ter.  For your case, you likely
>need to be using -ter and selecting "None."

When I input the command "pdb2gmx_mpi -ff amber99sb -water tip3p -f 1NB1.pdb -o conf.gro -p topol.top -i kalata.itp -ter", I received the same message as before without any kind of interactivity for selecting options. Because of this, I tried the command "pdb2gmx_mpi -ff amber99sb -water tip3p -f 1NB1.pdb -o conf.gro -p topol.top -i kalata.itp -inter -ter", which asked about the protonation of residue 24 (arginine), and when I answer "none", I get this message:

"Processing chain 1 'A' (376 atoms, 29 residues)
Which ARGININE type do you want for residue 24
0. Not protonated (charge 0) (-)
1. Protonated (charge +1) (ARG)

Type a number:none

Program pdb2gmx_mpi, VERSION 4.6.2
Source code file: /home/msaum/apps/gromacs-4.6.2/src/kernel/pdb2gmx.c, line: 120

Fatal error:
Answer me for res ARGININE 24!
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors

If I  repeat the command and answer it, I am then asked about the protonation of residue 3 (glumatic acid). When I answer it, I then get the original error message again (the first message).

>> However, I could not find an explanation of how to use this file in my
>> command, and I was hoping someone could guide me on how to use specbond.dat
>> in order to specify a particular bond. (I have tried both "-chainsep
>> specbond.dat" and "-merge specbond.dat".)

>specbond.dat is a standard file that is always read by pdb2gmx; it is not an
>input option.  It lives in $GMXLIB but can be overridden by a copy in the local
>directory (as with all Gromacs library files).

I put a copy in the local directory.

>> I have been able to make a specbond.dat file that thoroughly explains what
>> special bonds I need pdb2gmx to account for, but when I have tried to type a
>> command file to use it, pdb2gmx generates topology files nearly identical to
>> the topologies of the non-cyclic peptides when I had used the "-ignh" option
>> rather than generating a cyclic peptide topol.top.
>> If more specific information is needed, I can provide that; however, I am not
>> looking for an answer that is specific to my own issue, but rather a tutorial
>> into generating a topology for gromacs from a pdb file of a cyclic peptide by
>> using specbond.dat with pdb2gmx. (My version of pdb2gmx is 4.6.2.)

>Post the contents of your version of specbond.dat, the exact command you're
>issuing, and any screen output related to parsing of specbond.dat.  

Here are the contents of my version of specbond.dat:

CYS     N       1       VAL     C       1       0.25    CYS2    VAL2
CYS     SG      1       CYS     SG      1       0.25    CYS2    CYS2
CYM     SG      1       CYM     SG      1       0.25    CYS2    CYS2

The exact command I am issuing is "pdb2gmx_mpi -ff amber99sb -water tip3p -f 1NB1.pdb -o conf.gro -p topol.top -i kalata.itp".

This is the screen output related to parsing of specbond.dat:

"3 out of 3 lines of specbond.dat converted successfully
Special Atom Distance matrix:
                    CYS1    CYS1    CYS5    CYS5    VAL6   CYS10   CYS10
                      N1     SG6     N47    SG52     C59    N101   SG106
    CYS1     SG6   0.302
    CYS5     N47   0.988   0.900
    CYS5    SG52   1.070   0.917   0.333
    VAL6     C59   1.498   1.396   0.510   0.585
   CYS10    N101   1.122   0.901   0.614   0.609   0.835
   CYS10   SG106   0.879   0.655   0.490   0.422   0.857   0.322
   CYS15    N160   0.784   0.488   1.066   0.984   1.477   0.774   0.620
   CYS15   SG165   0.477   0.202   0.960   0.961   1.432   0.823   0.625
   CYS17    N184   1.030   0.784   0.721   0.460   1.004   0.617   0.387
   CYS17   SG189   1.224   1.039   0.514   0.204   0.625   0.621   0.487
   VAL21    C245   0.756   0.685   0.408   0.388   0.847   0.812   0.531
   CYS22    N259   0.657   0.574   0.404   0.432   0.887   0.742   0.458
   CYS22   SG264   0.755   0.588   0.371   0.430   0.827   0.436   0.203
   VAL29    C363   0.132   0.394   1.117   1.193   1.626   1.248   1.004
                   CYS15   CYS15   CYS17   CYS17   VAL21   CYS22   CYS22
                    N160   SG165    N184   SG189    C245    N259   SG264
   CYS15   SG165   0.320
   CYS17    N184   0.661   0.767
   CYS17   SG189   1.017   1.056   0.412
   VAL21    C245   0.933   0.802   0.602   0.569
   CYS22    N259   0.836   0.686   0.594   0.612   0.133
   CYS22   SG264   0.701   0.606   0.533   0.568   0.419   0.321
   VAL29    C363   0.854   0.555   1.137   1.345   0.869   0.777   0.886

>You should only need to add one bond to the one that is provided by default.  Your version
>of specbond.dat is either in $GMXLIB or the working directory, correct?

After adding my 1 line to a copy of the default specbond.dat file, I removed any line that I knew was irrelevant for this peptide.

And yes, my specbond.dat is in the working directory now, but as mentioned above, I am still having the same issues.

-Matthew Stancea

Also, I would like to thank you for you help so far, Dr. Lemkul. I did not previously mention this, but I am an undergraduate working on a molecular dynamics research project with my adviser, and he has actually told me that he has been benefited from your solutions to issues that others have had and your gromacs support over the years.

More information about the gromacs.org_gmx-users mailing list