Skip to content

L
libcifpp
Commit 9c636544 by maarten
Options

new Symmetry 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@484 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent 19714ecb
Showing with 452 additions and 220 deletions
include/cif++/DistanceMap.h
...@@ -53,118 +53,4 @@ class DistanceMap ...@@ -53,118 +53,4 @@ class DistanceMap
std::vector<clipper::RTop_orth> mRtOrth; std::vector<clipper::RTop_orth> mRtOrth;
}; };
// --------------------------------------------------------------------
class SymmetryAtomIteratorFactory
{
public:
// SymmetryAtomIteratorFactory(const Structure& p);
SymmetryAtomIteratorFactory(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell)
: mD(DistanceMap::CalculateOffsetForCell(p, spacegroup, cell))
, mRtOrth(DistanceMap::AlternativeSites(spacegroup, cell))
, mSpacegroup(spacegroup)
, mCell(cell) {}
SymmetryAtomIteratorFactory(const SymmetryAtomIteratorFactory&) = delete;
SymmetryAtomIteratorFactory& operator=(const SymmetryAtomIteratorFactory&) = delete;
class SymmetryAtomIterator : public std::iterator<std::forward_iterator_tag, const Atom>
{
public:
typedef std::iterator<std::forward_iterator_tag, const Atom> baseType;
typedef typename baseType::pointer pointer;
typedef typename baseType::reference reference;
SymmetryAtomIterator(const SymmetryAtomIteratorFactory& factory, const Atom& atom)
: m_f(&factory), m_i(0), m_a(atom), m_c(atom) {}
SymmetryAtomIterator(const SymmetryAtomIteratorFactory& factory, const Atom& atom, int)
: SymmetryAtomIterator(factory, atom)
{
m_i = m_f->mRtOrth.size();
}
SymmetryAtomIterator(const SymmetryAtomIterator& iter)
: m_f(iter.m_f), m_i(iter.m_i), m_a(iter.m_a), m_c(iter.m_c) {}
SymmetryAtomIterator& operator=(const SymmetryAtomIterator& iter)
{
if (this != &iter)
{
m_f = iter.m_f;
m_i = iter.m_i;
m_a = iter.m_a;
m_c = iter.m_c;
}
return *this;
}
reference operator*() { return m_c; }
pointer operator->() { return &m_c; }
SymmetryAtomIterator operator++()
{
if (++m_i < m_f->mRtOrth.size())
m_c = m_a.symmetryCopy(m_f->mD, m_f->mRtOrth[m_i]);
return *this;
}
SymmetryAtomIterator operator++(int)
{
SymmetryAtomIterator result(*this);
this->operator++();
return result;
}
bool operator==(const SymmetryAtomIterator& iter) const
{
return m_f == iter.m_f and m_i == iter.m_i;
}
bool operator!=(const SymmetryAtomIterator& iter) const
{
return m_f != iter.m_f or m_i != iter.m_i;
}
private:
const SymmetryAtomIteratorFactory* m_f;
size_t m_i;
Atom m_a, m_c;
};
class SymmetryAtomIteratorRange
{
public:
SymmetryAtomIteratorRange(const SymmetryAtomIteratorFactory& f, const Atom& a)
: m_f(f), m_a(a) {}
SymmetryAtomIterator begin()
{
return SymmetryAtomIterator(m_f, m_a);
}
SymmetryAtomIterator end()
{
return SymmetryAtomIterator(m_f, m_a, 1);
}
private:
const SymmetryAtomIteratorFactory& m_f;
Atom m_a;
};
SymmetryAtomIteratorRange operator()(const Atom& a) const
{
return SymmetryAtomIteratorRange(*this, a);
}
std::string symop_mmcif(const Atom& a) const;
private:
Point mD; // needed to move atoms to center
std::vector<clipper::RTop_orth> mRtOrth;
clipper::Spacegroup mSpacegroup;
clipper::Cell mCell;
};
} }
include/cif++/Symmetry.h 0 → 100644
// copyright
#pragma once
#include "cif++/Structure.h"
namespace mmcif
{
// --------------------------------------------------------------------
// Functions to use when working with symmetry stuff
clipper::Coord_orth CalculateOffsetForCell(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell);
std::vector<clipper::RTop_orth> AlternativeSites(const clipper::Spacegroup& spacegroup, const clipper::Cell& cell);
int32_t GetRotationalIndexNumber(const clipper::Spacegroup& spacegroup, const clipper::Cell& cell,
const clipper::RTop_orth& rt);
// --------------------------------------------------------------------
// To iterate over all near symmetry copies of an atom
class SymmetryAtomIteratorFactory
{
public:
// SymmetryAtomIteratorFactory(const Structure& p);
SymmetryAtomIteratorFactory(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell)
: mD(CalculateOffsetForCell(p, spacegroup, cell))
, mRtOrth(AlternativeSites(spacegroup, cell))
, mSpacegroup(spacegroup)
, mCell(cell) {}
SymmetryAtomIteratorFactory(const SymmetryAtomIteratorFactory&) = delete;
SymmetryAtomIteratorFactory& operator=(const SymmetryAtomIteratorFactory&) = delete;
class SymmetryAtomIterator : public std::iterator<std::forward_iterator_tag, const Atom>
{
public:
typedef std::iterator<std::forward_iterator_tag, const Atom> baseType;
typedef typename baseType::pointer pointer;
typedef typename baseType::reference reference;
SymmetryAtomIterator(const SymmetryAtomIteratorFactory& factory, const Atom& atom)
: m_f(&factory), m_i(0), m_a(atom), m_c(atom) {}
SymmetryAtomIterator(const SymmetryAtomIteratorFactory& factory, const Atom& atom, int)
: SymmetryAtomIterator(factory, atom)
{
m_i = m_f->mRtOrth.size();
}
SymmetryAtomIterator(const SymmetryAtomIterator& iter)
: m_f(iter.m_f), m_i(iter.m_i), m_a(iter.m_a), m_c(iter.m_c) {}
SymmetryAtomIterator& operator=(const SymmetryAtomIterator& iter)
{
if (this != &iter)
{
m_f = iter.m_f;
m_i = iter.m_i;
m_a = iter.m_a;
m_c = iter.m_c;
}
return *this;
}
reference operator*() { return m_c; }
pointer operator->() { return &m_c; }
SymmetryAtomIterator operator++()
{
if (++m_i < m_f->mRtOrth.size())
m_c = m_a.symmetryCopy(m_f->mD, m_f->mRtOrth[m_i]);
return *this;
}
SymmetryAtomIterator operator++(int)
{
SymmetryAtomIterator result(*this);
this->operator++();
return result;
}
bool operator==(const SymmetryAtomIterator& iter) const
{
return m_f == iter.m_f and m_i == iter.m_i;
}
bool operator!=(const SymmetryAtomIterator& iter) const
{
return m_f != iter.m_f or m_i != iter.m_i;
}
private:
const SymmetryAtomIteratorFactory* m_f;
size_t m_i;
Atom m_a, m_c;
};
class SymmetryAtomIteratorRange
{
public:
SymmetryAtomIteratorRange(const SymmetryAtomIteratorFactory& f, const Atom& a)
: m_f(f), m_a(a) {}
SymmetryAtomIterator begin()
{
return SymmetryAtomIterator(m_f, m_a);
}
SymmetryAtomIterator end()
{
return SymmetryAtomIterator(m_f, m_a, 1);
}
private:
const SymmetryAtomIteratorFactory& m_f;
Atom m_a;
};
SymmetryAtomIteratorRange operator()(const Atom& a) const
{
return SymmetryAtomIteratorRange(*this, a);
}
std::string symop_mmcif(const Atom& a) const;
private:
Point mD; // needed to move atoms to center
std::vector<clipper::RTop_orth> mRtOrth;
clipper::Spacegroup mSpacegroup;
clipper::Cell mCell;
};
}
src/DistanceMap.cpp
...@@ -450,60 +450,4 @@ vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const ...@@ -450,60 +450,4 @@ vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const
return result; return result;
} }
// -----------------------------------------------------------------------
string SymmetryAtomIteratorFactory::symop_mmcif(const Atom& a) const
{
string result;
if (not a.isSymmetryCopy())
result = "1_555";
else
{
auto rtop_o = a.symop();
for (int i = 0; i < mSpacegroup.num_symops(); ++i)
{
const auto& symop = mSpacegroup.symop(i);
for (int u: { -1, 0, 1})
for (int v: { -1, 0, 1})
for (int w: { -1, 0, 1})
{
// if (i == 0 and u == 0 and v == 0 and w == 0)
// continue;
auto rtop = clipper::RTop_frac(
symop.rot(), symop.trn() + clipper::Vec3<>(u, v, w)
).rtop_orth(mCell);
if (rtop.rot().equals(rtop_o.rot(), 0.00001) and rtop.trn().equals(rtop_o.trn(), 0.000001))
{
// gotcha
auto rtop_f = rtop.rtop_frac(mCell);
uint32_t t[3] = {
static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[0]))),
static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[1]))),
static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[2])))
};
if (t[0] > 9 or t[1] > 9 or t[2] > 9)
throw runtime_error("Symmetry operation has an out-of-range translation.");
result += to_string(1 + i) + "_"
+ to_string(t[0])
+ to_string(t[1])
+ to_string(t[2]);
return result;
}
}
}
}
return result;
}
} }
src/Symmetry.cpp 0 → 100644
// copyright
#include "cif++/Config.h"
#include <atomic>
#include <mutex>
#include <boost/thread.hpp>
#include "cif++/Symmetry.h"
#include "cif++/CifUtils.h"
using namespace std;
namespace mmcif
{
// --------------------------------------------------------------------
clipper::Coord_orth CalculateOffsetForCell(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell)
{
auto& atoms = p.atoms();
size_t dim = atoms.size();
vector<clipper::Coord_orth> locations;
locations.reserve(dim);
// bounding box
Point pMin(numeric_limits<float>::max(), numeric_limits<float>::max(), numeric_limits<float>::max()),
pMax(numeric_limits<float>::min(), numeric_limits<float>::min(), numeric_limits<float>::min());
for (auto& atom: atoms)
{
auto p = atom.location();
locations.push_back(p);
if (pMin.mX > p.mX)
pMin.mX = p.mX;
if (pMin.mY > p.mY)
pMin.mY = p.mY;
if (pMin.mZ > p.mZ)
pMin.mZ = p.mZ;
if (pMax.mX < p.mX)
pMax.mX = p.mX;
if (pMax.mY < p.mY)
pMax.mY = p.mY;
if (pMax.mZ < p.mZ)
pMax.mZ = p.mZ;
};
// correct locations so that the median of x, y and z are inside the cell
vector<float> c(dim);
auto median = [&]()
{
return dim % 1 == 0
? c[dim / 2]
: (c[dim / 2 - 1] + c[dim / 2]) / 2;
};
transform(locations.begin(), locations.end(), c.begin(), [](auto& l) { return l[0]; });
sort(c.begin(), c.end());
float mx = median();
transform(locations.begin(), locations.end(), c.begin(), [](auto& l) { return l[1]; });
sort(c.begin(), c.end());
float my = median();
transform(locations.begin(), locations.end(), c.begin(), [](auto& l) { return l[2]; });
sort(c.begin(), c.end());
float mz = median();
if (cif::VERBOSE > 1)
cerr << "median position of atoms: " << Point(mx, my, mz) << endl;
auto calculateD = [&](float m, float c)
{
float d = 0;
while (m + d < -(c / 2))
d += c;
while (m + d > (c / 2))
d -= c;
return d;
};
Point D;
D.mX = calculateD(mx, cell.a());
D.mY = calculateD(my, cell.b());
D.mZ = calculateD(mz, cell.c());
if (D.mX != 0 or D.mY != 0 or D.mZ != 0)
{
if (cif::VERBOSE)
cerr << "moving coorinates by " << D.mX << ", " << D.mY << " and " << D.mZ << endl;
}
return D;
}
// --------------------------------------------------------------------
vector<clipper::RTop_orth> AlternativeSites(const clipper::Spacegroup& spacegroup,
const clipper::Cell& cell)
{
vector<clipper::RTop_orth> result;
// to make the operation at index 0 equal to identity
result.push_back(clipper::RTop_orth::identity());
for (int i = 0; i < spacegroup.num_symops(); ++i)
{
const auto& symop = spacegroup.symop(i);
for (int u: { -1, 0, 1})
for (int v: { -1, 0, 1})
for (int w: { -1, 0, 1})
{
if (i == 0 and u == 0 and v == 0 and w == 0)
continue;
auto rtop = clipper::RTop_frac(
symop.rot(), symop.trn() + clipper::Vec3<>(u, v, w)
).rtop_orth(cell);
result.push_back(move(rtop));
}
}
return result;
}
// -----------------------------------------------------------------------
string SymmetryAtomIteratorFactory::symop_mmcif(const Atom& a) const
{
string result;
if (not a.isSymmetryCopy())
result = "1_555";
else
{
auto rtop_o = a.symop();
for (int i = 0; i < mSpacegroup.num_symops(); ++i)
{
const auto& symop = mSpacegroup.symop(i);
for (int u: { -1, 0, 1})
for (int v: { -1, 0, 1})
for (int w: { -1, 0, 1})
{
// if (i == 0 and u == 0 and v == 0 and w == 0)
// continue;
auto rtop = clipper::RTop_frac(
symop.rot(), symop.trn() + clipper::Vec3<>(u, v, w)
).rtop_orth(mCell);
if (rtop.rot().equals(rtop_o.rot(), 0.00001) and rtop.trn().equals(rtop_o.trn(), 0.000001))
{
// gotcha
auto rtop_f = rtop.rtop_frac(mCell);
uint32_t t[3] = {
static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[0]))),
static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[1]))),
static_cast<uint32_t>(5 + static_cast<int>(rint(rtop_f.trn()[2])))
};
if (t[0] > 9 or t[1] > 9 or t[2] > 9)
throw runtime_error("Symmetry operation has an out-of-range translation.");
result += to_string(1 + i) + "_"
+ to_string(t[0])
+ to_string(t[1])
+ to_string(t[2]);
return result;
}
}
}
}
return result;
}
}
tools/symop-map-generator.cpp
...@@ -8,6 +8,7 @@ clang++ -I ~/my-clipper/include -L ~/my-clipper/lib -o symop-map-generator symop ...@@ -8,6 +8,7 @@ clang++ -I ~/my-clipper/include -L ~/my-clipper/lib -o symop-map-generator symop
#include <cassert> #include <cassert>
#include <iostream> #include <iostream>
#include <iomanip>
#include <fstream> #include <fstream>
#include <regex> #include <regex>
...@@ -22,94 +23,141 @@ std::regex kNameRx(R"(^(\d+) +(\d+) +(\d+) +(\S+) +(\S+) +(\S+) +'([^']+)'( +'([ ...@@ -22,94 +23,141 @@ std::regex kNameRx(R"(^(\d+) +(\d+) +(\d+) +(\S+) +(\S+) +(\S+) +'([^']+)'( +'([
class SymopParser class SymopParser
{ {
public: public:
SymopParser(const std::string& s) SymopParser() {}
: m_s(s)
, m_p(s.begin())
, m_e(s.end())
{
m_lookahead = next_char();
}
void parse() array<int,15> parse(const string& s)
{ {
m_p = s.begin();
m_e = s.end();
m_lookahead = next_token();
parsepart(0); parsepart(0);
match(','); match((Token)',');
parsepart(1); parsepart(1);
match(','); match((Token)',');
parsepart(2); parsepart(2);
if (m_lookahead != 0)
if (m_lookahead != 0 or m_p != m_e)
throw runtime_error("symmetry expression contains more data than expected"); throw runtime_error("symmetry expression contains more data than expected");
return {
m_rot[0][0], m_rot[0][1], m_rot[0][2],
m_rot[1][0], m_rot[1][1], m_rot[1][2],
m_rot[2][0], m_rot[2][1], m_rot[2][2],
m_trn[0][0], m_trn[0][1],
m_trn[1][0], m_trn[1][1],
m_trn[2][0], m_trn[2][1]
};
} }
private: private:
char next_char() enum Token : int { Eof = 0, Number = 256, XYZ };
string to_string(Token t)
{
switch (t)
{
case Eof: return "end of expression";
case Number: return "number";
case XYZ: return "'x', 'y' or 'z'";
default:
if (isprint(t))
return string({'\'', static_cast<char>(t), '\''});
return "invalid character " + std::to_string(static_cast<int>(t));
}
}
Token next_token()
{ {
char result = 0; Token result = Eof;
while (m_p != m_e) while (m_p != m_e)
{ {
char ch = *m_p++; char ch = *m_p++;
if (ch == ' ') if (ch == ' ')
continue; continue;
result = ch; switch (ch)
{
case 'x':
case 'X':
result = XYZ;
m_nr = 0;
break;
case 'y':
case 'Y':
result = XYZ;
m_nr = 1;
break;
case 'z':
case 'Z':
result = XYZ;
m_nr = 2;
break;
default:
if (isdigit(ch))
{
m_nr = ch - '0';
result = Number;
}
else
result = (Token)ch;
break;
}
break; break;
} }
return result; return result;
} }
void retract() void match(Token token)
{
--m_p;
}
void match(char c)
{ {
if (m_lookahead != c) if (m_lookahead != token)
throw runtime_error("Unexpected character " + m_lookahead + " expected " + c); throw runtime_error("Unexpected character " + to_string(m_lookahead) + " expected " + to_string(token));
m_lookahead = next_char(); m_lookahead = next_token();
} }
void parsepart(int xyz) void parsepart(int row)
{ {
if (isdigit(m_lookahead)) bool first = false;
do
{ {
parserational(xyz); int sign = m_lookahead == '-' ? -1 : 1;
if (m_lookahead == '-' or m_lookahead == '+')
match(m_lookahead);
switch (m_lookahead) if (m_lookahead == Number)
{ {
case '-': m_trn[row][0] = sign * m_nr;
match('-'); match(Number);
break;
}
} match((Token)'/');
}
if (m_lookahead == '+')
{
match('+');
}
else if (m_lookahead = '-')
{
m_trn[row][1] = m_nr;
match(Number);
}
else
{
m_rot[row][m_nr] = sign;
match(XYZ);
}
} }
while (m_lookahead == '+' or m_lookahead == '-');
} }
tuple<int,int> parserational(); Token m_lookahead;
char m_nr;
char m_lookahead;
string m_s; string m_s;
string::iterator m_p, m_e; string::const_iterator m_p, m_e;
int m_rot[3][3] = {}; int m_rot[3][3] = {};
int m_trn[3] = {}; int m_trn[3][2] = {};
}; };
int main() int main()
...@@ -129,6 +177,21 @@ int main() ...@@ -129,6 +177,21 @@ int main()
clipper::Spacegroup spacegroup; clipper::Spacegroup spacegroup;
enum class State { Initial, InSpacegroup, Error } state = State::Initial; enum class State { Initial, InSpacegroup, Error } state = State::Initial;
int symopnr = 0;
cout << R"(
// This file was generated from $CLIBD/symop.lib
struct SymopNrTable
{
uint8_t spacegroupNr;
uint8_t rotationalNr:7;
int8_t rot[3][3]:2;
uint8_t trn[3][2]:4;
};
kSymopNrTable[] = {
)" << endl;
while (getline(file, line)) while (getline(file, line))
{ {
...@@ -144,7 +207,21 @@ int main() ...@@ -144,7 +207,21 @@ int main()
if (state == State::Error) if (state == State::Error)
continue; continue;
try
{
SymopParser p;
auto symop = p.parse(line);
cout << " { " << setw(3) << spacegroup.spacegroup_number() << '\t' << symopnr++;
for (auto i: symop)
cout << ',' << setw(2) << i;
cout << " }," << endl;
}
catch (const exception& e)
{
cerr << line << endl
<< e.what() << endl;
}
continue; continue;
...@@ -170,7 +247,10 @@ int main() ...@@ -170,7 +247,10 @@ int main()
try try
{ {
spacegroup = clipper::Spacegroup(clipper::Spgr_descr(stoi(m[1].str()))); spacegroup = clipper::Spacegroup(clipper::Spgr_descr(stoi(m[1].str())));
symopnr = 1;
cout << " // " << spacegroup.symbol_hm() << endl;
if (spacegroup.num_symops() != stoi(m[2])) if (spacegroup.num_symops() != stoi(m[2]))
{ {
cerr << line << endl cerr << line << endl
...@@ -191,7 +271,7 @@ int main() ...@@ -191,7 +271,7 @@ int main()
} }
} }
} }
cout << "};" << endl;
} }
catch (const exception& ex) catch (const exception& ex)
{ {
......
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