Skip to content

L
libcifpp
Commit d4de02a2 by maarten
Options

backup 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@251 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent b60face9
Showing with 449 additions and 498 deletions
include/cif++/AtomShape.h
......@@ -4,7 +4,7 @@
#include "cif++/Structure.h"
namespace libcif
namespace mmcif
{
// --------------------------------------------------------------------
......
include/cif++/AtomType.h
......@@ -7,7 +7,7 @@
#include <boost/filesystem/operations.hpp>
#include <boost/math/quaternion.hpp>
namespace libcif
namespace mmcif
{
enum AtomType : uint8
......
include/cif++/BondMap.h
......@@ -6,7 +6,7 @@
#include "cif++/Structure.h"
namespace libcif
namespace mmcif
{
class BondMap
......
include/cif++/Compound.h
......@@ -10,7 +10,7 @@
#include "cif++/AtomType.h"
#include "cif++/Cif++.h"
namespace libcif
namespace mmcif
{
// --------------------------------------------------------------------
......@@ -24,32 +24,33 @@ namespace libcif
class Composition;
class Entity;
class Compound;
class Link;
struct CompoundAtom;
enum BondType { singleBond, doubleBond, tripleBond, delocalizedBond };
// --------------------------------------------------------------------
// struct containing information about an atom in a chemical compound
// This information comes from the CCP4 monomer library.
struct CompoundAtom
{
std::string id;
AtomType typeSymbol;
std::string typeEnergy;
float partialCharge;
std::string id;
AtomType typeSymbol;
std::string typeEnergy;
float partialCharge;
};
// --------------------------------------------------------------------
// struct containing information about the bonds
// This information comes from the CCP4 monomer library.
enum CompoundBondType { singleBond, doubleBond, tripleBond, delocalizedBond };
struct CompoundBond
{
std::string atomID[2];
CompoundBondType type;
float distance;
float esd;
std::string atomID[2];
BondType type;
float distance;
float esd;
};
// --------------------------------------------------------------------
......@@ -58,9 +59,9 @@ struct CompoundBond
struct CompoundAngle
{
std::string atomID[3];
float angle;
float esd;
std::string atomID[3];
float angle;
float esd;
};
// --------------------------------------------------------------------
......@@ -80,7 +81,7 @@ struct CompoundPlane
enum ChiralVolumeSign { negativ, positiv, both };
struct ChiralCentre
struct CompoundChiralCentre
{
std::string id;
std::string atomIDCentre;
......@@ -123,7 +124,8 @@ class Compound
std::vector<CompoundAtom> atoms() const { return mAtoms; }
std::vector<CompoundBond> bonds() const { return mBonds; }
std::vector<CompoundAngle> angles() const { return mAngles; }
std::vector<ChiralCentre> chiralCentres() const { return mChiralCentres; }
std::vector<CompoundChiralCentre> chiralCentres() const
{ return mChiralCentres; }
std::vector<CompoundPlane> planes() const { return mPlanes; }
CompoundAtom getAtomById(const std::string& atomId) const;
......@@ -155,12 +157,96 @@ class Compound
std::vector<CompoundAtom> mAtoms;
std::vector<CompoundBond> mBonds;
std::vector<CompoundAngle> mAngles;
std::vector<ChiralCentre> mChiralCentres;
std::vector<CompoundChiralCentre>
mChiralCentres;
std::vector<CompoundPlane> mPlanes;
};
// --------------------------------------------------------------------
// Factory class for Compound objects
// struct containing information about the bonds
// This information comes from the CCP4 monomer library.
struct LinkAtom
{
int compID;
std::string atomID;
};
struct LinkBond
{
LinkAtom atom[2];
BondType type;
float distance;
float esd;
};
// --------------------------------------------------------------------
// struct containing information about the bond-angles
// This information comes from the CCP4 monomer library.
struct LinkAngle
{
LinkAtom atom[3];
float angle;
float esd;
};
// --------------------------------------------------------------------
// struct containing information about the bond-angles
// This information comes from the CCP4 monomer library.
struct LinkPlane
{
std::string id;
std::vector<LinkAtom> atoms;
float esd;
};
// --------------------------------------------------------------------
// struct containing information about a chiral centre
// This information comes from the CCP4 monomer library.
struct LinkChiralCentre
{
std::string id;
LinkAtom atomCentre;
LinkAtom atom[3];
ChiralVolumeSign volumeSign;
};
// --------------------------------------------------------------------
// a class that contains information about a chemical link between compounds.
// This information is derived from the ccp4 monomer library by default.
class Link
{
public:
Link(cif::Datablock& db);
// Factory method.
static const Link& create(const std::string& id);
// accessors
std::string id() const { return mId; }
std::vector<LinkBond> bonds() const { return mBonds; }
std::vector<LinkAngle> angles() const { return mAngles; }
std::vector<LinkChiralCentre> chiralCentres() const { return mChiralCentres; }
std::vector<LinkPlane> planes() const { return mPlanes; }
private:
~Link();
std::string mId;
std::vector<LinkBond> mBonds;
std::vector<LinkAngle> mAngles;
std::vector<LinkChiralCentre> mChiralCentres;
std::vector<LinkPlane> mPlanes;
};
// --------------------------------------------------------------------
// Factory class for Compound and Link objects
extern const std::map<std::string,char> kAAMap, kBaseMap;
......@@ -180,6 +266,9 @@ class CompoundFactory
const Compound* get(std::string id);
const Compound* create(std::string id);
const Link* getLink(std::string id);
const Link* createLink(std::string id);
private:
......
include/cif++/DistanceMap.h
......@@ -8,7 +8,7 @@
#include "cif++/Structure.h"
namespace libcif
namespace mmcif
{
class DistanceMap
......
include/cif++/Edia.h
......@@ -11,7 +11,7 @@
#include "cif++/DistanceMap.h"
#include "cif++/BondMap.h"
namespace libcif
namespace mmcif
{
// --------------------------------------------------------------------
......
include/cif++/MapMaker.h
......@@ -5,7 +5,7 @@
#include "cif++/Structure.h"
namespace libcif
namespace mmcif
{
template<typename FTYPE>
......
include/cif++/PeptideDB.h deleted 100644 → 0
#pragma once
#include <map>
#include <string>
#include <boost/filesystem/path.hpp>
extern const std::map<std::string,char> kAAMap, kBaseMap;
class PeptideDB
{
public:
static PeptideDB& Instance();
void pushDictionary(boost::filesystem::path dict);
void popDictionary();
bool isKnownPeptide(const std::string& res_name) const;
bool isKnownBase(const std::string& res_name) const;
std::string nameForResidue(const std::string& res_name) const;
std::string formulaForResidue(const std::string& res_name) const;
std::string unalias(const std::string& res_name) const;
private:
PeptideDB();
~PeptideDB();
PeptideDB(const PeptideDB&) = delete;
PeptideDB& operator=(const PeptideDB&) = delete;
struct PeptideDBImpl* mImpl;
static PeptideDB* sInstance;
};
include/cif++/Point.h
......@@ -9,7 +9,7 @@
#include "clipper/core/coords.h"
namespace libcif
namespace mmcif
{
typedef boost::math::quaternion<float> quaternion;
......@@ -26,69 +26,71 @@ const long double
// float x, y, z;
// tie(x, y, z) = atom.loc();
struct Point : public std::tuple<float,float,float>
struct Point
{
typedef std::tuple<float,float,float> base_type;
float mX, mY, mZ;
Point() : base_type(0.f, 0.f, 0.f) {}
Point(float x, float y, float z) : base_type(x, y, z) {}
Point(const clipper::Coord_orth& pt): base_type(pt[0], pt[1], pt[2]) {}
Point() : mX(0), mY(0), mZ(0) {}
Point(float x, float y, float z) : mX(x), mY(y), mZ(z) {}
Point(const clipper::Coord_orth& pt): mX(pt[0]), mY(pt[1]), mZ(pt[2]) {}
Point& operator=(const clipper::Coord_orth& rhs)
{
setX(rhs[0]);
setY(rhs[1]);
setZ(rhs[2]);
mX = rhs[0];
mY = rhs[1];
mZ = rhs[2];
return *this;
}
float& getX() { return std::get<0>(*this); }
float getX() const { return std::get<0>(*this); }
void setX(float x) { std::get<0>(*this) = x; }
float& getX() { return mX; }
float getX() const { return mX; }
void setX(float x) { mX = x; }
float& getY() { return std::get<1>(*this); }
float getY() const { return std::get<1>(*this); }
void setY(float y) { std::get<1>(*this) = y; }
float& getY() { return mY; }
float getY() const { return mY; }
void setY(float y) { mY = y; }
float& getZ() { return std::get<2>(*this); }
float getZ() const { return std::get<2>(*this); }
void setZ(float z) { std::get<2>(*this) = z; }
float& getZ() { return mZ; }
float getZ() const { return mZ; }
void setZ(float z) { mZ = z; }
Point& operator+=(const Point& rhs)
{
getX() += rhs.getX();
getY() += rhs.getY();
getZ() += rhs.getZ();
mX += rhs.mX;
mY += rhs.mY;
mZ += rhs.mZ;
return *this;
}
Point& operator-=(const Point& rhs)
{
getX() -= rhs.getX();
getY() -= rhs.getY();
getZ() -= rhs.getZ();
mX -= rhs.mX;
mY -= rhs.mY;
mZ -= rhs.mZ;
return *this;
}
Point& operator*=(float rhs)
{
getX() *= rhs;
getY() *= rhs;
getZ() *= rhs;
mX *= rhs;
mY *= rhs;
mZ *= rhs;
return *this;
}
Point& operator/=(float rhs)
{
getX() *= rhs;
getY() *= rhs;
getZ() *= rhs;
mX *= rhs;
mY *= rhs;
mZ *= rhs;
return *this;
}
float normalize()
{
auto length = getX() * getX() + getY() * getY() + getZ() * getZ();
auto length = mX * mX + mY * mY + mZ * mZ;
if (length > 0)
{
length = std::sqrt(length);
......@@ -99,50 +101,55 @@ struct Point : public std::tuple<float,float,float>
void rotate(const boost::math::quaternion<float>& q)
{
boost::math::quaternion<float> p(0, getX(), getY(), getZ());
boost::math::quaternion<float> p(0, mX, mY, mZ);
p = q * p * boost::math::conj(q);
getX() = p.R_component_2();
getY() = p.R_component_3();
getZ() = p.R_component_4();
mX = p.R_component_2();
mY = p.R_component_3();
mZ = p.R_component_4();
}
operator clipper::Coord_orth() const
{
return clipper::Coord_orth(getX(), getY(), getZ());
return clipper::Coord_orth(mX, mY, mZ);
}
operator std::tuple<float,float,float>() const
{
return std::make_tuple(mX, mY, mZ);
}
};
inline std::ostream& operator<<(std::ostream& os, const Point& pt)
{
os << '(' << pt.getX() << ',' << pt.getY() << ',' << pt.getZ() << ')';
os << '(' << pt.mX << ',' << pt.mY << ',' << pt.mZ << ')';
return os;
}
inline Point operator+(const Point& lhs, const Point& rhs)
{
return Point(lhs.getX() + rhs.getX(), lhs.getY() + rhs.getY(), lhs.getZ() + rhs.getZ());
return Point(lhs.mX + rhs.mX, lhs.mY + rhs.mY, lhs.mZ + rhs.mZ);
}
inline Point operator-(const Point& lhs, const Point& rhs)
{
return Point(lhs.getX() - rhs.getX(), lhs.getY() - rhs.getY(), lhs.getZ() - rhs.getZ());
return Point(lhs.mX - rhs.mX, lhs.mY - rhs.mY, lhs.mZ - rhs.mZ);
}
inline Point operator-(const Point& pt)
{
return Point(-pt.getX(), -pt.getY(), -pt.getZ());
return Point(-pt.mX, -pt.mY, -pt.mZ);
}
inline Point operator*(const Point& pt, float f)
{
return Point(pt.getX() * f, pt.getY() * f, pt.getZ() * f);
return Point(pt.mX * f, pt.mY * f, pt.mZ * f);
}
inline Point operator/(const Point& pt, float f)
{
return Point(pt.getX() / f, pt.getY() / f, pt.getZ() / f);
return Point(pt.mX / f, pt.mY / f, pt.mZ / f);
}
// --------------------------------------------------------------------
......@@ -151,29 +158,29 @@ inline Point operator/(const Point& pt, float f)
inline double DistanceSquared(const Point& a, const Point& b)
{
return
(a.getX() - b.getX()) * (a.getX() - b.getX()) +
(a.getY() - b.getY()) * (a.getY() - b.getY()) +
(a.getZ() - b.getZ()) * (a.getZ() - b.getZ());
(a.mX - b.mX) * (a.mX - b.mX) +
(a.mY - b.mY) * (a.mY - b.mY) +
(a.mZ - b.mZ) * (a.mZ - b.mZ);
}
inline double Distance(const Point& a, const Point& b)
{
return sqrt(
(a.getX() - b.getX()) * (a.getX() - b.getX()) +
(a.getY() - b.getY()) * (a.getY() - b.getY()) +
(a.getZ() - b.getZ()) * (a.getZ() - b.getZ()));
(a.mX - b.mX) * (a.mX - b.mX) +
(a.mY - b.mY) * (a.mY - b.mY) +
(a.mZ - b.mZ) * (a.mZ - b.mZ));
}
inline float DotProduct(const Point& a, const Point& b)
{
return a.getX() * b.getX() + a.getY() * b.getY() + a.getZ() * b.getZ();
return a.mX * b.mX + a.mY * b.mY + a.mZ * b.mZ;
}
inline Point CrossProduct(const Point& a, const Point& b)
{
return Point(a.getY() * b.getZ() - b.getY() * a.getZ(),
a.getZ() * b.getX() - b.getZ() * a.getX(),
a.getX() * b.getY() - b.getX() * a.getY());
return Point(a.mY * b.mZ - b.mY * a.mZ,
a.mZ * b.mX - b.mZ * a.mX,
a.mX * b.mY - b.mX * a.mY);
}
float DihedralAngle(const Point& p1, const Point& p2, const Point& p3, const Point& p4);
......@@ -238,9 +245,9 @@ class SphericalDots
double lat = std::asin((2.0 * i) / P);
double lon = std::fmod(i, kGoldenRatio) * 2 * kPI / kGoldenRatio;
p->setX(sin(lon) * cos(lat));
p->setY(cos(lon) * cos(lat));
p->setZ( sin(lat));
p->mX = sin(lon) * cos(lat);
p->mY = cos(lon) * cos(lat);
p->mZ = sin(lat);
++p;
}
......@@ -254,5 +261,4 @@ class SphericalDots
typedef SphericalDots<50> SphericalDots_50;
}
include/cif++/ResolutionCalculator.h
......@@ -8,7 +8,7 @@
#include <cmath>
namespace libcif
namespace mmcif
{
// --------------------------------------------------------------------
......
include/cif++/Structure.h
......@@ -45,7 +45,7 @@
// class RowSet;
//};
namespace libcif
namespace mmcif
{
class Atom;
......@@ -88,7 +88,7 @@ class Atom
Atom& operator=(const Atom& rhs);
std::string id() const;
const std::string& id() const;
AtomType type() const;
Point location() const;
......@@ -144,6 +144,7 @@ typedef std::vector<Atom> AtomView;
class Residue
{
public:
Residue();
Residue(const Residue& rhs);
Residue& operator=(const Residue& rhs);
......@@ -184,6 +185,7 @@ class Residue
class Monomer : public Residue
{
public:
Monomer();
Monomer(const Monomer& rhs);
Monomer& operator=(const Monomer& rhs);
......
src/AtomShape.cpp
......@@ -8,7 +8,7 @@
using namespace std;
namespace libcif
namespace mmcif
{
// --------------------------------------------------------------------
......@@ -118,7 +118,7 @@ double sineIntegration(double x)
double gx = (sn / sd) / (x * x);
result = libcif::kPI / 2 - fx * cos(x) - gx * sin(x);
result = mmcif::kPI / 2 - fx * cos(x) - gx * sin(x);
}
return result;
......@@ -194,7 +194,7 @@ DensityIntegration::DensityIntegration(float resolutionLow, float resolutionHigh
{
double z, zo, dp;
z = cos(libcif::kPI * (i - 0.25) / (N + 0.5));
z = cos(mmcif::kPI * (i - 0.25) / (N + 0.5));
do
{
......@@ -237,7 +237,7 @@ double DensityIntegration::integrateDensity(double r, int ks, const vector<doubl
if (rt > 1e-10)
{
double t = 4 * libcif::kPI * rt;
double t = 4 * mmcif::kPI * rt;
y = 0;
for (size_t i = 0; i < mST.size(); ++i)
......@@ -252,7 +252,7 @@ double DensityIntegration::integrateDensity(double r, int ks, const vector<doubl
double DensityIntegration::integrateRadius(float perc, float occupancy, double yi, const vector<double>& fst) const
{
double yt = perc * 0.25 * libcif::kPI * occupancy * yi;
double yt = perc * 0.25 * mmcif::kPI * occupancy * yi;
double initialValue = 0.25;
......
src/AtomType.cpp
......@@ -5,7 +5,7 @@
using namespace std;
namespace libcif
namespace mmcif
{
namespace data
......
src/BondMap.cpp
......@@ -9,7 +9,7 @@
using namespace std;
namespace libcif
namespace mmcif
{
// --------------------------------------------------------------------
......@@ -128,7 +128,7 @@ BondMap::BondMap(const Structure& p)
for (auto c: compounds)
{
auto* compound = libcif::Compound::create(c);
auto* compound = mmcif::Compound::create(c);
if (not compound)
{
cerr << "Missing compound information for " << c << endl;
......
src/Cif2PDB.cpp
......@@ -13,9 +13,9 @@
#include <boost/iostreams/concepts.hpp> // output_filter
#include <boost/iostreams/operations.hpp> // put
#include "cif++/PeptideDB.h"
#include "cif++/Cif2PDB.h"
#include "cif++/AtomType.h"
#include "cif++/Compound.h"
using namespace std;
namespace ba = boost::algorithm;
......@@ -2666,7 +2666,7 @@ int WriteHeterogen(ostream& pdbFile, Datablock& db)
cif::tie(entity_id, seqNum, comp_id, chain_id, iCode, modelNr) =
r.get("label_entity_id", "auth_seq_id", "auth_comp_id", "auth_asym_id", "pdbx_PDB_ins_code", "pdbx_PDB_model_num");
if (kAAMap.count(comp_id) or kBaseMap.count(comp_id))
if (mmcif::kAAMap.count(comp_id) or mmcif::kBaseMap.count(comp_id))
continue;
if (chain_id.length() != 1)
......
src/Compound.cpp
......@@ -21,7 +21,7 @@ namespace ba = boost::algorithm;
namespace fs = boost::filesystem;
namespace zx = zeep::xml;
namespace libcif
namespace mmcif
{
// --------------------------------------------------------------------
......@@ -141,10 +141,10 @@ const vector<string>& IsomerDB::operator[](const string& compound) const
struct Node
{
string id;
AtomType symbol;
vector<tuple<size_t,CompoundBondType>> links;
size_t hydrogens = 0;
string id;
AtomType symbol;
vector<tuple<size_t,BondType>> links;
size_t hydrogens = 0;
};
// Check to see if the nodes a[iA] and b[iB] are the start of a similar sub structure
......@@ -179,7 +179,7 @@ bool SubStructuresAreIsomeric(
for (size_t i = 0; result and i < N; ++i)
{
size_t lA, lB;
CompoundBondType typeA, typeB;
BondType typeA, typeB;
tie(lA, typeA) = na.links[i]; assert(lA < a.size());
tie(lB, typeB) = nb.links[ilb[i]]; assert(lB < b.size());
......@@ -372,7 +372,7 @@ Compound::Compound(const fs::path& file, const std::string& id,
for (auto row: compChir)
{
ChiralCentre cc;
CompoundChiralCentre cc;
string volumeSign;
cif::tie(cc.id, cc.atomIDCentre, cc.atomID[0],
......@@ -788,6 +788,125 @@ float Compound::chiralVolume(const string& centreID) const
}
// --------------------------------------------------------------------
Link::Link(cif::Datablock& db)
{
mId = db.getName();
auto& linkBonds = db["chem_link_bond"];
for (auto row: linkBonds)
{
LinkBond b;
string type, aromatic;
cif::tie(b.atom[0].compID, b.atom[0].atomID,
b.atom[1].compID, b.atom[1].atomID, type, b.distance, b.esd) =
row.get("atom_1_comp_id", "atom_id_1", "atom_2_comp_id", "atom_id_2",
"type", "value_dist", "value_dist_esd");
using cif::iequals;
if (iequals(type, "single") or iequals(type, "sing")) b.type = singleBond;
else if (iequals(type, "double") or iequals(type, "doub")) b.type = doubleBond;
else if (iequals(type, "triple") or iequals(type, "trip")) b.type = tripleBond;
else if (iequals(type, "deloc") or iequals(type, "aromat") or iequals(type, "aromatic"))
b.type = delocalizedBond;
else
{
if (VERBOSE)
cerr << "Unimplemented chem_link_bond.type " << type << " in " << mId << endl;
b.type = singleBond;
}
// if (b.atom[0] > b.atom[1])
// swap(b.atom[0], b.atom[1]);
mBonds.push_back(b);
}
// sort(mBonds.begin(), mBonds.end(), LinkBondLess());
auto& linkAngles = db["chem_link_angle"];
for (auto row: linkAngles)
{
LinkAngle a;
cif::tie(a.atom[0].compID, a.atom[0].atomID, a.atom[1].compID, a.atom[1].atomID,
a.atom[2].compID, a.atom[2].atomID, a.angle, a.esd) =
row.get("atom_1_comp_id", "atom_id_1", "atom_2_comp_id", "atom_id_2",
"atom_3_comp_id", "atom_id_3", "value_angle", "value_angle_esd");
mAngles.push_back(a);
}
auto& linkChir = db["chem_link_chir"];
for (auto row: linkChir)
{
LinkChiralCentre cc;
string volumeSign;
cif::tie(cc.id, cc.atomCentre.compID, cc.atomCentre.atomID,
cc.atom[0].compID, cc.atom[0].atomID,
cc.atom[1].compID, cc.atom[1].atomID,
cc.atom[2].compID, cc.atom[2].atomID,
volumeSign) =
row.get("id", "atom_centre_comp_id", "atom_id_centre",
"atom_1_comp_id", "atom_id_1", "atom_2_comp_id", "atom_id_2",
"atom_3_comp_id", "atom_id_3", "volume_sign");
if (volumeSign == "negativ" or volumeSign == "negative")
cc.volumeSign = negativ;
else if (volumeSign == "positiv" or volumeSign == "positive")
cc.volumeSign = positiv;
else if (volumeSign == "both")
cc.volumeSign = both;
else
{
if (VERBOSE)
cerr << "Unimplemented chem_link_chir.volume_sign " << volumeSign << " in " << mId << endl;
continue;
}
mChiralCentres.push_back(cc);
}
auto& linkPlanes = db["chem_link_plane"];
for (auto row: linkPlanes)
{
int compID;
string atomID, planeID;
float esd;
cif::tie(planeID, compID, atomID, esd) = row.get("plane_id", "atom_comp_id", "atom_id", "dist_esd");
auto i = find_if(mPlanes.begin(), mPlanes.end(), [&](auto& p) { return p.id == planeID;});
if (i == mPlanes.end())
{
vector<LinkAtom> atoms{LinkAtom{ compID, atomID }};
mPlanes.emplace_back(LinkPlane{ planeID, move(atoms), esd });
}
else
i->atoms.push_back({ compID, atomID });
}
}
const Link& Link::create(const std::string& id)
{
auto result = CompoundFactory::instance().getLink(id);
if (result == nullptr)
result = CompoundFactory::instance().createLink(id);
if (result == nullptr)
throw runtime_error("Link with id " + id + " not found");
return *result;
}
Link::~Link()
{
}
// --------------------------------------------------------------------
// a factory class to generate compounds
const map<string,char> kAAMap{
......@@ -841,6 +960,9 @@ class CompoundFactoryImpl
Compound* get(string id);
Compound* create(string id);
const Link* getLink(std::string id);
const Link* createLink(std::string id);
CompoundFactoryImpl* pop()
{
......@@ -889,8 +1011,10 @@ class CompoundFactoryImpl
fs::path mPath;
vector<Compound*> mCompounds;
vector<const Link*> mLinks;
/*unordered_*/set<string> mKnownPeptides;
set<string> mKnownBases;
set<string> mMissing;
cif::File mFile;
cif::Category& mChemComp;
CompoundFactoryImpl* mNext;
......@@ -1020,7 +1144,7 @@ Compound* CompoundFactoryImpl::create(std::string id)
ba::to_upper(id);
Compound* result = get(id);
if (result == nullptr)
if (result == nullptr and mMissing.count(id) == 0)
{
boost::upgrade_lock<boost::shared_mutex> lock(mMutex);
......@@ -1042,12 +1166,15 @@ Compound* CompoundFactoryImpl::create(std::string id)
{
auto clibd_mon = fs::path(getenv("CLIBD_MON"));
fs::path resFile = clibd_mon / ba::to_lower_copy(name.substr(0, 1)) / (name + ".cif");
fs::path resFile = clibd_mon / ba::to_lower_copy(name.substr(0, 1)) / (id + ".cif");
if (not fs::exists(resFile) and (name == "COM" or name == "CON" or "PRN")) // seriously...
resFile = clibd_mon / ba::to_lower_copy(name.substr(0, 1)) / (name + '_' + name + ".cif");
if (not fs::exists(resFile) and (id == "COM" or id == "CON" or "PRN")) // seriously...
resFile = clibd_mon / ba::to_lower_copy(id.substr(0, 1)) / (id + '_' + id + ".cif");
mCompounds.push_back(new Compound(resFile, id, name, group));
if (not fs::exists(resFile))
mMissing.insert(id);
else
mCompounds.push_back(new Compound(resFile, id, name, group));
}
else
mCompounds.push_back(new Compound(mPath, id, name, group));
......@@ -1062,6 +1189,53 @@ Compound* CompoundFactoryImpl::create(std::string id)
return result;
}
const Link* CompoundFactoryImpl::getLink(string id)
{
boost::shared_lock<boost::shared_mutex> lock(mMutex);
ba::to_upper(id);
const Link* result = nullptr;
for (auto link: mLinks)
{
if (link->id() == id)
{
result = link;
break;
}
}
if (result == nullptr and mNext != nullptr)
result = mNext->getLink(id);
return result;
}
const Link* CompoundFactoryImpl::createLink(std::string id)
{
ba::to_upper(id);
const Link* result = getLink(id);
if (result == nullptr)
{
boost::upgrade_lock<boost::shared_mutex> lock(mMutex);
auto db = mFile.get("link_" + id);
if (db != nullptr)
{
result = new Link(*db);
mLinks.push_back(result);
}
if (result == nullptr and mNext != nullptr)
result = mNext->createLink(id);
}
return result;
}
// --------------------------------------------------------------------
CompoundFactory* CompoundFactory::sInstance;
......@@ -1119,6 +1293,16 @@ const Compound* CompoundFactory::create(std::string id)
return mImpl->create(id);
}
const Link* CompoundFactory::getLink(std::string id)
{
return mImpl->getLink(id);
}
const Link* CompoundFactory::createLink(std::string id)
{
return mImpl->createLink(id);
}
bool CompoundFactory::isKnownPeptide(const string& resName) const
{
return mImpl->isKnownPeptide(resName);
......
src/DistanceMap.cpp
......@@ -11,7 +11,7 @@
using namespace std;
namespace libcif
namespace mmcif
{
// --------------------------------------------------------------------
......
src/Edia.cpp
......@@ -13,7 +13,7 @@ using namespace std;
extern int VERBOSE;
namespace libcif
namespace mmcif
{
// --------------------------------------------------------------------
......@@ -163,7 +163,7 @@ float CalculateEDIA(const Atom& atom, const clipper::Xmap<float>& xmap,
float radius = kAtomRadius(atom.type(), 0, resolution);
float x, y, z;
tie(x, y, z) = atom.location();
tie(x, y, z) = (tuple<float,float,float>)atom.location();
// calculate min and max orthogonal coordinates first, based on atom position, radius and bfactor
clipper::Coord_orth oMin = { x - radius * 2, y - radius * 2, z - radius * 2 },
......
src/MapMaker.cpp
......@@ -23,7 +23,7 @@ namespace ba = boost::algorithm;
extern int VERBOSE;
namespace libcif
namespace mmcif
{
// --------------------------------------------------------------------
......
src/PDB2Cif.cpp
......@@ -14,7 +14,6 @@
#include "cif++/PDB2Cif.h"
#include "cif++/AtomType.h"
#include "cif++/Compound.h"
#include "cif++/PeptideDB.h"
#include "cif++/PDB2CifRemark3.h"
#include "cif++/CifUtils.h"
......@@ -26,6 +25,7 @@ using cif::Category;
using cif::Row;
using cif::Key;
using cif::iequals;
using mmcif::CompoundFactory;
// --------------------------------------------------------------------
// attempt to come up with better error handling
......@@ -2369,14 +2369,14 @@ void PDBFileParser::ParseRemarks()
int seq = vI(22, 25);
char iCode = vC(26);
auto compound = libcif::Compound::create(res);
auto compound = mmcif::Compound::create(res);
if (compound == nullptr)
continue;
vector<string> atoms;
for (auto atom: compound->atoms())
{
if (atom.typeSymbol != libcif::H)
if (atom.typeSymbol != mmcif::H)
atoms.push_back(atom.id);
}
......@@ -3281,7 +3281,7 @@ void PDBFileParser::ConstructEntities()
PDBChain::AtomRes ar{ resName, resSeq, iCode };
if ((chain.mResiduesSeen.empty() or chain.mResiduesSeen.back() != ar) and
(PeptideDB::Instance().isKnownPeptide(resName) or PeptideDB::Instance().isKnownBase(resName)))
(CompoundFactory::instance().isKnownPeptide(resName) or CompoundFactory::instance().isKnownBase(resName)))
{
chain.mResiduesSeen.push_back(ar);
}
......@@ -3361,10 +3361,10 @@ void PDBFileParser::ConstructEntities()
{
string resName = chain.mResiduesSeen[ix].mMonId;
if (kAAMap.count(resName) or
kBaseMap.count(resName) or
PeptideDB::Instance().isKnownPeptide(resName) or
PeptideDB::Instance().isKnownBase(resName))
if (mmcif::kAAMap.count(resName) or
mmcif::kBaseMap.count(resName) or
CompoundFactory::instance().isKnownPeptide(resName) or
CompoundFactory::instance().isKnownBase(resName))
{
chain.mTerIndex = ix + 1;
}
......@@ -3795,14 +3795,14 @@ void PDBFileParser::ConstructEntities()
if (mMod2parent.count(res.mMonId))
stdRes = mMod2parent.at(res.mMonId);
if (kAAMap.count(res.mMonId))
if (mmcif::kAAMap.count(res.mMonId))
{
letter = kAAMap.at(res.mMonId);
letter = mmcif::kAAMap.at(res.mMonId);
mightBeDNA = false;
}
else if (kBaseMap.count(res.mMonId))
else if (mmcif::kBaseMap.count(res.mMonId))
{
letter = kBaseMap.at(res.mMonId);
letter = mmcif::kBaseMap.at(res.mMonId);
mightBePolyPeptide = false;
}
else
......@@ -3811,7 +3811,7 @@ void PDBFileParser::ConstructEntities()
letter = '(' + res.mMonId + ')';
// sja...
auto compound = libcif::Compound::create(stdRes.empty() ? res.mMonId : stdRes);
auto compound = mmcif::Compound::create(stdRes.empty() ? res.mMonId : stdRes);
if (compound != nullptr and
not iequals(compound->type(), "L-peptide linking") and
not iequals(compound->type(), "RNA linking"))
......@@ -3831,10 +3831,10 @@ void PDBFileParser::ConstructEntities()
if (letter.length() > 1)
{
if (not stdRes.empty() and kAAMap.count(stdRes))
letter = kAAMap.at(stdRes);
else if (kBaseMap.count(res.mMonId))
letter = kBaseMap.at(res.mMonId);
if (not stdRes.empty() and mmcif::kAAMap.count(stdRes))
letter = mmcif::kAAMap.at(stdRes);
else if (mmcif::kBaseMap.count(res.mMonId))
letter = mmcif::kBaseMap.at(res.mMonId);
else
letter = 'X';
}
......@@ -3961,7 +3961,11 @@ void PDBFileParser::ConstructEntities()
else
{
if (mHetnams[hetID].empty())
mHetnams[hetID] = PeptideDB::Instance().nameForResidue(hetID);
{
auto compound = mmcif::Compound::create(hetID);
if (compound != nullptr)
mHetnams[hetID] = compound->name();
}
getCategory("entity")->emplace({
{ "id", entityId },
......@@ -4095,7 +4099,7 @@ void PDBFileParser::ConstructEntities()
for (auto cc: mChemComp)
{
auto compound = libcif::Compound::create(
auto compound = mmcif::Compound::create(
mMod2parent.count(cc) ? mMod2parent[cc] : cc
);
......@@ -4447,10 +4451,10 @@ static bool IsMetal(const string& resName, const string& atomID)
try
{
auto compound = libcif::Compound::create(resName);
auto compound = mmcif::Compound::create(resName);
if (compound != nullptr)
{
auto at = libcif::AtomTypeTraits(compound->getAtomById(atomID).typeSymbol);
auto at = mmcif::AtomTypeTraits(compound->getAtomById(atomID).typeSymbol);
result = at.isMetal();
}
}
......
src/PDB2CifRemark3.cpp
......@@ -11,7 +11,6 @@
#include "cif++/AtomType.h"
#include "cif++/Compound.h"
#include "cif++/PDB2CifRemark3.h"
#include "cif++/PeptideDB.h"
using namespace std;
namespace ba = boost::algorithm;
......
src/PeptideDB.cpp deleted 100644 → 0
#include "cif++/Config.h"
#include <set>
#include <map>
#include <unordered_set>
#include <boost/filesystem/operations.hpp>
#include <boost/filesystem/fstream.hpp>
#include <boost/algorithm/string.hpp>
#include "cif++/Cif++.h"
#include "cif++/PeptideDB.h"
using namespace std;
namespace fs = boost::filesystem;
namespace ba = boost::algorithm;
const map<string,char> kAAMap{
{ "ALA", 'A' },
{ "ARG", 'R' },
{ "ASN", 'N' },
{ "ASP", 'D' },
{ "CYS", 'C' },
{ "GLN", 'Q' },
{ "GLU", 'E' },
{ "GLY", 'G' },
{ "HIS", 'H' },
{ "ILE", 'I' },
{ "LEU", 'L' },
{ "LYS", 'K' },
{ "MET", 'M' },
{ "PHE", 'F' },
{ "PRO", 'P' },
{ "SER", 'S' },
{ "THR", 'T' },
{ "TRP", 'W' },
{ "TYR", 'Y' },
{ "VAL", 'V' },
{ "GLX", 'Z' },
{ "ASX", 'B' }
};
const map<string,char> kBaseMap{
{ "A", 'A' },
{ "C", 'C' },
{ "G", 'G' },
{ "T", 'T' },
{ "U", 'U' },
{ "DA", 'A' },
{ "DC", 'C' },
{ "DG", 'G' },
{ "DT", 'T' }
};
// --------------------------------------------------------------------
class PeptideDBImpl
{
public:
PeptideDBImpl(istream& data, PeptideDBImpl* next);
~PeptideDBImpl()
{
delete mNext;
}
PeptideDBImpl* pop()
{
auto result = mNext;
mNext = nullptr;
delete this;
return result;
}
string nameFor(const string& resName) const
{
string result;
for (auto& chemComp: mChemComp)
{
if (ba::iequals(chemComp["three_letter_code"].as<string>(), resName) == false)
continue;
result = chemComp["name"].as<string>();
ba::trim(result);
break;
}
if (result.empty() and mNext)
result = mNext->nameFor(resName);
return result;
}
string formulaFor(string resName) const;
string unalias(const string& resName) const
{
string result = resName;
auto& e = const_cast<cif::File&>(mFile)["comp_synonym_list"];
for (auto& synonym: e["chem_comp_synonyms"])
{
if (ba::iequals(synonym["comp_alternative_id"].as<string>(), resName) == false)
continue;
result = synonym["comp_id"].as<string>();
ba::trim(result);
break;
}
if (result.empty() and mNext)
result = mNext->unalias(resName);
return result;
}
bool isKnownPeptide(const string& resName)
{
return mKnownPeptides.count(resName) or
(mNext != nullptr and mNext->isKnownPeptide(resName));
}
bool isKnownBase(const string& resName)
{
return mKnownBases.count(resName) or
(mNext != nullptr and mNext->isKnownBase(resName));
}
private:
/*unordered_*/set<string> mKnownPeptides;
set<string> mKnownBases;
cif::File mFile;
cif::Category& mChemComp;
PeptideDBImpl* mNext;
};
PeptideDBImpl::PeptideDBImpl(istream& data, PeptideDBImpl* next)
: mFile(data), mChemComp(mFile.firstDatablock()["chem_comp"]), mNext(next)
{
const std::regex peptideRx("(?:[lmp]-)?peptide", std::regex::icase);
for (auto& chemComp: mChemComp)
{
string group, threeLetterCode;
cif::tie(group, threeLetterCode) = chemComp.get("group", "three_letter_code");
if (std::regex_match(group, peptideRx))
// if (ba::iequals(group, "peptide") or ba::iequals(group, "M-peptide") or ba::iequals(group, "P-peptide"))
mKnownPeptides.insert(threeLetterCode);
else if (ba::iequals(group, "DNA") or ba::iequals(group, "RNA"))
mKnownBases.insert(threeLetterCode);
}
}
string PeptideDBImpl::formulaFor(string res) const
{
string result;
ba::to_upper(res);
for (auto& db: mFile)
{
if (db.getName() != "comp_" + res)
continue;
auto& cat = db["chem_comp_atom"];
map<string,uint32> atoms;
for (auto r: cat)
atoms[r["type_symbol"].as<string>()] += 1;
for (auto a: atoms)
{
if (not result.empty())
result += ' ';
result += a.first;
if (a.second > 1)
result += to_string(a.second);
}
}
if (result.empty())
{
if (mNext != nullptr)
result = mNext->formulaFor(res);
else
{
const char* clibdMon = getenv("CLIBD_MON");
if (clibdMon == nullptr)
throw runtime_error("Cannot locate peptide list, please souce the CCP4 environment");
fs::path resFile = fs::path(clibdMon) / ba::to_lower_copy(res.substr(0, 1)) / (res + ".cif");
if (fs::exists(resFile))
{
fs::ifstream file(resFile);
if (file.is_open())
{
try
{
cif::File cf(file);
auto& cat = cf["comp_" + res]["chem_comp_atom"];
map<string,uint32> atoms;
for (auto r: cat)
atoms[r["type_symbol"].as<string>()] += 1;
for (auto a: atoms)
{
if (not result.empty())
result += ' ';
result += a.first;
if (a.second > 1)
result += to_string(a.second);
}
}
catch (exception& ex)
{
if (VERBOSE)
cerr << ex.what();
result.clear();
}
}
}
}
}
return result;
}
// --------------------------------------------------------------------
PeptideDB* PeptideDB::sInstance;
PeptideDB& PeptideDB::Instance()
{
if (sInstance == nullptr)
sInstance = new PeptideDB();
return *sInstance;
}
PeptideDB::PeptideDB()
: mImpl(nullptr)
{
const char* clibdMon = getenv("CLIBD_MON");
if (clibdMon == nullptr)
throw runtime_error("Cannot locate peptide list, please souce the CCP4 environment");
fs::path db = fs::path(clibdMon) / "list" / "mon_lib_list.cif";
pushDictionary(db);
sInstance = this;
}
void PeptideDB::pushDictionary(boost::filesystem::path dict)
{
if (not fs::exists(dict))
throw runtime_error("file not found: " + dict.string());
fs::ifstream file(dict);
if (not file.is_open())
throw runtime_error("Could not open peptide list " + dict.string());
mImpl = new PeptideDBImpl(file, mImpl);
}
void PeptideDB::popDictionary()
{
if (mImpl != nullptr)
mImpl = mImpl->pop();
}
PeptideDB::~PeptideDB()
{
delete mImpl;
}
bool PeptideDB::isKnownPeptide(const string& resName) const
{
return mImpl->isKnownPeptide(resName);
}
bool PeptideDB::isKnownBase(const string& resName) const
{
return mImpl->isKnownBase(resName);
}
string PeptideDB::nameForResidue(const string& resName) const
{
return mImpl->nameFor(resName);
}
string PeptideDB::formulaForResidue(const string& resName) const
{
return mImpl->formulaFor(resName);
}
string PeptideDB::unalias(const string& resName) const
{
return mImpl->unalias(resName);
}
src/Point.cpp
......@@ -6,7 +6,7 @@
using namespace std;
namespace libcif
namespace mmcif
{
// --------------------------------------------------------------------
......@@ -101,20 +101,20 @@ Point CenterPoints(vector<Point>& Points)
for (Point& pt : Points)
{
t.getX() += pt.getX();
t.getY() += pt.getY();
t.getZ() += pt.getZ();
t.mX += pt.mX;
t.mY += pt.mY;
t.mZ += pt.mZ;
}
t.getX() /= Points.size();
t.getY() /= Points.size();
t.getZ() /= Points.size();
t.mX /= Points.size();
t.mY /= Points.size();
t.mZ /= Points.size();
for (Point& pt : Points)
{
pt.getX() -= t.getX();
pt.getY() -= t.getY();
pt.getZ() -= t.getZ();
pt.mX -= t.mX;
pt.mY -= t.mY;
pt.mZ -= t.mZ;
}
return t;
......@@ -139,9 +139,9 @@ double RMSd(const vector<Point>& a, const vector<Point>& b)
{
valarray<double> d(3);
d[0] = b[i].getX() - a[i].getX();
d[1] = b[i].getY() - a[i].getY();
d[2] = b[i].getZ() - a[i].getZ();
d[0] = b[i].mX - a[i].mX;
d[1] = b[i].mY - a[i].mY;
d[2] = b[i].mZ - a[i].mZ;
d *= d;
......@@ -200,9 +200,9 @@ double LargestDepressedQuarticSolution(double a, double b, double c)
// const Point& a = pa[i];
// const Point& b = pb[i];
//
// M(0, 0) += a.getX() * b.getX(); M(0, 1) += a.getX() * b.getY(); M(0, 2) += a.getX() * b.getZ();
// M(1, 0) += a.getY() * b.getX(); M(1, 1) += a.getY() * b.getY(); M(1, 2) += a.getY() * b.getZ();
// M(2, 0) += a.getZ() * b.getX(); M(2, 1) += a.getZ() * b.getY(); M(2, 2) += a.getZ() * b.getZ();
// M(0, 0) += a.mX * b.mX; M(0, 1) += a.mX * b.mY; M(0, 2) += a.mX * b.mZ;
// M(1, 0) += a.mY * b.mX; M(1, 1) += a.mY * b.mY; M(1, 2) += a.mY * b.mZ;
// M(2, 0) += a.mZ * b.mX; M(2, 1) += a.mZ * b.mY; M(2, 2) += a.mZ * b.mZ;
// }
//
// // Now calculate N, a symmetric 4x4 matrix
......
src/ResolutionCalculator.cpp
......@@ -6,7 +6,7 @@
using namespace std;
namespace libcif
namespace mmcif
{
ResolutionCalculator::ResolutionCalculator(const clipper::Cell& cell)
......
src/Structure.cpp
......@@ -21,7 +21,7 @@ namespace io = boost::iostreams;
extern int VERBOSE;
namespace libcif
namespace mmcif
{
// --------------------------------------------------------------------
......@@ -173,12 +173,11 @@ struct AtomImpl
cif::tie(x, y, z) = mRow.get("Cartn_x", "Cartn_y", "Cartn_z");
mLocation = Point(x, y, z);
string compId;
cif::tie(compId) = mRow.get("label_comp_id");
try
{
comp();
}
catch (...) {}
mCompound = Compound::create(compId);
}
clipper::Atom toClipper() const
......@@ -362,7 +361,7 @@ float Atom::property<float>(const string& name) const
return stof(mImpl->property(name));
}
string Atom::id() const
const string& Atom::id() const
{
return mImpl->mId;
}
......@@ -486,6 +485,11 @@ float Atom::radius() const
// --------------------------------------------------------------------
// residue
Residue::Residue()
: mStructure(nullptr), mSeqID(0)
{
}
Residue::Residue(const Residue& rhs)
: mStructure(rhs.mStructure)
, mCompoundID(rhs.mCompoundID), mAsymID(rhs.mAsymID), mAltID(rhs.mAltID), mSeqID(rhs.mSeqID)
......@@ -602,6 +606,11 @@ bool Residue::isEntity() const
// --------------------------------------------------------------------
// monomer
Monomer::Monomer()
: mPolymer(nullptr), mIndex(0)
{
}
Monomer::Monomer(const Monomer& rhs)
: Residue(rhs), mPolymer(rhs.mPolymer), mIndex(rhs.mIndex)
{
......
src/TlsParser.cpp
......@@ -44,7 +44,7 @@ using namespace std;
namespace po = boost::program_options;
namespace ba = boost::algorithm;
namespace fs = boost::filesystem;
namespace c = libcif;
namespace c = mmcif;
namespace cif
{
......
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