Dear All,
I just want to share some recent experiences with sleap and hopefully
suggest some better feedback to users and highlight some behaviour
differences between it and {x,t}leap.
Recently I have been following a pretty standard protocol to prepare a
prmtop/inpcrd pair. Initially, I had used an external tool to determine
the protonation state of various residues in a PDB file
(
http://www.poissonboltzmann.org/pdb2pqr) and then, using these results,
preprocessed the raw PDB file via a simple sed script translating HIS to
HIP, HID, or HIE. Some other residues (crystal waters) were also being
stripped in this stage.
Issue #1
========
This PDB triggered an assert fail when loaded via the loadpdb command:
sleap: second.cpp:73: void mort::pdbent::apply_second(const
mort::pdbent::second_s&, mort::molecule_t&): Assertion `start !=
mol.resd_end() && end != mol.resd_end()' failed.
and after quite a lot of debugging, I found it out what the cause was.
PDB files have records pertaining to secondary structure[1]:
HELIX 1 1 GLY A 31 LEU A 36 1
This states that there a helix secondary structure 1 in the residue
range 31 to 36. Sleap parses this information and the assert was being
triggered by the fact that within the HELIX records, one of the HIS
residues lay on a boundary, but did not match its renamed form (HID) in
the corresponding ATOM records. This is a valid error and is essentially
my fault for not being rigorous in my renaming script. However,
{t,x}leap accepted this fine and I can only assume that they do not
process and HELIX/SHEET records.
This extra level of checking is nice, but it would be useful to output
an informative error message to the user stating that there is some form
of residue naming mismatch between the ATOM and HELIX/SHEET records.
Issue #2
========
Having correct for above, and attempting to save the prmtop/inpcrd pair,
an error was generated about a missing bond, which was being caused by a
bond being created between the first atom in the system and an atom in
a HEATM record. It was unclear why this bond was being generated.
Again, after lots of debugging, it turned out that the stripping process
(in the preamble) had removed one ATOM record that a CONECT record [2]
pair was referring to. As a result mort::pdbent::read_bond() was bonding
the present ATOM to a non-existent atom and it defaulted to ATOM serial
number "1" if the stated CONECT atom did not exist, resulting in a
spurious bond.
Again, this was my fault for the hacky modification of the PDB, however,
mort::pdbent::read_bond() should throw a warning if it cannot find the
respective atoms in a CONECT record and *not* proceed to create a bond
with atom 1.
I am not (yet) familiar with the sleap code yet to put these fixed in
myself, but nor do I want to skip over this, since I think this may bite
someone in the future.
Regards,
Mark
Refs
====
[1]
http://www.wwpdb.org/documentation/format33/sect5.html
[2]
http://www.wwpdb.org/documentation/format33/sect10.html
_______________________________________________
AMBER-Developers mailing list
AMBER-Developers.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber-developers
Received on Wed Nov 02 2011 - 10:00:05 PDT
