Skip to content

L
libcifpp
Commit 83965b9a by maarten
Options

snellere bondmap 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@349 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent e5bd42b4
Showing with 53 additions and 205 deletions
src/BondMap.cpp
...@@ -57,11 +57,21 @@ BondMap::BondMap(const Structure& p) ...@@ -57,11 +57,21 @@ BondMap::BondMap(const Structure& p)
if (compounds.count(c)) if (compounds.count(c))
continue; continue;
if (VERBOSE) if (VERBOSE > 1)
cerr << "Warning: mon_id " << c << " is missing in the chem_comp category" << endl; cerr << "Warning: mon_id " << c << " is missing in the chem_comp category" << endl;
compounds.insert(c); compounds.insert(c);
} }
cif::Progress progress(compounds.size(), "Creating bond map");
// some helper indices to speed things up a bit
map<tuple<string,int,string>,string> atomMapByAsymSeqAndAtom;
for (auto& a: p.atoms())
{
auto key = make_tuple(a.labelAsymId(), a.labelSeqId(), a.labelAtomId());
atomMapByAsymSeqAndAtom[key] = a.id();
}
// first link all residues in a polyseq // first link all residues in a polyseq
string lastAsymID; string lastAsymID;
...@@ -80,16 +90,22 @@ BondMap::BondMap(const Structure& p) ...@@ -80,16 +90,22 @@ BondMap::BondMap(const Structure& p)
continue; continue;
} }
auto c = db["atom_site"].find(cif::Key("label_asym_id") == asymId and cif::Key("label_seq_id") == lastSeqID and cif::Key("label_atom_id") == "C"); auto c = atomMapByAsymSeqAndAtom[make_tuple(asymId, lastSeqID, "C")];
if (c.size() != 1 and VERBOSE > 1) auto n = atomMapByAsymSeqAndAtom[make_tuple(asymId, seqId, "N")];
cerr << "Unexpected number (" << c.size() << ") of atoms with atom ID C in asym_id " << asymId << " with seq id " << lastSeqID << endl;
auto n = db["atom_site"].find(cif::Key("label_asym_id") == asymId and cif::Key("label_seq_id") == seqId and cif::Key("label_atom_id") == "N"); // auto c = db["atom_site"].find(cif::Key("label_asym_id") == asymId and cif::Key("label_seq_id") == lastSeqID and cif::Key("label_atom_id") == "C");
if (n.size() != 1 and VERBOSE > 1) // if (c.size() != 1 and VERBOSE > 1)
cerr << "Unexpected number (" << n.size() << ") of atoms with atom ID N in asym_id " << asymId << " with seq id " << seqId << endl; // cerr << "Unexpected number (" << c.size() << ") of atoms with atom ID C in asym_id " << asymId << " with seq id " << lastSeqID << endl;
//
// auto n = db["atom_site"].find(cif::Key("label_asym_id") == asymId and cif::Key("label_seq_id") == seqId and cif::Key("label_atom_id") == "N");
// if (n.size() != 1 and VERBOSE > 1)
// cerr << "Unexpected number (" << n.size() << ") of atoms with atom ID N in asym_id " << asymId << " with seq id " << seqId << endl;
//
// if (not (c.empty() or n.empty()))
// bindAtoms(c.front()["id"].as<string>(), n.front()["id"].as<string>());
if (not (c.empty() or n.empty())) if (not (c.empty() or n.empty()))
bindAtoms(c.front()["id"].as<string>(), n.front()["id"].as<string>()); bindAtoms(c, n);
lastSeqID = seqId; lastSeqID = seqId;
} }
...@@ -97,36 +113,39 @@ BondMap::BondMap(const Structure& p) ...@@ -97,36 +113,39 @@ BondMap::BondMap(const Structure& p)
for (auto l: db["struct_conn"]) for (auto l: db["struct_conn"])
{ {
string asym1, asym2, atomId1, atomId2; string asym1, asym2, atomId1, atomId2;
int seqId1, seqId2; int seqId1 = 0, seqId2 = 0;
cif::tie(asym1, asym2, atomId1, atomId2, seqId1, seqId2) = cif::tie(asym1, asym2, atomId1, atomId2, seqId1, seqId2) =
l.get("ptnr1_label_asym_id", "ptnr2_label_asym_id", l.get("ptnr1_label_asym_id", "ptnr2_label_asym_id",
"ptnr1_label_atom_id", "ptnr2_label_atom_id", "ptnr1_label_atom_id", "ptnr2_label_atom_id",
"ptnr1_label_seq_id", "ptnr2_label_seq_id"); "ptnr1_label_seq_id", "ptnr2_label_seq_id");
auto a = string a = atomMapByAsymSeqAndAtom[make_tuple(asym1, seqId1, atomId1)];
l["ptnr1_label_seq_id"].empty() ? string b = atomMapByAsymSeqAndAtom[make_tuple(asym2, seqId2, atomId2)];
db["atom_site"].find(cif::Key("label_asym_id") == asym1 and cif::Key("label_atom_id") == atomId1) :
db["atom_site"].find(cif::Key("label_asym_id") == asym1 and cif::Key("label_seq_id") == seqId1 and cif::Key("label_atom_id") == atomId1); // auto a =
// l["ptnr1_label_seq_id"].empty() ?
if (a.size() != 1 and VERBOSE > 1) // db["atom_site"].find(cif::Key("label_asym_id") == asym1 and cif::Key("label_atom_id") == atomId1) :
cerr << "Unexpected number (" << a.size() << ") of atoms for link with asym_id " << asym1 << " seq_id " << seqId1 << " atom_id " << atomId1 << endl; // db["atom_site"].find(cif::Key("label_asym_id") == asym1 and cif::Key("label_seq_id") == seqId1 and cif::Key("label_atom_id") == atomId1);
//
auto b = // if (a.size() != 1 and VERBOSE > 1)
l["ptnr2_label_seq_id"].empty() ? // cerr << "Unexpected number (" << a.size() << ") of atoms for link with asym_id " << asym1 << " seq_id " << seqId1 << " atom_id " << atomId1 << endl;
db["atom_site"].find(cif::Key("label_asym_id") == asym2 and cif::Key("label_atom_id") == atomId2) :
db["atom_site"].find(cif::Key("label_asym_id") == asym2 and cif::Key("label_seq_id") == seqId2 and cif::Key("label_atom_id") == atomId2); // auto b =
// l["ptnr2_label_seq_id"].empty() ?
if (b.size() != 1 and VERBOSE > 1) // db["atom_site"].find(cif::Key("label_asym_id") == asym2 and cif::Key("label_atom_id") == atomId2) :
cerr << "Unexpected number (" << b.size() << ") of atoms for link with asym_id " << asym2 << " seq_id " << seqId2 << " atom_id " << atomId2 << endl; // db["atom_site"].find(cif::Key("label_asym_id") == asym2 and cif::Key("label_seq_id") == seqId2 and cif::Key("label_atom_id") == atomId2);
//
// if (b.size() != 1 and VERBOSE > 1)
// cerr << "Unexpected number (" << b.size() << ") of atoms for link with asym_id " << asym2 << " seq_id " << seqId2 << " atom_id " << atomId2 << endl;
// if (not (a.empty() or b.empty()))
// linkAtoms(a.front()["id"].as<string>(), b.front()["id"].as<string>());
if (not (a.empty() or b.empty())) if (not (a.empty() or b.empty()))
linkAtoms(a.front()["id"].as<string>(), b.front()["id"].as<string>()); linkAtoms(a, b);
} }
// 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);
......
src/DistanceMap.cpp
...@@ -159,12 +159,10 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -159,12 +159,10 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
AddDistancesForAtoms(r, r, dist, 0); AddDistancesForAtoms(r, r, dist, 0);
} }
cif::Progress progress(residues.size(), "Creating distance map"); cif::Progress progress(residues.size() * residues.size(), "Creating distance map");
for (size_t i = 0; i + 1 < residues.size(); ++i) for (size_t i = 0; i + 1 < residues.size(); ++i)
{ {
progress.consumed(1);
auto& ri = *residues[i]; auto& ri = *residues[i];
Point centerI; Point centerI;
...@@ -173,6 +171,8 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -173,6 +171,8 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
for (size_t j = i + 1; j < residues.size(); ++j) for (size_t j = i + 1; j < residues.size(); ++j)
{ {
progress.consumed(1);
auto& rj = *residues[j]; auto& rj = *residues[j];
// first case, no symmetry operations // first case, no symmetry operations
...@@ -211,69 +211,9 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -211,69 +211,9 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
} }
if (minR2 < mMaxDistance) if (minR2 < mMaxDistance)
{
//cout << ri.labelID() << " en " << rj.labelID() << " liggen dicht bij elkaar na symmetrie operatie: " << kbest << endl;
AddDistancesForAtoms(ri, rj, dist, kbest); AddDistancesForAtoms(ri, rj, dist, kbest);
} }
} }
}
// cif::Progress progress(locations.size() - 1, "Creating distance map");
//
// boost::thread_group t;
// size_t N = boost::thread::hardware_concurrency();
// atomic<size_t> next(0);
// mutex m;
//
// for (size_t i = 0; i < N; ++i)
// t.create_thread([&]()
// {
// for (;;)
// {
// size_t i = next++;
//
// if (i >= locations.size())
// break;
//
// clipper::Coord_orth pi = locations[i];
//
// for (size_t j = i + 1; j < locations.size(); ++j)
// {
// // find nearest location based on spacegroup/cell
// double minR2 = numeric_limits<double>::max();
//
// size_t kb = 0;
// for (size_t k = 0; k < mRtOrth.size(); ++k)
// {
// auto& rt = mRtOrth[k];
//
// auto pj = locations[j];
//
// pj = pj.transform(rt);
// double r2 = (pi - pj).lengthsq();
//
// if (minR2 > r2)
// {
// minR2 = r2;
// kb = k;
// }
// }
//
// if (minR2 < mMaxDistanceSQ)
// {
// float d = sqrt(minR2);
// auto key = make_tuple(i, j);
//
// lock_guard<mutex> lock(m);
// dist[key] = make_tuple(d, kb);
// }
// }
//
// progress.consumed(1);
// }
// });
//
// t.join_all();
// Store as a sparse CSR compressed matrix // Store as a sparse CSR compressed matrix
...@@ -303,77 +243,6 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -303,77 +243,6 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
for (size_t ri = lastR; ri < dim; ++ri) for (size_t ri = lastR; ri < dim; ++ri)
mIA.push_back(mA.size()); mIA.push_back(mA.size());
//// debug code
//
// assert(mIA.size() == dim + 1);
// assert(mA.size() == nnz);
// assert(mJA.size() == nnz);
//
//cerr << "nnz: " << nnz << endl;
//
// auto get = [&](size_t i, size_t j)
// {
// if (i > j)
// std::swap(i, j);
//
// for (size_t cix = mIA[i]; cix < mIA[i + 1]; ++cix)
// {
// if (mJA[cix] == j)
// return std::get<0>(mA[cix]);
// }
//
// return 100.f;
// };
//
// auto get2 = [&](size_t ixa, size_t ixb)
// {
// if (ixb < ixa)
// std::swap(ixa, ixb);
//
// int32 L = mIA[ixa];
// int32 R = mIA[ixa + 1] - 1;
//
// while (L <= R)
// {
// int32 i = (L + R) / 2;
//
// if (mJA[i] == ixb)
// return std::get<0>(mA[i]);
//
// if (mJA[i] < ixb)
// L = i + 1;
// else
// R = i - 1;
// }
//
// return 100.f;
// };
//
// // test all values
// for (size_t i = 0; i < dim; ++i)
// for (size_t j = 0; j < dim; ++j)
// {
// float a = get(i, j);
//
// auto ixa = i, ixb = j;
// if (ixb < ixa)
// std::swap(ixa, ixb);
//
// tuple<size_t,size_t> k{ ixa, ixb };
//
// auto ii = dist.find(k);
//
// float b = 100;
//
// if (ii != dist.end())
// b = std::get<0>(ii->second);
//
// assert(a == b);
//
// float c = get2(i, j);
// assert(a == c);
// }
} }
// -------------------------------------------------------------------- // --------------------------------------------------------------------
...@@ -483,52 +352,12 @@ vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const ...@@ -483,52 +352,12 @@ vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const
size_t ixb = mJA[i]; size_t ixb = mJA[i];
Atom b = structure.getAtomById(rIndex.at(ixb)); Atom b = structure.getAtomById(rIndex.at(ixb));
clipper::RTop_orth rt = clipper::RTop_orth::identity();
if (rti > 0) if (rti > 0)
{ result.emplace_back(b.symmetryCopy(mD, mRtOrth.at(rti)));
rt = mRtOrth.at(rti);
result.emplace_back(b.symmetryCopy(mD, rt));
}
else if (rti < 0) else if (rti < 0)
{ result.emplace_back(b.symmetryCopy(mD, mRtOrth.at(-rti).inverse()));
rt = mRtOrth.at(-rti).inverse();
result.emplace_back(b.symmetryCopy(mD, rt));
}
else else
result.emplace_back(b); result.emplace_back(b);
#if 1 //DEBUG
// if (rti != 0)
// cerr << "symmetrie contact " << a.labelID() << " en " << result.back().labelID()
// << " d: " << d
// << " rti: " << rti
// << endl;
auto d2 = Distance(a, result.back());
if (abs(d2 - d) > 0.01)
{
cerr << "Voor a: " << a.location() << " en b: " << b.location() << " => " << result.back().location() << endl;
cerr << "Afstand " << a.labelID() << " en " << result.back().labelID()
<< " is niet gelijk aan verwachtte waarde:"
<< "d: " << d
<< " d2: " << d2
<< " rti: " << rti
<< endl;
rt = rt.inverse();
result.back() = b.symmetryCopy(mD, rt);
d2 = Distance(a, result.back());
cerr << "inverse b: " << result.back().location() << endl;
if (abs(d2 - d) < 0.01)
cerr << "==> But the inverse is correct" << endl;
}
#endif
} }
return result; return result;
......
src/PDB2Cif.cpp
...@@ -4650,7 +4650,7 @@ void PDBFileParser::ParseConnectivtyAnnotation() ...@@ -4650,7 +4650,7 @@ void PDBFileParser::ParseConnectivtyAnnotation()
catch (const invalid_argument&) catch (const invalid_argument&)
{ {
if (VERBOSE) if (VERBOSE)
cerr << "Distance value '" << distance << "' is not a valid float" << endl; cerr << "Distance value '" << distance << "' is not a valid float in LINK record" << endl;
distance.clear(); distance.clear();
} }
} }
......
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