From 754d26ff0fd1845a096d59afa0cbb4476a148f67 Mon Sep 17 00:00:00 2001
From: Marco Biasini <marco.biasini@unibas.ch>
Date: Thu, 3 May 2012 22:07:33 +0200
Subject: [PATCH] fix BZDNG318
---
modules/io/src/mol/pdb_reader.cc | 14 +++---
modules/io/src/mol/pdb_writer.cc | 4 +-
modules/io/tests/test_io_pdb.cc | 19 +++++++
modules/io/tests/testfiles/pdb/bzdng-318.pdb | 4 ++
modules/mol/base/src/atom_base.cc | 11 ++++
modules/mol/base/src/atom_base.hh | 2 +
modules/mol/base/src/editor_base.cc | 13 +++--
modules/mol/base/src/editor_base.hh | 6 ++-
modules/mol/base/src/impl/atom_group.hh | 7 ++-
modules/mol/base/src/impl/residue_impl.cc | 53 +++++++++++++++++---
modules/mol/base/src/impl/residue_impl.hh | 10 ++--
11 files changed, 116 insertions(+), 27 deletions(-)
create mode 100644 modules/io/tests/testfiles/pdb/bzdng-318.pdb
diff --git a/modules/io/src/mol/pdb_reader.cc b/modules/io/src/mol/pdb_reader.cc
index c52e5cf31..6f5ac4633 100644
--- a/modules/io/src/mol/pdb_reader.cc
+++ b/modules/io/src/mol/pdb_reader.cc
@@ -808,12 +808,14 @@ void PDBReader::ParseAndAddAtom(const StringRef& line, int line_num,
return;
}
}
+ Real b=temp.first ? temp.second : 0.0;
+ Real o=occ.first ? occ.second : 1.0;
if (!profile_.quack_mode && alt_loc!=' ') {
// Check if there is already a atom with the same name.
mol::AtomHandle me=curr_residue_.FindAtom(aname);
if (me.IsValid()) {
try {
- editor.AddAltAtomPos(String(1, alt_loc), me, apos);
+ editor.AddAltAtomPos(String(1, alt_loc), me, apos, o, b);
} catch (Error) {
LOG_INFO("Ignoring atom alt location since there is already an atom "
"with name " << aname << ", but without an alt loc");
@@ -822,7 +824,7 @@ void PDBReader::ParseAndAddAtom(const StringRef& line, int line_num,
return;
} else {
ah=editor.InsertAltAtom(curr_residue_, aname,
- String(1, alt_loc), apos, s_ele);
+ String(1, alt_loc), apos, s_ele, o, b);
++atom_count_;
}
} else {
@@ -844,12 +846,8 @@ void PDBReader::ParseAndAddAtom(const StringRef& line, int line_num,
ah.SetRadius(radius.second);
}
}
- if (temp.first) {
- ah.SetBFactor(temp.second);
- }
- if (occ.first) {
- ah.SetOccupancy(occ.second);
- }
+ ah.SetBFactor(b);
+ ah.SetOccupancy(o);
if (charge.first) {
ah.SetCharge(charge.second);
}
diff --git a/modules/io/src/mol/pdb_writer.cc b/modules/io/src/mol/pdb_writer.cc
index f3c863594..fdae17a51 100644
--- a/modules/io/src/mol/pdb_writer.cc
+++ b/modules/io/src/mol/pdb_writer.cc
@@ -172,8 +172,8 @@ void write_atom(std::ostream& ostr, FormattedLine& line,
line(54, 6)=fmt::LPaddedFloat(atom.GetCharge(), 2);
line(60, 6)=fmt::LPaddedFloat(atom.GetRadius(), 2);
} else {
- line(54, 6)=fmt::LPaddedFloat(atom.GetOccupancy(), 2);
- Real bfac=atom.GetBFactor();
+ line(54, 6)=fmt::LPaddedFloat(atom.GetAltOcc(*i), 2);
+ Real bfac=atom.GetAltBFactor(*i);
if (bfac>999.99) {
line(60, 6)=fmt::LPaddedFloat(999.99, 2);
} else {
diff --git a/modules/io/tests/test_io_pdb.cc b/modules/io/tests/test_io_pdb.cc
index 1e7ecde2f..cb40a2c04 100644
--- a/modules/io/tests/test_io_pdb.cc
+++ b/modules/io/tests/test_io_pdb.cc
@@ -424,6 +424,25 @@ BOOST_AUTO_TEST_CASE(deuterium_import)
BOOST_CHECK(ent.FindResidue("A", 297).IsPeptideLinking());
}
+BOOST_AUTO_TEST_CASE(bzdng_318)
+{
+ String fname("testfiles/pdb/bzdng-318.pdb");
+ PDBReader reader(fname, IOProfile());
+ mol::EntityHandle ent=mol::CreateEntity();
+ reader.Import(ent);
+ // we use conopology to mark amino acids as peptide-linking.
+ conop::Conopology& conop_inst=conop::Conopology::Instance();
+ conop_inst.ConnectAll(conop_inst.GetBuilder(), ent);
+ {
+ PDBWriter writer(std::string("testfiles/pdb/bzdng-318-out.pdb"),
+ IOProfile());
+ writer.Write(ent);
+ }
+
+ BOOST_CHECK(compare_files("testfiles/pdb/bzdng-318.pdb",
+ "testfiles/pdb/bzdng-318-out.pdb"));
+}
+
BOOST_AUTO_TEST_CASE(faulty_lines)
{
String fname("testfiles/pdb/faulty.pdb");
diff --git a/modules/io/tests/testfiles/pdb/bzdng-318.pdb b/modules/io/tests/testfiles/pdb/bzdng-318.pdb
new file mode 100644
index 000000000..8008e84b1
--- /dev/null
+++ b/modules/io/tests/testfiles/pdb/bzdng-318.pdb
@@ -0,0 +1,4 @@
+ATOM 1 CA ALYS A 9 31.209 -9.896 12.154 0.62 22.96 C
+ATOM 1 CA BLYS A 9 31.213 -9.902 12.145 0.38 24.64 C
+TER 2 LYS A 9
+END
\ No newline at end of file
diff --git a/modules/mol/base/src/atom_base.cc b/modules/mol/base/src/atom_base.cc
index 642ea5aa9..c818fbc1c 100644
--- a/modules/mol/base/src/atom_base.cc
+++ b/modules/mol/base/src/atom_base.cc
@@ -69,6 +69,17 @@ geom::Vec3 AtomBase::GetAltPos(const String& alt_group) const
return impl_->GetResidue()->GetAltAtomPos(Impl(), alt_group);
}
+Real AtomBase::GetAltBFactor(const String& alt_group) const
+{
+ this->CheckValidity();
+ return impl_->GetResidue()->GetAltAtomBFactor(Impl(), alt_group);
+}
+Real AtomBase::GetAltOcc(const String& alt_group) const
+{
+ this->CheckValidity();
+ return impl_->GetResidue()->GetAltAtomOcc(Impl(), alt_group);
+}
+
std::vector<String> AtomBase::GetAltGroupNames() const
{
this->CheckValidity();
diff --git a/modules/mol/base/src/atom_base.hh b/modules/mol/base/src/atom_base.hh
index 8d19514af..137953115 100644
--- a/modules/mol/base/src/atom_base.hh
+++ b/modules/mol/base/src/atom_base.hh
@@ -87,6 +87,8 @@ public:
const geom::Vec3& GetOriginalPos() const;
/// \brief get alternative atom position
geom::Vec3 GetAltPos(const String& alt_group) const;
+ Real GetAltBFactor(const String& alt_group) const;
+ Real GetAltOcc(const String& alt_group) const;
std::vector<String> GetAltGroupNames() const;
diff --git a/modules/mol/base/src/editor_base.cc b/modules/mol/base/src/editor_base.cc
index ba822f0f5..d75827806 100644
--- a/modules/mol/base/src/editor_base.cc
+++ b/modules/mol/base/src/editor_base.cc
@@ -119,21 +119,26 @@ AtomHandle EditorBase::InsertAtom(ResidueHandle res, const String& name,
AtomHandle EditorBase::InsertAltAtom(ResidueHandle res, const String& name,
const String& alt_group,
const geom::Vec3& pos,
- const String& ele)
+ const String& ele, Real occ,
+ Real b_factor)
{
CheckHandleValidity(res);
ent_.Impl()->MarkTraceDirty();
- AtomHandle atom(res.Impl()->InsertAltAtom(name, alt_group, pos, ele));
+ AtomHandle atom(res.Impl()->InsertAltAtom(name, alt_group, pos,
+ ele, occ, b_factor));
this->UpdateTrace();
return atom;
}
void EditorBase::AddAltAtomPos(const String& group,
const AtomHandle& atom,
- const geom::Vec3& position)
+ const geom::Vec3& position, Real occ,
+ Real b_factor)
{
CheckHandleValidity(atom);
- atom.GetResidue().Impl()->AddAltAtomPos(group, atom.Impl(), position);
+ atom.GetResidue().Impl()->AddAltAtomPos(group, atom.Impl(),
+ position, occ, b_factor);
+
}
void EditorBase::DeleteChain(const ChainHandle& chain)
diff --git a/modules/mol/base/src/editor_base.hh b/modules/mol/base/src/editor_base.hh
index c6309cce3..504311d54 100644
--- a/modules/mol/base/src/editor_base.hh
+++ b/modules/mol/base/src/editor_base.hh
@@ -124,7 +124,8 @@ public:
/// \sa EditorBase::AddAltAtomPos(), ResidueHandle
AtomHandle InsertAltAtom(ResidueHandle residue, const String& name,
const String& alt_group, const geom::Vec3& pos,
- const String& ele="");
+ const String& ele="", Real occ=1.0,
+ Real b_factor=0.0);
/// \brief Add alternative atom position
/// \param group is the name of the alternative atom position group. If no
/// group of that name exists, it will be created.
@@ -136,7 +137,8 @@ public:
/// is the alternative position
/// \sa EditorBase::InsertAltAtom(), ResidueHandle
void AddAltAtomPos(const String& group, const AtomHandle& atom,
- const geom::Vec3& position);
+ const geom::Vec3& position, Real occ=1.0,
+ Real b_factor=0.0);
//\}
/// \brief connect two atoms with bond
diff --git a/modules/mol/base/src/impl/atom_group.hh b/modules/mol/base/src/impl/atom_group.hh
index 934fda241..50a7c6338 100644
--- a/modules/mol/base/src/impl/atom_group.hh
+++ b/modules/mol/base/src/impl/atom_group.hh
@@ -29,10 +29,13 @@ namespace ost { namespace mol {namespace impl {
struct AtomGroupEntry {
AtomGroupEntry() {}
- AtomGroupEntry(const AtomImplPtr& a, const geom::Vec3& p)
- : atom(a), pos(p) {}
+ AtomGroupEntry(const AtomImplPtr& a, const geom::Vec3& p,
+ Real o, Real b)
+ : atom(a), pos(p), occ(o), b_factor(b) {}
impl::AtomImplW atom;
geom::Vec3 pos;
+ Real occ;
+ Real b_factor;
};
typedef std::vector<AtomGroupEntry> AtomGroupEntryList;
diff --git a/modules/mol/base/src/impl/residue_impl.cc b/modules/mol/base/src/impl/residue_impl.cc
index eb633c5be..1e0fa1d46 100644
--- a/modules/mol/base/src/impl/residue_impl.cc
+++ b/modules/mol/base/src/impl/residue_impl.cc
@@ -91,7 +91,8 @@ Real ResidueImpl::GetAverageBFactor() const
}
void ResidueImpl::AddAltAtom(const String& group, const AtomImplPtr& atom,
- const geom::Vec3& position) {
+ const geom::Vec3& position,
+ Real occ, Real b_factor) {
if (group.length()==0)
throw Error("alt atom group name can't be empty String");
@@ -106,7 +107,7 @@ void ResidueImpl::AddAltAtom(const String& group, const AtomImplPtr& atom,
curr_group_=group;
}
}
- i->second.atoms.push_back(AtomGroupEntry(atom, position));
+ i->second.atoms.push_back(AtomGroupEntry(atom, position, occ, b_factor));
}
geom::Vec3 ResidueImpl::GetAltAtomPos(const AtomImplPtr& atom,
@@ -127,12 +128,47 @@ geom::Vec3 ResidueImpl::GetAltAtomPos(const AtomImplPtr& atom,
" does not have alternative atom position '"+group+"'");
}
+Real ResidueImpl::GetAltAtomOcc(const AtomImplPtr& atom,
+ const String& group) const
+{
+ AtomEntryGroups::const_iterator i=alt_groups_.find(group);
+ if (i==alt_groups_.end()) {
+ throw Error("No alt atom group '"+group+"'");
+ }
+ const AtomGroup& g=i->second;
+ for (AtomGroupEntryList::const_iterator j=g.atoms.begin(),
+ e=g.atoms.end(); j!=e; ++j) {
+ if (atom==j->atom.lock()) {
+ return j->occ;
+ }
+ }
+ throw Error(atom->GetQualifiedName()+
+ " does not have alternative atom position '"+group+"'");
+}
+Real ResidueImpl::GetAltAtomBFactor(const AtomImplPtr& atom,
+ const String& group) const
+{
+ AtomEntryGroups::const_iterator i=alt_groups_.find(group);
+ if (i==alt_groups_.end()) {
+ throw Error("No alt atom group '"+group+"'");
+ }
+ const AtomGroup& g=i->second;
+ for (AtomGroupEntryList::const_iterator j=g.atoms.begin(),
+ e=g.atoms.end(); j!=e; ++j) {
+ if (atom==j->atom.lock()) {
+ return j->b_factor;
+ }
+ }
+ throw Error(atom->GetQualifiedName()+
+ " does not have alternative atom position '"+group+"'");
+}
AtomImplPtr ResidueImpl::InsertAltAtom(const String& name,
const String& alt_group,
const geom::Vec3& pos,
- const String& ele) {
+ const String& ele,
+ Real occupancy, Real b_factor) {
AtomImplPtr atom=this->InsertAtom(name, pos, ele);
- this->AddAltAtom(alt_group, atom, pos);
+ this->AddAltAtom(alt_group, atom, pos, occupancy, b_factor);
return atom;
}
@@ -506,7 +542,8 @@ geom::Vec3 ResidueImpl::GetCenterOfMass() const
void ResidueImpl::AddAltAtomPos(const String& group,
const AtomImplPtr& atom,
- const geom::Vec3& position) {
+ const geom::Vec3& position,
+ Real occ, Real b_factor) {
// Make sure atom is already registered for having an alternative position.
// Bail out, if this is not the case.
AtomEntryGroups::iterator i=alt_groups_.begin();
@@ -522,7 +559,7 @@ void ResidueImpl::AddAltAtomPos(const String& group,
}
}
if (found)
- this->AddAltAtom(group, atom, position);
+ this->AddAltAtom(group, atom, position, occ, b_factor);
else {
String m="Definition of alternative position without prior call to "
"InsertAltAtom is not allowed";
@@ -582,6 +619,8 @@ bool ResidueImpl::SwitchAtomPos(const String& group) {
AtomGroupEntry& entry=*k;
assert(!entry.atom.expired());
entry.pos=entry.atom.lock()->OriginalPos();
+ entry.occ=entry.atom.lock()->GetOccupancy();
+ entry.b_factor=entry.atom.lock()->GetBFactor();
}
}
AtomGroup& agr=i->second;
@@ -595,6 +634,8 @@ bool ResidueImpl::SwitchAtomPos(const String& group) {
geom::Mat4 transf_matrix = ent.GetTransformationMatrix();
geom::Vec3 transf_pos = geom::Vec3(transf_matrix*geom::Vec4(entry.pos));
entry.atom.lock()->TransformedPos()=transf_pos;
+ entry.atom.lock()->SetBFactor(j->b_factor);
+ entry.atom.lock()->SetOccupancy(j->occ);
}
curr_group_=group;
return true;
diff --git a/modules/mol/base/src/impl/residue_impl.hh b/modules/mol/base/src/impl/residue_impl.hh
index b52d26fc7..24f829c69 100644
--- a/modules/mol/base/src/impl/residue_impl.hh
+++ b/modules/mol/base/src/impl/residue_impl.hh
@@ -58,7 +58,8 @@ public:
/// no bonds
AtomImplPtr InsertAtom(const AtomImplPtr& atom);
AtomImplPtr InsertAltAtom(const String& name, const String& alt_group,
- const geom::Vec3& pos, const String& ele);
+ const geom::Vec3& pos, const String& ele,
+ Real occ, Real b_factor);
const ResNum& GetNumber() const {return num_;}
void SetNumber(const ResNum& num) {num_=num;}
@@ -167,8 +168,11 @@ public:
void AddAltAtomPos(const String& group,
const AtomImplPtr& atom,
- const geom::Vec3& position);
+ const geom::Vec3& position,
+ Real occ, Real b_factor);
geom::Vec3 GetAltAtomPos(const AtomImplPtr& atom, const String& group) const;
+ Real GetAltAtomOcc(const AtomImplPtr& atom, const String& group) const;
+ Real GetAltAtomBFactor(const AtomImplPtr& atom, const String& group) const;
const String& GetCurrentAltGroupName() const {
@@ -221,7 +225,7 @@ public:
void SetIsLigand(bool flag) { ligand_=flag; }
private:
void AddAltAtom(const String& group, const AtomImplPtr& atom,
- const geom::Vec3& position);
+ const geom::Vec3& position, Real occ, Real b_factor);
void RemoveAltPositionsForAtom(const AtomImplPtr& atom);
String curr_group_;
--
GitLab