Skip to content

L
libcifpp
Commit f7b98c05 by Maarten L. Hekkelman
Options

refactored AtomImpl

parent c4f3b1cd
Showing with 331 additions and 621 deletions
include/cif++/Structure.hpp
...@@ -60,30 +60,101 @@ class File; ...@@ -60,30 +60,101 @@ class File;
class Atom class Atom
{ {
private:
struct AtomImpl : public std::enable_shared_from_this<AtomImpl>
{
AtomImpl(cif::Datablock &db, const std::string &id, cif::Row row);
// constructor for a symmetry copy of an atom
AtomImpl(const AtomImpl &impl, const Point &loc, const std::string &sym_op);
AtomImpl(const AtomImpl &i) = default;
void prefetch();
int compare(const AtomImpl &b) const;
bool getAnisoU(float anisou[6]) const;
void moveTo(const Point &p);
const Compound &comp() const;
const std::string get_property(const std::string_view name) const;
void set_property(const std::string_view name, const std::string &value);
const cif::Datablock &mDb;
std::string mID;
AtomType mType;
std::string mAtomID;
std::string mCompID;
std::string mAsymID;
int mSeqID;
std::string mAltID;
std::string mAuthSeqID;
Point mLocation;
int mRefcount;
cif::Row mRow;
mutable std::vector<std::tuple<std::string,cif::detail::ItemReference>> mCachedRefs;
mutable const Compound *mCompound = nullptr;
bool mSymmetryCopy = false;
bool mClone = false;
std::string mSymmetryOperator = "1_555";
};
public: public:
Atom();
Atom(struct AtomImpl *impl); Atom() {}
Atom(const Atom &rhs);
Atom(std::shared_ptr<AtomImpl> impl)
: mImpl(impl) {}
Atom(const Atom &rhs)
: mImpl(rhs.mImpl) {}
Atom(cif::Datablock &db, cif::Row &row); Atom(cif::Datablock &db, cif::Row &row);
// a special constructor to create symmetry copies // a special constructor to create symmetry copies
Atom(const Atom &rhs, const Point &symmmetry_location, const std::string &symmetry_operation); Atom(const Atom &rhs, const Point &symmmetry_location, const std::string &symmetry_operation);
~Atom(); explicit operator bool() const { return (bool)mImpl; }
explicit operator bool() const { return mImpl_ != nullptr; }
// return a copy of this atom, with data copied instead of referenced // return a copy of this atom, with data copied instead of referenced
Atom clone() const; Atom clone() const
{
auto copy = std::make_shared<AtomImpl>(*mImpl);
copy->mClone = true;
return Atom(copy);
}
Atom &operator=(const Atom &rhs) = default;
template <typename T>
T get_property(const std::string_view name) const;
Atom &operator=(const Atom &rhs); void set_property(const std::string_view name, const std::string &value)
{
mImpl->set_property(name, value);
}
template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
void property(const std::string_view name, const T &value)
{
set_property(name, std::to_string(value));
}
const std::string &id() const; const std::string &id() const { return mImpl->mID; }
AtomType type() const; AtomType type() const { return mImpl->mType; }
Point location() const; Point location() const { return mImpl->mLocation; }
void location(Point p); void location(Point p) { mImpl->moveTo(p); }
/// \brief Translate the position of this atom by \a t /// \brief Translate the position of this atom by \a t
void translate(Point t); void translate(Point t);
...@@ -92,46 +163,33 @@ class Atom ...@@ -92,46 +163,33 @@ class Atom
void rotate(Quaternion q); void rotate(Quaternion q);
// for direct access to underlying data, be careful! // for direct access to underlying data, be careful!
const cif::Row getRow() const; const cif::Row getRow() const { return mImpl->mRow; }
const cif::Row getRowAniso() const; const cif::Row getRowAniso() const;
// Atom symmetryCopy(const Point& d, const clipper::RTop_orth& rt); bool isSymmetryCopy() const { return mImpl->mSymmetryCopy; }
bool isSymmetryCopy() const; std::string symmetry() const { return mImpl->mSymmetryOperator; }
std::string symmetry() const;
// const clipper::RTop_orth& symop() const;
const Compound &comp() const; const Compound &comp() const { return mImpl->comp(); }
bool isWater() const; bool isWater() const { return mImpl->mCompID == "HOH" or mImpl->mCompID == "H2O" or mImpl->mCompID == "WAT"; }
int charge() const; int charge() const;
float uIso() const; float uIso() const;
bool getAnisoU(float anisou[6]) const; bool getAnisoU(float anisou[6]) const { return mImpl->getAnisoU(anisou); }
float occupancy() const; float occupancy() const;
template <typename T>
T property(const std::string_view name) const;
void property(const std::string_view name, const std::string &value);
template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
void property(const std::string_view name, const T &value)
{
property(name, std::to_string(value));
}
// specifications // specifications
const std::string& labelAtomID() const { return mAtomID; } const std::string& labelAtomID() const { return mImpl->mAtomID; }
const std::string& labelCompID() const { return mCompID; } const std::string& labelCompID() const { return mImpl->mCompID; }
const std::string& labelAsymID() const { return mAsymID; } const std::string& labelAsymID() const { return mImpl->mAsymID; }
std::string labelEntityID() const; std::string labelEntityID() const;
int labelSeqID() const { return mSeqID; } int labelSeqID() const { return mImpl->mSeqID; }
const std::string& labelAltID() const { return mAltID; } const std::string& labelAltID() const { return mImpl->mAltID; }
bool isAlternate() const { return not mAltID.empty(); } bool isAlternate() const { return not mImpl->mAltID.empty(); }
std::string authAtomID() const; std::string authAtomID() const;
std::string authCompID() const; std::string authCompID() const;
std::string authAsymID() const; std::string authAsymID() const;
const std::string& authSeqID() const { return mAuthSeqID; } const std::string& authSeqID() const { return mImpl->mAuthSeqID; }
std::string pdbxAuthInsCode() const; std::string pdbxAuthInsCode() const;
std::string pdbxAuthAltID() const; std::string pdbxAuthAltID() const;
...@@ -140,13 +198,6 @@ class Atom ...@@ -140,13 +198,6 @@ class Atom
bool operator==(const Atom &rhs) const; bool operator==(const Atom &rhs) const;
// // get clipper format Atom
// clipper::Atom toClipper() const;
// Radius calculation based on integrating the density until perc of electrons is found
void calculateRadius(float resHigh, float resLow, float perc);
float radius() const;
// access data in compound for this atom // access data in compound for this atom
// convenience routine // convenience routine
...@@ -158,16 +209,10 @@ class Atom ...@@ -158,16 +209,10 @@ class Atom
void swap(Atom &b) void swap(Atom &b)
{ {
std::swap(mImpl_, b.mImpl_); std::swap(mImpl, b.mImpl);
std::swap(mAtomID, b.mAtomID);
std::swap(mCompID, b.mCompID);
std::swap(mAsymID, b.mAsymID);
std::swap(mSeqID, b.mSeqID);
std::swap(mAltID, b.mAltID);
std::swap(mAuthSeqID, b.mAuthSeqID);
} }
int compare(const Atom &b) const; int compare(const Atom &b) const { return mImpl->compare(*b.mImpl); }
bool operator<(const Atom &rhs) const bool operator<(const Atom &rhs) const
{ {
...@@ -178,21 +223,30 @@ class Atom ...@@ -178,21 +223,30 @@ class Atom
private: private:
friend class Structure; friend class Structure;
void setID(int id); void setID(int id);
AtomImpl *impl(); std::shared_ptr<AtomImpl> mImpl;
const AtomImpl *impl() const; };
template <>
inline std::string Atom::get_property<std::string>(const std::string_view name) const
{
return mImpl->get_property(name);
}
struct AtomImpl *mImpl_; template <>
inline int Atom::get_property<int>(const std::string_view name) const
{
auto v = mImpl->get_property(name);
return v.empty() ? 0 : stoi(v);
}
// cached values template <>
std::string mAtomID; inline float Atom::get_property<float>(const std::string_view name) const
std::string mCompID; {
std::string mAsymID; return stof(mImpl->get_property(name));
int mSeqID; }
std::string mAltID;
std::string mAuthSeqID;
};
inline void swap(mmcif::Atom &a, mmcif::Atom &b) inline void swap(mmcif::Atom &a, mmcif::Atom &b)
{ {
...@@ -216,19 +270,16 @@ typedef std::vector<Atom> AtomView; ...@@ -216,19 +270,16 @@ typedef std::vector<Atom> AtomView;
class Residue class Residue
{ {
public: public:
// constructors should be private, but that's not possible for now (needed in emplace) // constructor
// constructor for waters
Residue(const Structure &structure, const std::string &compoundID,
const std::string &asymID, const std::string &authSeqID);
// constructor for a residue without a sequence number
Residue(const Structure &structure, const std::string &compoundID, Residue(const Structure &structure, const std::string &compoundID,
const std::string &asymID); const std::string &asymID, int seqID = 0, const std::string &authSeqID = {})
: mStructure(&structure)
// constructor for a residue with a sequence number , mCompoundID(compoundID)
Residue(const Structure &structure, const std::string &compoundID, , mAsymID(asymID)
const std::string &asymID, int seqID, const std::string &authSeqID); , mSeqID(seqID)
, mAuthSeqID(authSeqID)
{
}
Residue(const Residue &rhs) = delete; Residue(const Residue &rhs) = delete;
Residue &operator=(const Residue &rhs) = delete; Residue &operator=(const Residue &rhs) = delete;
......
src/Structure.cpp
/*- /*-
* SPDX-License-Identifier: BSD-2-Clause * SPDX-License-Identifier: BSD-2-Clause
* *
* Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
* *
* Redistribution and use in source and binary forms, with or without * Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met: * modification, are permitted provided that the following conditions are met:
* *
* 1. Redistributions of source code must retain the above copyright notice, this * 1. Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer * list of conditions and the following disclaimer
* 2. Redistributions in binary form must reproduce the above copyright notice, * 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation * this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution. * and/or other materials provided with the distribution.
* *
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
...@@ -202,443 +202,243 @@ void FileImpl::save(const std::filesystem::path &path) ...@@ -202,443 +202,243 @@ void FileImpl::save(const std::filesystem::path &path)
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// Atom // Atom
struct AtomImpl Atom::AtomImpl::AtomImpl(cif::Datablock &db, const std::string &id, cif::Row row)
: mDb(db)
, mID(id)
, mRefcount(1)
, mRow(row)
, mCompound(nullptr)
{ {
AtomImpl(const AtomImpl &i) prefetch();
: mDb(i.mDb) }
, mID(i.mID)
, mType(i.mType)
, mLocation(i.mLocation)
, mRefcount(1)
, mRow(i.mRow)
, mCachedRefs(i.mCachedRefs)
, mCompound(i.mCompound)
, mRadius(i.mRadius)
, mSymmetryCopy(i.mSymmetryCopy)
, mClone(true)
{
}
AtomImpl(cif::Datablock &db, const std::string &id)
: mDb(db)
, mID(id)
, mRefcount(1)
, mCompound(nullptr)
{
auto &cat = db["atom_site"];
mRow = cat[cif::Key("id") == mID];
}
AtomImpl(cif::Datablock &db, cif::Row &row)
: mDb(db)
, mID(row["id"].as<std::string>())
, mRefcount(1)
, mRow(row)
, mCompound(nullptr)
{
}
AtomImpl(cif::Datablock &db, const std::string &id, cif::Row row)
: mDb(db)
, mID(id)
, mRefcount(1)
, mRow(row)
, mCompound(nullptr)
{
}
AtomImpl(const AtomImpl &impl, const Point &loc, const std::string &sym_op) // constructor for a symmetry copy of an atom
: mDb(impl.mDb) Atom::AtomImpl::AtomImpl(const AtomImpl &impl, const Point &loc, const std::string &sym_op)
, mID(impl.mID) : mDb(impl.mDb)
, mType(impl.mType) , mID(impl.mID)
, mLocation(loc) , mType(impl.mType)
, mRefcount(1)
, mRow(impl.mRow)
, mCachedRefs(impl.mCachedRefs)
, mCompound(impl.mCompound)
, mRadius(impl.mRadius)
, mSymmetryCopy(true)
, mSymmetryOperator(sym_op)
{
}
void prefetch(std::string &atomID, std::string &compID, std::string &asymID, int &seqID, std::string &altID, std::string &authSeqID) , mAtomID(impl.mAtomID)
{ , mCompID(impl.mCompID)
// Prefetch some data , mAsymID(impl.mAsymID)
std::string symbol; , mSeqID(impl.mSeqID)
cif::tie(symbol, atomID, compID, asymID, seqID, altID, authSeqID) = , mAltID(impl.mAltID)
mRow.get("type_symbol", "label_atom_id", "label_comp_id", "label_asym_id", "label_seq_id", "label_alt_id", "auth_seq_id"); , mAuthSeqID(impl.mAuthSeqID)
if (symbol != "X") , mLocation(loc)
mType = AtomTypeTraits(symbol).type(); , mRefcount(1)
, mRow(impl.mRow)
, mCachedRefs(impl.mCachedRefs)
, mCompound(impl.mCompound)
, mSymmetryCopy(true)
, mSymmetryOperator(sym_op)
{
}
float x, y, z; void Atom::AtomImpl::prefetch()
cif::tie(x, y, z) = mRow.get("Cartn_x", "Cartn_y", "Cartn_z"); {
// Prefetch some data
std::string symbol;
cif::tie(symbol, mAtomID, mCompID, mAsymID, mSeqID, mAltID, mAuthSeqID) =
mRow.get("type_symbol", "label_atom_id", "label_comp_id", "label_asym_id", "label_seq_id", "label_alt_id", "auth_seq_id");
mLocation = Point(x, y, z); if (symbol != "X")
} mType = AtomTypeTraits(symbol).type();
void reference() float x, y, z;
{ cif::tie(x, y, z) = mRow.get("Cartn_x", "Cartn_y", "Cartn_z");
++mRefcount;
}
void release() mLocation = Point(x, y, z);
{ }
if (--mRefcount <= 0)
delete this;
}
bool getAnisoU(float anisou[6]) const int Atom::AtomImpl::compare(const AtomImpl &b) const
{ {
bool result = false; int d = mAsymID.compare(b.mAsymID);
if (d == 0)
d = mSeqID - b.mSeqID;
if (d == 0)
d = mAtomID.compare(b.mAtomID);
if (d == 0)
d = mAuthSeqID.compare(b.mAuthSeqID);
auto cat = mDb.get("atom_site_anisotrop"); return d;
if (cat) }
{
try
{
auto r = cat->find1(cif::Key("id") == mID);
cif::tie(anisou[0], anisou[1], anisou[2], anisou[3], anisou[4], anisou[5]) =
r.get("U[1][1]", "U[1][2]", "U[1][3]", "U[2][2]", "U[2][3]", "U[3][3]");
result = true;
}
catch(const std::exception& e)
{
}
}
return result; bool Atom::AtomImpl::getAnisoU(float anisou[6]) const
} {
bool result = false;
void moveTo(const Point &p) auto cat = mDb.get("atom_site_anisotrop");
if (cat)
{ {
assert(not mSymmetryCopy); try
if (mSymmetryCopy)
throw std::runtime_error("Moving symmetry copy");
if (not mClone)
{ {
property("Cartn_x", std::to_string(p.getX())); auto r = cat->find1(cif::Key("id") == mID);
property("Cartn_y", std::to_string(p.getY())); cif::tie(anisou[0], anisou[1], anisou[2], anisou[3], anisou[4], anisou[5]) =
property("Cartn_z", std::to_string(p.getZ())); r.get("U[1][1]", "U[1][2]", "U[1][3]", "U[2][2]", "U[2][3]", "U[3][3]");
result = true;
} }
catch (const std::exception &e)
// boost::format kPosFmt("%.3f");
//
// mRow["Cartn_x"] = (kPosFmt % p.getX()).str();
// mRow["Cartn_y"] = (kPosFmt % p.getY()).str();
// mRow["Cartn_z"] = (kPosFmt % p.getZ()).str();
mLocation = p;
}
const Compound &comp() const
{
if (mCompound == nullptr)
{ {
std::string compID;
cif::tie(compID) = mRow.get("label_comp_id");
mCompound = CompoundFactory::instance().create(compID);
if (cif::VERBOSE and mCompound == nullptr)
std::cerr << "Compound not found: '" << compID << '\'' << std::endl;
} }
if (mCompound == nullptr)
throw std::runtime_error("no compound");
return *mCompound;
} }
float radius() const return result;
{ }
return mRadius;
}
const std::string property(const std::string_view name) const
{
for (auto &&[tag, ref] : mCachedRefs)
{
if (tag == name)
return ref.as<std::string>();
}
mCachedRefs.emplace_back(name, mRow[name]); void Atom::AtomImpl::moveTo(const Point &p)
return std::get<1>(mCachedRefs.back()).as<std::string>(); {
} assert(not mSymmetryCopy);
if (mSymmetryCopy)
throw std::runtime_error("Moving symmetry copy");
void property(const std::string_view name, const std::string &value) if (not mClone)
{ {
for (auto &&[tag, ref] : mCachedRefs) set_property("Cartn_x", std::to_string(p.getX()));
{ set_property("Cartn_y", std::to_string(p.getY()));
if (tag != name) set_property("Cartn_z", std::to_string(p.getZ()));
continue;
ref = value;
return;
}
mCachedRefs.emplace_back(name, mRow[name]);
std::get<1>(mCachedRefs.back()) = value;
} }
const cif::Datablock &mDb; // boost::format kPosFmt("%.3f");
std::string mID; //
AtomType mType; // mRow["Cartn_x"] = (kPosFmt % p.getX()).str();
// mRow["Cartn_y"] = (kPosFmt % p.getY()).str();
Point mLocation; // mRow["Cartn_z"] = (kPosFmt % p.getZ()).str();
int mRefcount;
cif::Row mRow;
mutable std::vector<std::tuple<std::string,cif::detail::ItemReference>> mCachedRefs;
mutable const Compound *mCompound = nullptr; mLocation = p;
float mRadius = std::nanf("4");
bool mSymmetryCopy = false;
bool mClone = false;
std::string mSymmetryOperator = "1_555";
// clipper::RTop_orth mRTop;
// Point mD;
// int32_t mRTix;
};
//Atom::Atom(const File& f, const std::string& id)
// : mImpl(new AtomImpl(f, id))
//{
//}
//
Atom::Atom()
: Atom(nullptr)
{
}
Atom::Atom(AtomImpl *impl)
: mImpl_(impl)
{
if (mImpl_)
mImpl_->prefetch(mAtomID, mCompID, mAsymID, mSeqID, mAltID, mAuthSeqID);
} }
Atom::Atom(cif::Datablock &db, cif::Row &row) const Compound &Atom::AtomImpl::comp() const
: Atom(new AtomImpl(db, row))
{ {
} if (mCompound == nullptr)
{
std::string compID;
cif::tie(compID) = mRow.get("label_comp_id");
Atom::Atom(const Atom &rhs, const Point &loc, const std::string &sym_op) mCompound = CompoundFactory::instance().create(compID);
: Atom(new AtomImpl(*rhs.mImpl_, loc, sym_op))
{
}
Atom::Atom(const Atom &rhs) if (cif::VERBOSE and mCompound == nullptr)
: mImpl_(rhs.mImpl_) std::cerr << "Compound not found: '" << compID << '\'' << std::endl;
, mAtomID(rhs.mAtomID) }
, mCompID(rhs.mCompID)
, mAsymID(rhs.mAsymID)
, mSeqID(rhs.mSeqID)
, mAltID(rhs.mAltID)
, mAuthSeqID(rhs.mAuthSeqID)
{
if (mImpl_)
mImpl_->reference();
}
Atom::~Atom() if (mCompound == nullptr)
{ throw std::runtime_error("no compound");
if (mImpl_)
mImpl_->release();
}
AtomImpl *Atom::impl() return *mCompound;
{
if (mImpl_ == nullptr)
throw std::runtime_error("atom is not set");
return mImpl_;
} }
const AtomImpl *Atom::impl() const const std::string Atom::AtomImpl::get_property(const std::string_view name) const
{ {
if (mImpl_ == nullptr) for (auto &&[tag, ref] : mCachedRefs)
throw std::runtime_error("atom is not set"); {
return mImpl_; if (tag == name)
} return ref.as<std::string>();
}
Atom Atom::clone() const mCachedRefs.emplace_back(name, mRow[name]);
{ return std::get<1>(mCachedRefs.back()).as<std::string>();
return Atom(mImpl_ ? new AtomImpl(*mImpl_) : nullptr);
} }
Atom &Atom::operator=(const Atom &rhs) void Atom::AtomImpl::set_property(const std::string_view name, const std::string &value)
{ {
if (this != &rhs) for (auto &&[tag, ref] : mCachedRefs)
{ {
if (mImpl_) if (tag != name)
mImpl_->release(); continue;
mImpl_ = rhs.mImpl_;
mAtomID = rhs.mAtomID;
mCompID = rhs.mCompID;
mAsymID = rhs.mAsymID;
mSeqID = rhs.mSeqID;
mAltID = rhs.mAltID;
mAuthSeqID = rhs.mAuthSeqID;
if (mImpl_) ref = value;
mImpl_->reference(); return;
} }
return *this; mCachedRefs.emplace_back(name, mRow[name]);
std::get<1>(mCachedRefs.back()) = value;
} }
const cif::Row Atom::getRow() const Atom::Atom(cif::Datablock &db, cif::Row &row)
: Atom(std::make_shared<AtomImpl>(db, row["id"].as<std::string>(), row))
{
}
Atom::Atom(const Atom &rhs, const Point &loc, const std::string &sym_op)
: Atom(std::make_shared<AtomImpl>(*rhs.mImpl, loc, sym_op))
{ {
return mImpl_->mRow;
} }
const cif::Row Atom::getRowAniso() const const cif::Row Atom::getRowAniso() const
{ {
auto &db = mImpl_->mDb; auto &db = mImpl->mDb;
auto cat = db.get("atom_site_anisotrop"); auto cat = db.get("atom_site_anisotrop");
if (not cat) if (not cat)
return {}; return {};
else else
return cat->find1(cif::Key("id") == mImpl_->mID); return cat->find1(cif::Key("id") == mImpl->mID);
}
template <>
std::string Atom::property<std::string>(const std::string_view name) const
{
return impl()->property(name);
}
template <>
int Atom::property<int>(const std::string_view name) const
{
auto v = impl()->property(name);
return v.empty() ? 0 : stoi(v);
}
template <>
float Atom::property<float>(const std::string_view name) const
{
return stof(impl()->property(name));
}
void Atom::property(const std::string_view name, const std::string &value)
{
impl()->property(name, value);
}
const std::string &Atom::id() const
{
return impl()->mID;
}
AtomType Atom::type() const
{
return impl()->mType;
}
int Atom::charge() const
{
return property<int>("pdbx_formal_charge");
} }
float Atom::uIso() const float Atom::uIso() const
{ {
float result; float result;
if (not property<std::string>("U_iso_or_equiv").empty()) if (not get_property<std::string>("U_iso_or_equiv").empty())
result = property<float>("U_iso_or_equiv"); result = get_property<float>("U_iso_or_equiv");
else if (not property<std::string>("B_iso_or_equiv").empty()) else if (not get_property<std::string>("B_iso_or_equiv").empty())
result = property<float>("B_iso_or_equiv") / static_cast<float>(8 * kPI * kPI); result = get_property<float>("B_iso_or_equiv") / static_cast<float>(8 * kPI * kPI);
else else
throw std::runtime_error("Missing B_iso or U_iso"); throw std::runtime_error("Missing B_iso or U_iso");
return result; return result;
} }
bool Atom::getAnisoU(float anisou[6]) const std::string Atom::labelID() const
{
return impl()->getAnisoU(anisou);
}
float Atom::occupancy() const
{ {
return property<float>("occupancy"); return mImpl->mCompID + '_' + mImpl->mAsymID + '_' + std::to_string(mImpl->mSeqID) + ':' + mImpl->mAtomID;
} }
// const std::string& Atom::labelAtomID() const std::string Atom::pdbID() const
// {
// return impl()->mAtomID;
// }
// const std::string& Atom::labelCompID() const
// {
// return impl()->mCompID;
// }
// const std::string& Atom::labelAsymID() const
// {
// return impl()->mAsymID;
// }
std::string Atom::labelEntityID() const
{ {
return property<std::string>("label_entity_id"); return get_property<std::string>("auth_comp_id") + '_' +
get_property<std::string>("auth_asym_id") + '_' +
get_property<std::string>("auth_seq_id") +
get_property<std::string>("pdbx_PDB_ins_code");
} }
std::string Atom::authAsymID() const int Atom::charge() const
{ {
return property<std::string>("auth_asym_id"); return get_property<int>("pdbx_formal_charge");
} }
std::string Atom::authAtomID() const float Atom::occupancy() const
{ {
return property<std::string>("auth_atom_id"); return get_property<float>("occupancy");
} }
std::string Atom::pdbxAuthAltID() const std::string Atom::labelEntityID() const
{ {
return property<std::string>("pdbx_auth_alt_id"); return get_property<std::string>("label_entity_id");
} }
std::string Atom::pdbxAuthInsCode() const std::string Atom::authAtomID() const
{ {
return property<std::string>("pdbx_PDB_ins_code"); return get_property<std::string>("auth_atom_id");
} }
std::string Atom::authCompID() const std::string Atom::authCompID() const
{ {
return property<std::string>("auth_comp_id"); return get_property<std::string>("auth_comp_id");
} }
std::string Atom::labelID() const std::string Atom::authAsymID() const
{
return mCompID + '_' + mAsymID + '_' + std::to_string(mSeqID) + ':' + mAtomID;
}
std::string Atom::pdbID() const
{ {
return property<std::string>("auth_comp_id") + '_' + return get_property<std::string>("auth_asym_id");
property<std::string>("auth_asym_id") + '_' +
property<std::string>("auth_seq_id") +
property<std::string>("pdbx_PDB_ins_code");
} }
Point Atom::location() const std::string Atom::pdbxAuthInsCode() const
{ {
return impl()->mLocation; return get_property<std::string>("pdbx_PDB_ins_code");
} }
void Atom::location(Point p) std::string Atom::pdbxAuthAltID() const
{ {
impl()->moveTo(p); return get_property<std::string>("pdbx_auth_alt_id");
} }
void Atom::translate(Point t) void Atom::translate(Point t)
...@@ -655,89 +455,15 @@ void Atom::rotate(Quaternion q) ...@@ -655,89 +455,15 @@ void Atom::rotate(Quaternion q)
location(loc); location(loc);
} }
// Atom Atom::symmetryCopy(const Point& d, const clipper::RTop_orth& rt)
// {
// return Atom(new AtomImpl(*impl(), d, rt));
// }
// bool Atom::isSymmetryCopy() const
// {
// return impl()->mSymmetryCopy;
// }
// std::string Atom::symmetry() const
// {
// return clipper::Symop(impl()->mRTop).format() + "\n" + impl()->mRTop.format();
// }
bool Atom::isSymmetryCopy() const
{
return mImpl_->mSymmetryCopy;
}
std::string Atom::symmetry() const
{
return mImpl_->mSymmetryOperator;
}
// const clipper::RTop_orth& Atom::symop() const
// {
// return impl()->mRTop;
// }
const Compound &Atom::comp() const
{
return impl()->comp();
}
bool Atom::isWater() const
{
(void)impl();
return mCompID == "HOH" or mCompID == "H2O" or mCompID == "WAT";
}
bool Atom::operator==(const Atom &rhs) const bool Atom::operator==(const Atom &rhs) const
{ {
return impl() == rhs.impl() or return mImpl == rhs.mImpl or
(&impl()->mDb == &rhs.impl()->mDb and impl()->mID == rhs.impl()->mID); (&mImpl->mDb == &rhs.mImpl->mDb and mImpl->mID == rhs.mImpl->mID);
}
// clipper::Atom Atom::toClipper() const
// {
// return impl()->toClipper();
// }
// void Atom::calculateRadius(float resHigh, float resLow, float perc)
// {
// AtomShape shape(*this, resHigh, resLow, false);
// impl()->mRadius = shape.radius();
// // verbose
// if (cif::VERBOSE > 1)
// cout << "Calculated radius for " << AtomTypeTraits(impl()->mType).name() << " with charge " << charge() << " is " << impl()->mRadius << std::endl;
// }
float Atom::radius() const
{
return impl()->mRadius;
}
int Atom::compare(const Atom &b) const
{
int d = mAsymID.compare(b.mAsymID);
if (d == 0)
d = mSeqID - b.mSeqID;
if (d == 0)
d = mAtomID.compare(b.mAtomID);
if (d == 0)
d = mAuthSeqID.compare(b.mAuthSeqID);
return d;
} }
void Atom::setID(int id) void Atom::setID(int id)
{ {
impl()->mID = std::to_string(id); mImpl->mID = std::to_string(id);
} }
std::ostream &operator<<(std::ostream &os, const Atom &atom) std::ostream &operator<<(std::ostream &os, const Atom &atom)
...@@ -755,66 +481,6 @@ std::ostream &operator<<(std::ostream &os, const Atom &atom) ...@@ -755,66 +481,6 @@ std::ostream &operator<<(std::ostream &os, const Atom &atom)
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// residue // residue
// First constructor used to be for waters only, but now accepts sugars as well.
Residue::Residue(const Structure &structure, const std::string &compoundID,
const std::string &asymID, const std::string &authSeqID)
: mStructure(&structure)
, mCompoundID(compoundID)
, mAsymID(asymID)
, mAuthSeqID(authSeqID)
{
// for (auto &a : mStructure->atoms())
// {
// if (a.labelAsymID() != mAsymID or
// a.labelCompID() != mCompoundID)
// continue;
// if (compoundID == "HOH")
// {
// if (not mAuthSeqID.empty() and a.authSeqID() != mAuthSeqID)
// continue;
// }
// else
// {
// if (mSeqID > 0 and a.labelSeqID() != mSeqID)
// continue;
// }
// mAtoms.push_back(a);
// }
// assert(not mAtoms.empty());
}
Residue::Residue(const Structure &structure, const std::string &compoundID, const std::string &asymID)
: Residue(structure, compoundID, asymID, 0, {})
{
}
Residue::Residue(const Structure &structure, const std::string &compoundID,
const std::string &asymID, int seqID, const std::string &authSeqID)
: mStructure(&structure)
, mCompoundID(compoundID)
, mAsymID(asymID)
, mSeqID(seqID)
, mAuthSeqID(authSeqID)
{
assert(mCompoundID != "HOH");
// for (auto &a : mStructure->atoms())
// {
// if (mSeqID > 0 and a.labelSeqID() != mSeqID)
// continue;
// if (a.labelAsymID() != mAsymID or
// a.labelCompID() != mCompoundID)
// continue;
// mAtoms.push_back(a);
// }
}
Residue::Residue(Residue &&rhs) Residue::Residue(Residue &&rhs)
: mStructure(rhs.mStructure) : mStructure(rhs.mStructure)
, mCompoundID(std::move(rhs.mCompoundID)) , mCompoundID(std::move(rhs.mCompoundID))
...@@ -823,13 +489,13 @@ Residue::Residue(Residue &&rhs) ...@@ -823,13 +489,13 @@ Residue::Residue(Residue &&rhs)
, mAuthSeqID(rhs.mAuthSeqID) , mAuthSeqID(rhs.mAuthSeqID)
, mAtoms(std::move(rhs.mAtoms)) , mAtoms(std::move(rhs.mAtoms))
{ {
//std::cerr << "move constructor residue" << std::endl; // std::cerr << "move constructor residue" << std::endl;
rhs.mStructure = nullptr; rhs.mStructure = nullptr;
} }
Residue &Residue::operator=(Residue &&rhs) Residue &Residue::operator=(Residue &&rhs)
{ {
//std::cerr << "move assignment residue" << std::endl; // std::cerr << "move assignment residue" << std::endl;
mStructure = rhs.mStructure; mStructure = rhs.mStructure;
rhs.mStructure = nullptr; rhs.mStructure = nullptr;
mCompoundID = std::move(rhs.mCompoundID); mCompoundID = std::move(rhs.mCompoundID);
...@@ -843,7 +509,7 @@ Residue &Residue::operator=(Residue &&rhs) ...@@ -843,7 +509,7 @@ Residue &Residue::operator=(Residue &&rhs)
Residue::~Residue() Residue::~Residue()
{ {
//std::cerr << "~Residue" << std::endl; // std::cerr << "~Residue" << std::endl;
} }
std::string Residue::entityID() const std::string Residue::entityID() const
...@@ -1598,7 +1264,7 @@ File::~File() ...@@ -1598,7 +1264,7 @@ File::~File()
delete mImpl; delete mImpl;
} }
cif::Datablock& File::createDatablock(const std::string_view name) cif::Datablock &File::createDatablock(const std::string_view name)
{ {
auto db = new cif::Datablock(name); auto db = new cif::Datablock(name);
...@@ -1692,7 +1358,7 @@ void Structure::loadAtomsForModel(StructureOpenOptions options) ...@@ -1692,7 +1358,7 @@ void Structure::loadAtomsForModel(StructureOpenOptions options)
if ((options bitand StructureOpenOptions::SkipHydrogen) and type_symbol == "H") if ((options bitand StructureOpenOptions::SkipHydrogen) and type_symbol == "H")
continue; continue;
mAtoms.emplace_back(new AtomImpl(*db, id, a)); mAtoms.emplace_back(std::make_shared<Atom::AtomImpl>(*db, id, a));
} }
} }
...@@ -1736,7 +1402,7 @@ void Structure::loadData() ...@@ -1736,7 +1402,7 @@ void Structure::loadData()
r.get("asym_id", "mon_id", "pdb_seq_num"); r.get("asym_id", "mon_id", "pdb_seq_num");
if (monID == "HOH") if (monID == "HOH")
mNonPolymers.emplace_back(*this, monID, asymID, pdbSeqNum); mNonPolymers.emplace_back(*this, monID, asymID, 0, pdbSeqNum);
else if (mNonPolymers.empty() or mNonPolymers.back().asymID() != asymID) else if (mNonPolymers.empty() or mNonPolymers.back().asymID() != asymID)
mNonPolymers.emplace_back(*this, monID, asymID); mNonPolymers.emplace_back(*this, monID, asymID);
} }
...@@ -1749,28 +1415,28 @@ void Structure::loadData() ...@@ -1749,28 +1415,28 @@ void Structure::loadData()
cif::tie(asymID, monID, num) = cif::tie(asymID, monID, num) =
r.get("asym_id", "mon_id", "num"); r.get("asym_id", "mon_id", "num");
mBranchResidues.emplace_back(*this, monID, asymID, num); mBranchResidues.emplace_back(*this, monID, asymID, 0, num);
} }
// place atoms in residues // place atoms in residues
using key_type = std::tuple<std::string,std::string,int>; using key_type = std::tuple<std::string, std::string, int>;
using map_type = std::map<key_type,Residue*>; using map_type = std::map<key_type, Residue *>;
map_type resMap; map_type resMap;
for (auto &poly : mPolymers) for (auto &poly : mPolymers)
{ {
for (auto &res : poly) for (auto &res : poly)
resMap[{ res.asymID(), res.compoundID(), res.seqID() }] = &res; resMap[{res.asymID(), res.compoundID(), res.seqID()}] = &res;
} }
for (auto &res : mNonPolymers) for (auto &res : mNonPolymers)
resMap[{ res.asymID(), res.compoundID(), (res.isWater() ? std::stoi(res.mAuthSeqID) : res.seqID()) }] = &res; resMap[{res.asymID(), res.compoundID(), (res.isWater() ? std::stoi(res.mAuthSeqID) : res.seqID())}] = &res;
for (auto &res : mBranchResidues) for (auto &res : mBranchResidues)
resMap[{ res.asymID(), res.compoundID(), res.seqID() }] = &res; resMap[{res.asymID(), res.compoundID(), res.seqID()}] = &res;
for (auto &atom: mAtoms) for (auto &atom : mAtoms)
{ {
key_type k(atom.labelAsymID(), atom.labelCompID(), atom.isWater() ? std::stoi(atom.authSeqID()) : atom.labelSeqID()); key_type k(atom.labelAsymID(), atom.labelCompID(), atom.isWater() ? std::stoi(atom.authSeqID()) : atom.labelSeqID());
auto ri = resMap.find(k); auto ri = resMap.find(k);
...@@ -1784,7 +1450,7 @@ void Structure::loadData() ...@@ -1784,7 +1450,7 @@ void Structure::loadData()
continue; continue;
} }
ri->second->addAtom(atom); ri->second->addAtom(atom);
} }
} }
...@@ -1836,7 +1502,7 @@ AtomView Structure::waters() const ...@@ -1836,7 +1502,7 @@ AtomView Structure::waters() const
for (auto &a : mAtoms) for (auto &a : mAtoms)
{ {
if (a.property<std::string>("label_entity_id") == waterEntityID) if (a.get_property<std::string>("label_entity_id") == waterEntityID)
result.push_back(a); result.push_back(a);
} }
...@@ -1909,7 +1575,7 @@ Atom Structure::getAtomByPositionAndType(Point p, std::string_view type, std::st ...@@ -1909,7 +1575,7 @@ Atom Structure::getAtomByPositionAndType(Point p, std::string_view type, std::st
if (a.labelCompID() != res_type) if (a.labelCompID() != res_type)
continue; continue;
if (a.labelAtomID() != type) if (a.labelAtomID() != type)
continue; continue;
...@@ -1965,7 +1631,7 @@ const Residue &Structure::getResidue(const std::string &asymID, const std::strin ...@@ -1965,7 +1631,7 @@ const Residue &Structure::getResidue(const std::string &asymID, const std::strin
Residue &Structure::getResidue(const std::string &asymID, const std::string &compID, int seqID) Residue &Structure::getResidue(const std::string &asymID, const std::string &compID, int seqID)
{ {
return const_cast<Residue&>(const_cast<Structure const&>(*this).getResidue(asymID, compID, seqID)); return const_cast<Residue &>(const_cast<Structure const &>(*this).getResidue(asymID, compID, seqID));
} }
const Residue &Structure::getResidue(const std::string &asymID) const const Residue &Structure::getResidue(const std::string &asymID) const
...@@ -1983,7 +1649,7 @@ const Residue &Structure::getResidue(const std::string &asymID) const ...@@ -1983,7 +1649,7 @@ const Residue &Structure::getResidue(const std::string &asymID) const
Residue &Structure::getResidue(const std::string &asymID) Residue &Structure::getResidue(const std::string &asymID)
{ {
return const_cast<Residue&>(const_cast<Structure const&>(*this).getResidue(asymID)); return const_cast<Residue &>(const_cast<Structure const &>(*this).getResidue(asymID));
} }
File &Structure::getFile() const File &Structure::getFile() const
...@@ -2169,8 +1835,7 @@ std::string Structure::insertCompound(const std::string &compoundID, bool isEnti ...@@ -2169,8 +1835,7 @@ std::string Structure::insertCompound(const std::string &compoundID, bool isEnti
auto r = chemComp.find(cif::Key("id") == compoundID); auto r = chemComp.find(cif::Key("id") == compoundID);
if (r.empty()) if (r.empty())
{ {
chemComp.emplace({ chemComp.emplace({{"id", compoundID},
{"id", compoundID},
{"name", compound->name()}, {"name", compound->name()},
{"formula", compound->formula()}, {"formula", compound->formula()},
{"formula_weight", compound->formulaWeight()}, {"formula_weight", compound->formulaWeight()},
...@@ -2186,19 +1851,17 @@ std::string Structure::insertCompound(const std::string &compoundID, bool isEnti ...@@ -2186,19 +1851,17 @@ std::string Structure::insertCompound(const std::string &compoundID, bool isEnti
{ {
entity_id = pdbxEntityNonpoly.find1<std::string>("comp_id"_key == compoundID, "entity_id"); entity_id = pdbxEntityNonpoly.find1<std::string>("comp_id"_key == compoundID, "entity_id");
} }
catch(const std::exception& ex) catch (const std::exception &ex)
{ {
auto &entity = db["entity"]; auto &entity = db["entity"];
entity_id = entity.getUniqueID(""); entity_id = entity.getUniqueID("");
entity.emplace({ entity.emplace({{"id", entity_id},
{"id", entity_id},
{"type", "non-polymer"}, {"type", "non-polymer"},
{"pdbx_description", compound->name()}, {"pdbx_description", compound->name()},
{"formula_weight", compound->formulaWeight()}}); {"formula_weight", compound->formulaWeight()}});
pdbxEntityNonpoly.emplace({ pdbxEntityNonpoly.emplace({{"entity_id", entity_id},
{"entity_id", entity_id},
{"name", compound->name()}, {"name", compound->name()},
{"comp_id", compoundID}}); {"comp_id", compoundID}});
} }
...@@ -2312,7 +1975,7 @@ void Structure::swapAtoms(Atom &a1, Atom &a2) ...@@ -2312,7 +1975,7 @@ void Structure::swapAtoms(Atom &a1, Atom &a2)
auto l2 = r2["label_atom_id"]; auto l2 = r2["label_atom_id"];
l1.swap(l2); l1.swap(l2);
std::swap(a1.mAtomID, a2.mAtomID); std::swap(a1.mImpl->mAtomID, a2.mImpl->mAtomID);
auto l3 = r1["auth_atom_id"]; auto l3 = r1["auth_atom_id"];
auto l4 = r2["auth_atom_id"]; auto l4 = r2["auth_atom_id"];
...@@ -2442,47 +2105,43 @@ std::string Structure::createNonpoly(const std::string &entity_id, const std::ve ...@@ -2442,47 +2105,43 @@ std::string Structure::createNonpoly(const std::string &entity_id, const std::ve
auto &struct_asym = db["struct_asym"]; auto &struct_asym = db["struct_asym"];
std::string asym_id = struct_asym.getUniqueID(); std::string asym_id = struct_asym.getUniqueID();
struct_asym.emplace({ struct_asym.emplace({{"id", asym_id},
{ "id", asym_id }, {"pdbx_blank_PDB_chainid_flag", "N"},
{ "pdbx_blank_PDB_chainid_flag", "N" }, {"pdbx_modified", "N"},
{ "pdbx_modified", "N" }, {"entity_id", entity_id},
{ "entity_id", entity_id }, {"details", "?"}});
{ "details", "?" }
});
std::string comp_id = db["pdbx_entity_nonpoly"].find1<std::string>("entity_id"_key == entity_id, "comp_id"); std::string comp_id = db["pdbx_entity_nonpoly"].find1<std::string>("entity_id"_key == entity_id, "comp_id");
auto &atom_site = db["atom_site"]; auto &atom_site = db["atom_site"];
for (auto& atom: atoms) for (auto &atom : atoms)
{ {
auto atom_id = atom_site.getUniqueID(""); auto atom_id = atom_site.getUniqueID("");
auto &&[row, inserted ] = atom_site.emplace({ auto &&[row, inserted] = atom_site.emplace({{"group_PDB", atom.get_property<std::string>("group_PDB")},
{ "group_PDB", atom.property<std::string>("group_PDB") }, {"id", atom_id},
{ "id", atom_id }, {"type_symbol", atom.get_property<std::string>("type_symbol")},
{ "type_symbol", atom.property<std::string>("type_symbol") }, {"label_atom_id", atom.get_property<std::string>("label_atom_id")},
{ "label_atom_id", atom.property<std::string>("label_atom_id") }, {"label_alt_id", atom.get_property<std::string>("label_alt_id")},
{ "label_alt_id", atom.property<std::string>("label_alt_id") }, {"label_comp_id", comp_id},
{ "label_comp_id", comp_id }, {"label_asym_id", asym_id},
{ "label_asym_id", asym_id }, {"label_entity_id", entity_id},
{ "label_entity_id", entity_id }, {"label_seq_id", "."},
{ "label_seq_id", "." }, {"pdbx_PDB_ins_code", ""},
{ "pdbx_PDB_ins_code", "" }, {"Cartn_x", atom.get_property<std::string>("Cartn_x")},
{ "Cartn_x", atom.property<std::string>("Cartn_x") }, {"Cartn_y", atom.get_property<std::string>("Cartn_y")},
{ "Cartn_y", atom.property<std::string>("Cartn_y") }, {"Cartn_z", atom.get_property<std::string>("Cartn_z")},
{ "Cartn_z", atom.property<std::string>("Cartn_z") }, {"occupancy", atom.get_property<std::string>("occupancy")},
{ "occupancy", atom.property<std::string>("occupancy") }, {"B_iso_or_equiv", atom.get_property<std::string>("B_iso_or_equiv")},
{ "B_iso_or_equiv", atom.property<std::string>("B_iso_or_equiv") }, {"pdbx_formal_charge", atom.get_property<std::string>("pdbx_formal_charge")},
{ "pdbx_formal_charge", atom.property<std::string>("pdbx_formal_charge") }, {"auth_seq_id", ""},
{ "auth_seq_id", "" }, {"auth_comp_id", comp_id},
{ "auth_comp_id", comp_id }, {"auth_asym_id", asym_id},
{ "auth_asym_id", asym_id }, {"auth_atom_id", atom.get_property<std::string>("label_atom_id")},
{ "auth_atom_id", atom.property<std::string>("label_atom_id") }, {"pdbx_PDB_model_num", 1}});
{ "pdbx_PDB_model_num", 1 }
}); mAtoms.emplace_back(std::make_shared<Atom::AtomImpl>(db, atom_id, row));
mAtoms.emplace_back(new AtomImpl(db, atom_id, row));
} }
mNonPolymers.emplace_back(*this, comp_id, asym_id); mNonPolymers.emplace_back(*this, comp_id, asym_id);
...@@ -2563,7 +2222,7 @@ void Structure::cleanupEmptyCategories() ...@@ -2563,7 +2222,7 @@ void Structure::cleanupEmptyCategories()
{ {
// is this correct? // is this correct?
std::set<std::string> asym_ids; std::set<std::string> asym_ids;
for (const auto &[ asym_id ] : db["pdbx_branch_scheme"].find<std::string>("entity_id"_key == id, "asym_id")) for (const auto &[asym_id] : db["pdbx_branch_scheme"].find<std::string>("entity_id"_key == id, "asym_id"))
asym_ids.insert(asym_id); asym_ids.insert(asym_id);
count = asym_ids.size(); count = asym_ids.size();
} }
...@@ -2574,13 +2233,13 @@ void Structure::cleanupEmptyCategories() ...@@ -2574,13 +2233,13 @@ void Structure::cleanupEmptyCategories()
void Structure::translate(Point t) void Structure::translate(Point t)
{ {
for (auto& a: mAtoms) for (auto &a : mAtoms)
a.translate(t); a.translate(t);
} }
void Structure::rotate(Quaternion q) void Structure::rotate(Quaternion q)
{ {
for (auto& a: mAtoms) for (auto &a : mAtoms)
a.rotate(q); a.rotate(q);
} }
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment