Skip to content

L
libcifpp
Commit 9c051207 by maarten
Options

vergeten! 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@311 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent 91abdc56
Showing with 206 additions and 103 deletions
include/cif++/BondMap.h
...@@ -19,23 +19,46 @@ class BondMap ...@@ -19,23 +19,46 @@ class BondMap
bool operator()(const Atom& a, const Atom& b) const bool operator()(const Atom& a, const Atom& b) const
{ {
return isBonded(index.at(a.id()), index.at(b.id())); return get(index.at(a.id()), index.at(b.id())) == 1;
} }
bool operator()(const std::string& id_a, const std::string& id_b) const bool operator()(const std::string& id_a, const std::string& id_b) const
{ {
return isBonded(index.at(id_a), index.at(id_b)); return get(index.at(id_a), index.at(id_b)) == 1;
} }
// bool is1_4(const Atom& a, const Atom& b) const
// {
// return get(index.at(a.id()), index.at(b.id())) == 3;
// }
//
// bool is1_4(const std::string& id_a, const std::string& id_b) const
// {
// return get(index.at(id_a), index.at(id_b)) == 3;
// }
bool is1_4(const Atom& a, const Atom& b) const; bool is1_4(const Atom& a, const Atom& b) const;
bool is1_4(const std::string& id_a, const std::string& id_b) const; bool is1_4(const std::string& id_a, const std::string& id_b) const;
private: private:
bool isBonded(size_t ai, size_t bi) const; uint8 get(size_t ia, size_t ib) const
{
uint8 result = 0;
if (ia != ib)
{
if (ib < ia)
std::swap(ia, ib);
size_t ix = ib + ia * dim - ia * (ia + 1) / 2;
assert(ix < bond.size());
result = bond[ix];
}
return result;
}
size_t dim; size_t dim;
std::vector<bool> bond; std::vector<uint8> bond;
std::unordered_map<std::string,size_t> index; std::unordered_map<std::string,size_t> index;
}; };
......
include/cif++/DistanceMap.h
...@@ -19,12 +19,12 @@ class DistanceMap ...@@ -19,12 +19,12 @@ class DistanceMap
DistanceMap(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell); DistanceMap(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell);
// simplified version for subsets of atoms (used in refining e.g.) // simplified version for subsets of atoms (used in refining e.g.)
DistanceMap(const Structure& p, const std::vector<Atom>& atoms); // DistanceMap(const Structure& p, const std::vector<Atom>& atoms);
DistanceMap(const DistanceMap&) = delete; DistanceMap(const DistanceMap&) = delete;
DistanceMap& operator=(const DistanceMap&) = delete; DistanceMap& operator=(const DistanceMap&) = delete;
float operator()(const Atom& a, const Atom& b) const; // float operator()(const Atom& a, const Atom& b) const;
std::vector<Atom> near(const Atom& a, float maxDistance = 3.5f) const; std::vector<Atom> near(const Atom& a, float maxDistance = 3.5f) const;
std::vector<Atom> near(const Point& p, float maxDistance = 3.5f) const; std::vector<Atom> near(const Point& p, float maxDistance = 3.5f) const;
...@@ -45,7 +45,9 @@ class DistanceMap ...@@ -45,7 +45,9 @@ class DistanceMap
} }
}; };
std::unordered_map<std::tuple<size_t,size_t>,float,key_hash> dist; std::unordered_map<std::tuple<size_t,size_t>,std::tuple<float,size_t>,key_hash> dist;
Point mD; // needed to move atoms to center
std::vector<clipper::RTop_orth> mRtOrth;
}; };
} }
include/cif++/Structure.h
...@@ -94,7 +94,9 @@ class Atom ...@@ -94,7 +94,9 @@ class Atom
Point location() const; Point location() const;
void location(Point p); void location(Point p);
Atom symmetryCopy(const Point& d, const clipper::RTop_orth& rt);
const Compound& comp() const; const Compound& comp() const;
bool isWater() const; bool isWater() const;
int charge() const; int charge() const;
...@@ -123,6 +125,7 @@ class Atom ...@@ -123,6 +125,7 @@ class Atom
std::string pdbxAuthInsCode() const; std::string pdbxAuthInsCode() const;
std::string authAltId() const; std::string authAltId() const;
std::string labelID() const;// label_comp_id + '_' + label_asym_id + '_' + label_seq_id
std::string pdbID() const; // auth_comp_id + '_' + auth_asym_id + '_' + auth_seq_id + pdbx_PDB_ins_code std::string pdbID() const; // auth_comp_id + '_' + auth_asym_id + '_' + auth_seq_id + pdbx_PDB_ins_code
bool operator==(const Atom& rhs) const; bool operator==(const Atom& rhs) const;
......
src/BondMap.cpp
...@@ -20,7 +20,7 @@ BondMap::BondMap(const Structure& p) ...@@ -20,7 +20,7 @@ BondMap::BondMap(const Structure& p)
auto atoms = p.atoms(); auto atoms = p.atoms();
dim = atoms.size(); dim = atoms.size();
bond = vector<bool>(dim * (dim - 1), false); bond = vector<uint8>(dim * (dim - 1), 0);
for (auto& atom: atoms) for (auto& atom: atoms)
{ {
...@@ -39,7 +39,7 @@ BondMap::BondMap(const Structure& p) ...@@ -39,7 +39,7 @@ BondMap::BondMap(const Structure& p)
size_t ix = ixb + ixa * dim - ixa * (ixa + 1) / 2; size_t ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
assert(ix < bond.size()); assert(ix < bond.size());
bond[ix] = true; bond[ix] = 1;
}; };
...@@ -57,7 +57,8 @@ BondMap::BondMap(const Structure& p) ...@@ -57,7 +57,8 @@ BondMap::BondMap(const Structure& p)
if (compounds.count(c)) if (compounds.count(c))
continue; continue;
cerr << "Warning: mon_id " << c << " is missing in the chem_comp category" << endl; if (VERBOSE)
cerr << "Warning: mon_id " << c << " is missing in the chem_comp category" << endl;
compounds.insert(c); compounds.insert(c);
} }
...@@ -124,8 +125,6 @@ BondMap::BondMap(const Structure& p) ...@@ -124,8 +125,6 @@ BondMap::BondMap(const Structure& p)
// then link all atoms in the compounds // then link all atoms in the compounds
cif::Progress progress(compounds.size(), "Creating bond map");
for (auto c: compounds) for (auto c: compounds)
{ {
auto* compound = mmcif::Compound::create(c); auto* compound = mmcif::Compound::create(c);
...@@ -191,30 +190,59 @@ BondMap::BondMap(const Structure& p) ...@@ -191,30 +190,59 @@ BondMap::BondMap(const Structure& p)
size_t ix = ixb + ixa * dim - ixa * (ixa + 1) / 2; size_t ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
assert(ix < bond.size()); assert(ix < bond.size());
bond[ix] = true; bond[ix] = 1;
} }
} }
} }
} }
progress.consumed(1);
} }
}
bool BondMap::isBonded(size_t ixa, size_t ixb) const
{
if (ixa == ixb)
return false;
if (ixa > ixb) // The next step is to fill bond with the next steps to other atoms
swap(ixa, ixb); // First for two steps and then for three
size_t ix = ixb + ixa * dim - ixa * (ixa + 1) / 2; for (size_t steps: { 2, 3 })
{
assert(ix < bond.size()); for (size_t i = 0; i + 1 < dim; ++i)
return bond[ix]; {
for (size_t j = i + 1; j < dim; ++j)
{
size_t ix = j + i * dim - i * (i + 1) / 2;
if (bond[ix])
continue;
for (size_t k = 0; k < dim; ++k)
{
if (k == i or k == j)
continue;
size_t ni = get(k, i);
size_t nj = get(k, j);
if (ni > 0 and nj > 0 and ni + nj == steps)
{
bond[ix] = steps;
break;
}
}
}
}
}
} }
//bool BondMap::isBonded(size_t ixa, size_t ixb) const
//{
// if (ixa == ixb)
// return false;
//
// if (ixa > ixb)
// swap(ixa, ixb);
//
// size_t ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
//
// assert(ix < bond.size());
// return bond[ix];
//}
bool BondMap::is1_4(const Atom& a, const Atom& b) const bool BondMap::is1_4(const Atom& a, const Atom& b) const
{ {
size_t ixa = index.at(a.id()); size_t ixa = index.at(a.id());
...@@ -227,19 +255,22 @@ bool BondMap::is1_4(const Atom& a, const Atom& b) const ...@@ -227,19 +255,22 @@ bool BondMap::is1_4(const Atom& a, const Atom& b) const
for (size_t ia = 0; result == false and ia + 1 < dim; ++ia) for (size_t ia = 0; result == false and ia + 1 < dim; ++ia)
{ {
if (ia == ixa or ia == ixb or not isBonded(ixa, ia)) if (ia == ixa or ia == ixb or get(ixa, ia) != 1)
continue; continue;
for (size_t ib = ia + 1; result == false and ib < dim; ++ib) for (size_t ib = ia + 1; result == false and ib < dim; ++ib)
{ {
if (ib == ixa or ib == ixb or not isBonded(ib, ixb)) if (ib == ixa or ib == ixb or get(ib, ixb) != 1)
continue; continue;
size_t ix = ib + ia * dim - ia * (ia + 1) / 2; size_t ix = ib + ia * dim - ia * (ia + 1) / 2;
result = bond[ix]; result = bond[ix] == 1;
} }
} }
if (result != (get(ixa, ixb) == 3))
cerr << "Verschil in 1-4 binding voor " << a.labelID() << " en " << b.labelID() << " (c = " << (int)get(ixa, ixb) << ")" << endl;
return result; return result;
} }
......
src/CifUtils.cpp
...@@ -442,7 +442,7 @@ void ProgressImpl::Run() ...@@ -442,7 +442,7 @@ void ProgressImpl::Run()
PrintProgress(); PrintProgress();
printedAny = true; printedAny = true;
boost::this_thread::sleep(boost::posix_time::seconds(0.5)); boost::this_thread::sleep(boost::posix_time::milliseconds(500));
} }
} }
catch (...) {} catch (...) {}
......
src/DistanceMap.cpp
...@@ -125,26 +125,25 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -125,26 +125,25 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
return d; return d;
}; };
float dx = calculateD(mx, cell.a()); mD.mX = calculateD(mx, cell.a());
float dy = calculateD(my, cell.b()); mD.mY = calculateD(my, cell.b());
float dz = calculateD(mz, cell.c()); mD.mZ = calculateD(mz, cell.c());
if (VERBOSE) if (VERBOSE)
cerr << "Calculated d: " << dx << ", " << dy << " and " << dz << endl; cerr << "Calculated d: " << mD.mX << ", " << mD.mY << " and " << mD.mZ << endl;
if (dx != 0 or dy != 0 or dz != 0) if (mD.mX != 0 or mD.mY != 0 or mD.mZ != 0)
{ {
if (VERBOSE) if (VERBOSE)
cerr << "moving coorinates by " << dx << ", " << dy << " and " << dz << endl; cerr << "moving coorinates by " << mD.mX << ", " << mD.mY << " and " << mD.mZ << endl;
for_each(locations.begin(), locations.end(), [&](auto& p) { p[0] += dx; p[1] += dy; p[2] += dz; }); for_each(locations.begin(), locations.end(), [&](auto& p) { p += mD; });
} }
pMin -= kMaxDistance; // extend bounding box pMin -= kMaxDistance; // extend bounding box
pMax += kMaxDistance; pMax += kMaxDistance;
vector<clipper::RTop_orth> mRtOrth = AlternativeSites(spacegroup, cell);
rtOrth = AlternativeSites(spacegroup, cell);
cif::Progress progress(locations.size() - 1, "Creating distance map"); cif::Progress progress(locations.size() - 1, "Creating distance map");
...@@ -170,15 +169,21 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -170,15 +169,21 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
// find nearest location based on spacegroup/cell // find nearest location based on spacegroup/cell
double minR2 = numeric_limits<double>::max(); double minR2 = numeric_limits<double>::max();
for (auto rt: rtOrth) size_t kb = 0;
for (size_t k = 0; k < mRtOrth.size(); ++k)
{ {
auto& rt = mRtOrth[k];
auto pj = locations[j]; auto pj = locations[j];
pj = pj.transform(rt); pj = pj.transform(rt);
double r2 = (pi - pj).lengthsq(); double r2 = (pi - pj).lengthsq();
if (minR2 > r2) if (minR2 > r2)
{
minR2 = r2; minR2 = r2;
kb = k;
}
#if defined(DEBUG_VOOR_BART) #if defined(DEBUG_VOOR_BART)
if (r2 < 3.5 * 3.5 and not rt.equals(clipper::RTop<>::identity(), 0.1)) if (r2 < 3.5 * 3.5 and not rt.equals(clipper::RTop<>::identity(), 0.1))
...@@ -191,10 +196,10 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -191,10 +196,10 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
if (minR2 < kMaxDistanceSQ) if (minR2 < kMaxDistanceSQ)
{ {
float d = sqrt(minR2); float d = sqrt(minR2);
auto k = make_tuple(i, j); auto key = make_tuple(i, j);
lock_guard<mutex> lock(m); lock_guard<mutex> lock(m);
dist[k] = d; dist[key] = make_tuple(d, kb);
} }
} }
...@@ -205,64 +210,63 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -205,64 +210,63 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
t.join_all(); t.join_all();
} }
DistanceMap::DistanceMap(const Structure& p, const vector<Atom>& atoms) //DistanceMap::DistanceMap(const Structure& p, const vector<Atom>& atoms)
: structure(p), dim(0) // : structure(p), dim(0)
{ //{
dim = atoms.size(); // dim = atoms.size();
//
for (auto& atom: atoms) // for (auto& atom: atoms)
{ // {
size_t ix = index.size(); // size_t ix = index.size();
index[atom.id()] = ix; // index[atom.id()] = ix;
}; // };
//
for (size_t i = 0; i < dim; ++i) // for (size_t i = 0; i < dim; ++i)
{ // {
for (size_t j = i + 1; j < dim; ++j) // for (size_t j = i + 1; j < dim; ++j)
{ // {
dist[make_tuple(i, j)] = Distance(atoms[i].location(), atoms[j].location()); // dist[make_tuple(i, j)] = Distance(atoms[i].location(), atoms[j].location());
} // }
} // }
} //}
float DistanceMap::operator()(const Atom& a, const Atom& b) const
{
size_t ixa, ixb;
try
{
ixa = index.at(a.id());
}
catch (const out_of_range& ex)
{
throw runtime_error("atom " + a.id() + " not found in distance map");
}
try
{
ixb = index.at(b.id());
}
catch (const out_of_range& ex)
{
throw runtime_error("atom " + b.id() + " not found in distance map");
}
if (ixb < ixa)
swap(ixa, ixb);
tuple<size_t,size_t> k{ ixa, ixb };
auto ii = dist.find(k);
float result = 100; //float DistanceMap::operator()(const Atom& a, const Atom& b) const
//{
if (ii != dist.end()) // size_t ixa, ixb;
result = ii->second; //
// try
return result; // {
} // ixa = index.at(a.id());
// }
// catch (const out_of_range& ex)
// {
// throw runtime_error("atom " + a.id() + " not found in distance map");
// }
//
// try
// {
// ixb = index.at(b.id());
// }
// catch (const out_of_range& ex)
// {
// throw runtime_error("atom " + b.id() + " not found in distance map");
// }
//
// if (ixb < ixa)
// swap(ixa, ixb);
//
// tuple<size_t,size_t> k{ ixa, ixb };
//
// auto ii = dist.find(k);
//
// float result = 100;
//
// if (ii != dist.end())
// result = ii->second;
//
// return result;
//}
#warning("this method should return symmetry reoriented atoms...")
vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const
{ {
vector<Atom> result; vector<Atom> result;
...@@ -289,8 +293,19 @@ vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const ...@@ -289,8 +293,19 @@ vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const
dist.find(make_tuple(ixa, ixb)) : dist.find(make_tuple(ixa, ixb)) :
dist.find(make_tuple(ixb, ixa)); dist.find(make_tuple(ixb, ixa));
if (ii != dist.end() and ii->second <= maxDistance) if (ii == dist.end())
result.push_back(structure.getAtomById(i.first)); continue;
float distance;
size_t rti;
tie(distance, rti) = ii->second;
if (distance > maxDistance)
continue;
Atom a = structure.getAtomById(i.first);
result.push_back(a.symmetryCopy(mD, mRtOrth.at(rti)));
} }
return result; return result;
......
src/Structure.cpp
...@@ -161,7 +161,20 @@ struct AtomImpl ...@@ -161,7 +161,20 @@ struct AtomImpl
{ {
prefetch(); prefetch();
} }
AtomImpl(const AtomImpl& impl, const Point& d, const clipper::RTop_orth& rt)
: mFile(impl.mFile), mId(impl.mId), mType(impl.mType), mAtomID(impl.mAtomID)
, mCompID(impl.mCompID), mAsymID(impl.mAsymID), mSeqID(impl.mSeqID)
, mAltID(impl.mAltID), mLocation(impl.mLocation), mRefcount(impl.mRefcount)
, mRow(impl.mRow), mCompound(impl.mCompound), mRadius(impl.mRadius)
, mCachedProperties(impl.mCachedProperties)
, mSymmetryCopy(true)
{
mLocation += d;
mLocation = ((clipper::Coord_orth)mLocation).transform(rt);
mLocation -= d;
}
void prefetch() void prefetch()
{ {
// Prefetch some data // Prefetch some data
...@@ -259,6 +272,10 @@ struct AtomImpl ...@@ -259,6 +272,10 @@ struct AtomImpl
void moveTo(const Point& p) void moveTo(const Point& p)
{ {
assert(not mSymmetryCopy);
if (mSymmetryCopy)
throw runtime_error("Moving symmetry copy");
boost::format kPosFmt("%.3f"); boost::format kPosFmt("%.3f");
mRow["Cartn_x"] = (kPosFmt % p.getX()).str(); mRow["Cartn_x"] = (kPosFmt % p.getX()).str();
...@@ -330,6 +347,8 @@ struct AtomImpl ...@@ -330,6 +347,8 @@ struct AtomImpl
const Compound* mCompound; const Compound* mCompound;
float mRadius = nan("4"); float mRadius = nan("4");
map<string,string> mCachedProperties; map<string,string> mCachedProperties;
bool mSymmetryCopy = false;
}; };
//Atom::Atom(const File& f, const string& id) //Atom::Atom(const File& f, const string& id)
...@@ -477,6 +496,11 @@ int Atom::authSeqId() const ...@@ -477,6 +496,11 @@ int Atom::authSeqId() const
return property<int>("auth_seq_id"); return property<int>("auth_seq_id");
} }
string Atom::labelID() const
{
return property<string>("label_comp_id") + '_' + mImpl->mAsymID + '_' + to_string(mImpl->mSeqID) + ':' + mImpl->mAtomID;
}
string Atom::pdbID() const string Atom::pdbID() const
{ {
return return
...@@ -496,6 +520,11 @@ void Atom::location(Point p) ...@@ -496,6 +520,11 @@ void Atom::location(Point p)
mImpl->moveTo(p); mImpl->moveTo(p);
} }
Atom Atom::symmetryCopy(const Point& d, const clipper::RTop_orth& rt)
{
return Atom(new AtomImpl(*mImpl, d, rt));
}
const Compound& Atom::comp() const const Compound& Atom::comp() const
{ {
return mImpl->comp(); return mImpl->comp();
......
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