[gmx-users] Using specbond.dat in pdb2gmx

Matthew Stancea mstancea at ggc.edu
Fri Jun 13 18:57:49 CEST 2014


Dr. Justin Lemkul,

Now I understand the message about removal of bonds being just duplicate bonds that have been removed, and, also, I understand how the -ter flag works now.

>> 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

>Here is precisely why your specbond.dat is not working - the bond criteria are
>not satisfied.  The Cys1(N)-Val29(C) distance is 0.132 nm.  You're specifing in
>specbond.dat (likely from simply copying the contents of the other lines) that
>the reference distance is 0.25 nm.  If Gromacs does not find atoms within ±10%
>of the value in specbond.dat, a bond won't be created.  Since 0.132 nm is way
>off, you won't get a bond.

>Also note that your final columns will try to rename the residues to CYS2 (which
>is specific for a disulfide cysteine, so unless that's true you'll get more
>fatal errors) and VAL2, which doesn't exist.

>See http://www.gromacs.org/Documentation/File_Formats/specbond.dat


Ah I see. I also noticed that in the original specbond.dat, the bond distance specification for the disulfide bonds between 2 CYS's or between 2 CYM's are .2 nm, not .25 nm, so I changed that as well. The only thing I still find a bit unclear is the final columns (residue rename). Is this what my specbond.dat file should look like?

3
CYS     N       1       VAL     C       1       0.132    CYS    VAL
CYS     SG      1       CYS     SG      1       0.20    CYS2    CYS2
CYM     SG      1       CYM     SG      1       0.20    CYS2    CYS2

-Matthew Stancea

(This email chain has gotten a bit long in my opinion, so I have put the previous email below. With your permission, I would prefer just using the above message in order to continue correspondence in shorter emails; however, if you would prefer to instead use the entire email, I completely understand.)
________________________________________
From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Justin Lemkul <jalemkul at vt.edu>
Sent: Thursday, June 12, 2014 5:01 PM
To: gmx-users at gromacs.org
Subject: Re: [gmx-users] Using specbond.dat in pdb2gmx

On 6/12/14, 4:53 PM, Matthew Stancea wrote:
>> 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:
>
> "3
> 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.
>

This is related to removing duplicate bonds after merging .rtp and .tdb entries.
  It's normal output; there are almost always duplicates.

>> 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
> -------------------------------------------------------
> "
>

The choices are 0 or 1, so "none" is not a valid response.

> 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).
>

This is what -ter should be doing - asking you about every single possible
protonation state.

>>> 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:
>
> "3
> 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

Here is precisely why your specbond.dat is not working - the bond criteria are
not satisfied.  The Cys1(N)-Val29(C) distance is 0.132 nm.  You're specifing in
specbond.dat (likely from simply copying the contents of the other lines) that
the reference distance is 0.25 nm.  If Gromacs does not find atoms within ±10%
of the value in specbond.dat, a bond won't be created.  Since 0.132 nm is way
off, you won't get a bond.

Also note that your final columns will try to rename the residues to CYS2 (which
is specific for a disulfide cysteine, so unless that's true you'll get more
fatal errors) and VAL2, which doesn't exist.

See http://www.gromacs.org/Documentation/File_Formats/specbond.dat

>                     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.
>

Glad to hear it!

-Justin

--
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================
--
Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.


More information about the gromacs.org_gmx-users mailing list