[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:
"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.
>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:
"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
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