[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