Skip to content

L
libcifpp
Commit 9aa8a223 by Maarten L. Hekkelman
Options

symmetry work

parent fb59adcf
Showing with 887 additions and 470 deletions
include/cif++/matrix.hpp 0 → 100644
/*-
* SPDX-License-Identifier: BSD-2-Clause
*
* Copyright (c) 2023 NKI/AVL, Netherlands Cancer Institute
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* 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
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#pragma once
#include <array>
#include <cassert>
#include <cstdint>
#include <ostream>
#include <type_traits>
#include <vector>
namespace cif
{
// --------------------------------------------------------------------
// We're using expression templates here
template <typename M>
class matrix_expression
{
public:
constexpr uint32_t dim_m() const { return static_cast<const M &>(*this).dim_m(); }
constexpr uint32_t dim_n() const { return static_cast<const M &>(*this).dim_n(); }
constexpr auto &operator()(uint32_t i, uint32_t j)
{
return static_cast<M &>(*this).operator()(i, j);
}
constexpr auto operator()(uint32_t i, uint32_t j) const
{
return static_cast<const M &>(*this).operator()(i, j);
}
friend std::ostream &operator<<(std::ostream &os, const matrix_expression &m)
{
os << '[';
for (size_t i = 0; i < m.dim_m(); ++i)
{
os << '[';
for (size_t j = 0; j < m.dim_n(); ++j)
{
os << m(i, j);
if (j + 1 < m.dim_n())
os << ", ";
}
if (i + 1 < m.dim_m())
os << ", ";
os << ']';
}
os << ']';
return os;
}
};
// --------------------------------------------------------------------
// matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
// element m i,j is mapped to [i * n + j] and thus storage is row major
template <typename F = float>
class matrix : public matrix_expression<matrix<F>>
{
public:
using value_type = F;
template <typename M2>
matrix(const matrix_expression<M2> &m)
: m_m(m.dim_m())
, m_n(m.dim_n())
, m_data(m_m * m_n)
{
for (uint32_t i = 0; i < m_m; ++i)
{
for (uint32_t j = 0; j < m_n; ++j)
operator()(i, j) = m(i, j);
}
}
matrix(size_t m, size_t n, value_type v = 0)
: m_m(m)
, m_n(n)
, m_data(m_m * m_n)
{
std::fill(m_data.begin(), m_data.end(), v);
}
matrix() = default;
matrix(matrix &&m) = default;
matrix(const matrix &m) = default;
matrix &operator=(matrix &&m) = default;
matrix &operator=(const matrix &m) = default;
constexpr size_t dim_m() const { return m_m; }
constexpr size_t dim_n() const { return m_n; }
constexpr value_type operator()(size_t i, size_t j) const
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
constexpr value_type &operator()(size_t i, size_t j)
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
private:
size_t m_m = 0, m_n = 0;
std::vector<value_type> m_data;
};
// --------------------------------------------------------------------
// special case, 3x3 matrix
template <typename F, size_t M, size_t N>
class matrix_fixed : public matrix_expression<matrix_fixed<F, M, N>>
{
public:
using value_type = F;
template <typename M2>
matrix_fixed(const M2 &m)
{
assert(M == m.dim_m() and N == m.dim_n());
for (uint32_t i = 0; i < M; ++i)
{
for (uint32_t j = 0; j < N; ++j)
operator()(i, j) = m(i, j);
}
}
matrix_fixed(value_type v = 0)
{
std::fill(m_data.begin(), m_data.end(), v);
}
matrix_fixed(matrix_fixed &&m) = default;
matrix_fixed(const matrix_fixed &m) = default;
matrix_fixed &operator=(matrix_fixed &&m) = default;
matrix_fixed &operator=(const matrix_fixed &m) = default;
constexpr size_t dim_m() const { return M; }
constexpr size_t dim_n() const { return N; }
constexpr value_type operator()(size_t i, size_t j) const
{
assert(i < M);
assert(j < N);
return m_data[i * N + j];
}
constexpr value_type &operator()(size_t i, size_t j)
{
assert(i < M);
assert(j < N);
return m_data[i * N + j];
}
private:
std::array<value_type, M * N> m_data;
};
template <typename F>
using matrix3x3 = matrix_fixed<F, 3, 3>;
template <typename F>
using matrix4x4 = matrix_fixed<F, 4, 4>;
// --------------------------------------------------------------------
template <typename F = float>
class symmetric_matrix : public matrix_expression<symmetric_matrix<F>>
{
public:
using value_type = F;
symmetric_matrix(uint32_t n, value_type v = 0)
: m_n(n)
, m_data((m_n * (m_n + 1)) / 2)
{
std::fill(m_data.begin(), m_data.end(), v);
}
symmetric_matrix() = default;
symmetric_matrix(symmetric_matrix &&m) = default;
symmetric_matrix(const symmetric_matrix &m) = default;
symmetric_matrix &operator=(symmetric_matrix &&m) = default;
symmetric_matrix &operator=(const symmetric_matrix &m) = default;
constexpr uint32_t dim_m() const { return m_n; }
constexpr uint32_t dim_n() const { return m_n; }
constexpr value_type operator()(uint32_t i, uint32_t j) const
{
return i < j
? m_data[(j * (j + 1)) / 2 + i]
: m_data[(i * (i + 1)) / 2 + j];
}
constexpr value_type &operator()(uint32_t i, uint32_t j)
{
if (i > j)
std::swap(i, j);
assert(j < m_n);
return m_data[(j * (j + 1)) / 2 + i];
}
private:
uint32_t m_n;
std::vector<value_type> m_data;
};
// --------------------------------------------------------------------
template <typename F, size_t M>
class symmetric_matrix_fixed : public matrix_expression<symmetric_matrix_fixed<F, M>>
{
public:
using value_type = F;
symmetric_matrix_fixed(value_type v = 0)
{
std::fill(m_data.begin(), m_data.end(), v);
}
symmetric_matrix_fixed(symmetric_matrix_fixed &&m) = default;
symmetric_matrix_fixed(const symmetric_matrix_fixed &m) = default;
symmetric_matrix_fixed &operator=(symmetric_matrix_fixed &&m) = default;
symmetric_matrix_fixed &operator=(const symmetric_matrix_fixed &m) = default;
constexpr uint32_t dim_m() const { return M; }
constexpr uint32_t dim_n() const { return M; }
constexpr value_type operator()(uint32_t i, uint32_t j) const
{
return i < j
? m_data[(j * (j + 1)) / 2 + i]
: m_data[(i * (i + 1)) / 2 + j];
}
constexpr value_type &operator()(uint32_t i, uint32_t j)
{
if (i > j)
std::swap(i, j);
assert(j < M);
return m_data[(j * (j + 1)) / 2 + i];
}
private:
std::array<value_type, (M * (M + 1)) / 2> m_data;
};
template <typename F>
using symmetric_matrix3x3 = symmetric_matrix_fixed<F, 3>;
template <typename F>
using symmetric_matrix4x4 = symmetric_matrix_fixed<F, 4>;
// --------------------------------------------------------------------
template <typename F = float>
class identity_matrix : public matrix_expression<identity_matrix<F>>
{
public:
using value_type = F;
identity_matrix(uint32_t n)
: m_n(n)
{
}
constexpr uint32_t dim_m() const { return m_n; }
constexpr uint32_t dim_n() const { return m_n; }
constexpr value_type operator()(uint32_t i, uint32_t j) const
{
return i == j ? 1 : 0;
}
private:
uint32_t m_n;
};
// --------------------------------------------------------------------
// matrix functions, implemented as expression templates
template <typename M1, typename M2>
class matrix_subtraction : public matrix_expression<matrix_subtraction<M1, M2>>
{
public:
matrix_subtraction(const M1 &m1, const M2 &m2)
: m_m1(m1)
, m_m2(m2)
{
assert(m_m1.dim_m() == m_m2.dim_m());
assert(m_m1.dim_n() == m_m2.dim_n());
}
constexpr uint32_t dim_m() const { return m_m1.dim_m(); }
constexpr uint32_t dim_n() const { return m_m1.dim_n(); }
constexpr auto operator()(uint32_t i, uint32_t j) const
{
return m_m1(i, j) - m_m2(i, j);
}
private:
const M1 &m_m1;
const M2 &m_m2;
};
template <typename M1, typename M2>
auto operator-(const matrix_expression<M1> &m1, const matrix_expression<M2> &m2)
{
return matrix_subtraction(m1, m2);
}
template <typename M, typename F>
class matrix_multiplication : public matrix_expression<matrix_multiplication<M, F>>
{
public:
using value_type = F;
matrix_multiplication(const M &m, value_type v)
: m_m(m)
, m_v(v)
{
}
constexpr uint32_t dim_m() const { return m_m.dim_m(); }
constexpr uint32_t dim_n() const { return m_m.dim_n(); }
constexpr auto operator()(uint32_t i, uint32_t j) const
{
return m_m(i, j) * m_v;
}
private:
const M &m_m;
value_type m_v;
};
template <typename M, typename F>
auto operator*(const matrix_expression<M> &m, F v)
{
return matrix_multiplication(m, v);
}
// --------------------------------------------------------------------
template <typename M>
auto determinant(const M &m);
template <typename F = float>
auto determinant(const matrix3x3<F> &m)
{
return (m(0, 0) * (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) +
m(0, 1) * (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) +
m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0)));
}
template <typename M>
M inverse(const M &m);
template <typename F = float>
matrix3x3<F> inverse(const matrix3x3<F> &m)
{
F det = determinant(m);
matrix3x3<F> result;
result(0, 0) = (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) / det;
result(1, 0) = (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) / det;
result(2, 0) = (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0)) / det;
result(0, 1) = (m(2, 1) * m(0, 2) - m(2, 2) * m(0, 1)) / det;
result(1, 1) = (m(2, 2) * m(0, 0) - m(2, 0) * m(0, 2)) / det;
result(2, 1) = (m(2, 0) * m(0, 1) - m(2, 1) * m(0, 0)) / det;
result(0, 2) = (m(0, 1) * m(1, 2) - m(0, 2) * m(1, 1)) / det;
result(1, 2) = (m(0, 2) * m(1, 0) - m(0, 0) * m(1, 2)) / det;
result(2, 2) = (m(0, 0) * m(1, 1) - m(0, 1) * m(1, 0)) / det;
return result;
}
// --------------------------------------------------------------------
template <typename M>
class matrix_cofactors : public matrix_expression<matrix_cofactors<M>>
{
public:
matrix_cofactors(const M &m)
: m_m(m)
{
}
constexpr uint32_t dim_m() const { return m_m.dim_m(); }
constexpr uint32_t dim_n() const { return m_m.dim_n(); }
constexpr auto operator()(uint32_t i, uint32_t j) const
{
const size_t ixs[4][3] = {
{ 1, 2, 3 },
{ 0, 2, 3 },
{ 0, 1, 3 },
{ 0, 1, 2 }
};
const size_t *ix = ixs[i];
const size_t *iy = ixs[j];
auto result =
m_m(ix[0], iy[0]) * m_m(ix[1], iy[1]) * m_m(ix[2], iy[2]) +
m_m(ix[0], iy[1]) * m_m(ix[1], iy[2]) * m_m(ix[2], iy[0]) +
m_m(ix[0], iy[2]) * m_m(ix[1], iy[0]) * m_m(ix[2], iy[1]) -
m_m(ix[0], iy[2]) * m_m(ix[1], iy[1]) * m_m(ix[2], iy[0]) -
m_m(ix[0], iy[1]) * m_m(ix[1], iy[0]) * m_m(ix[2], iy[2]) -
m_m(ix[0], iy[0]) * m_m(ix[1], iy[2]) * m_m(ix[2], iy[1]);
return i + j % 2 == 1 ? -result : result;
}
private:
const M &m_m;
};
} // namespace cif
include/cif++/point.hpp
......@@ -718,7 +718,6 @@ template <int N>
class spherical_dots
{
public:
constexpr static int P = 2 * N * 1;
using array_type = typename std::array<point, P>;
......
include/cif++/symmetry.hpp
......@@ -27,6 +27,8 @@
#pragma once
#include "cif++/exports.hpp"
#include "cif++/matrix.hpp"
#include "cif++/point.hpp"
#include <array>
#include <cstdint>
......@@ -60,20 +62,20 @@ extern CIFPP_EXPORT const std::size_t kNrOfSpaceGroups;
struct symop_data
{
constexpr symop_data(const std::array<int, 15> &data)
: m_packed((data[0] bitand 0x03ULL) << 34 bitor
(data[1] bitand 0x03ULL) << 32 bitor
(data[2] bitand 0x03ULL) << 30 bitor
(data[3] bitand 0x03ULL) << 28 bitor
(data[4] bitand 0x03ULL) << 26 bitor
(data[5] bitand 0x03ULL) << 24 bitor
(data[6] bitand 0x03ULL) << 22 bitor
(data[7] bitand 0x03ULL) << 20 bitor
(data[8] bitand 0x03ULL) << 18 bitor
(data[9] bitand 0x07ULL) << 15 bitor
: m_packed((data[0] bitand 0x03ULL) << 34 bitor
(data[1] bitand 0x03ULL) << 32 bitor
(data[2] bitand 0x03ULL) << 30 bitor
(data[3] bitand 0x03ULL) << 28 bitor
(data[4] bitand 0x03ULL) << 26 bitor
(data[5] bitand 0x03ULL) << 24 bitor
(data[6] bitand 0x03ULL) << 22 bitor
(data[7] bitand 0x03ULL) << 20 bitor
(data[8] bitand 0x03ULL) << 18 bitor
(data[9] bitand 0x07ULL) << 15 bitor
(data[10] bitand 0x07ULL) << 12 bitor
(data[11] bitand 0x07ULL) << 9 bitor
(data[12] bitand 0x07ULL) << 6 bitor
(data[13] bitand 0x07ULL) << 3 bitor
(data[11] bitand 0x07ULL) << 9 bitor
(data[12] bitand 0x07ULL) << 6 bitor
(data[13] bitand 0x07ULL) << 3 bitor
(data[14] bitand 0x07ULL) << 0)
{
}
......@@ -156,8 +158,189 @@ extern CIFPP_EXPORT const symop_datablock kSymopNrTable[];
extern CIFPP_EXPORT const std::size_t kSymopNrTableSize;
// --------------------------------------------------------------------
// Some more symmetry related stuff here.
class datablock;
class cell;
class spacegroup;
class rtop;
class sym_op;
// --------------------------------------------------------------------
// class cell
class cell
{
public:
cell(float a, float b, float c, float alpha = 90.f, float beta = 90.f, float gamma = 90.f);
cell(const datablock &db);
float get_a() const { return m_a; }
float get_b() const { return m_b; }
float get_c() const { return m_c; }
float get_alpha() const { return m_alpha; }
float get_beta() const { return m_beta; }
float get_gamma() const { return m_gamma; }
matrix3x3<float> get_orthogonal_matrix() const { return m_orthogonal; }
matrix3x3<float> get_fractional_matrix() const { return m_fractional; }
private:
void init();
float m_a, m_b, m_c, m_alpha, m_beta, m_gamma;
matrix3x3<float> m_orthogonal, m_fractional;
};
/// @brief A class that encapsulates the symmetry operations as used in PDB files, i.e. a rotational number and a translation vector
class sym_op
{
public:
sym_op(uint8_t nr = 1, uint8_t ta = 5, uint8_t tb = 5, uint8_t tc = 5)
: m_nr(nr)
, m_ta(ta)
, m_tb(tb)
, m_tc(tc)
{
}
explicit sym_op(std::string_view s);
sym_op(const sym_op &) = default;
sym_op(sym_op &&) = default;
sym_op &operator=(const sym_op &) = default;
sym_op &operator=(sym_op &&) = default;
constexpr bool is_identity() const
{
return m_nr == 1 and m_ta == 5 and m_tb == 5 and m_tc == 5;
}
explicit constexpr operator bool() const
{
return not is_identity();
}
std::string string() const;
friend class spacegroup;
private:
uint8_t m_nr;
uint8_t m_ta, m_tb, m_tc;
};
namespace literals
{
inline sym_op operator""_symop(const char *text, size_t length)
{
return sym_op({ text, length });
}
} // namespace literals
class transformation
{
public:
transformation(const symop_data &data);
transformation(const matrix3x3<float> &r, const cif::point &t)
: m_rotation(r)
, m_translation(t)
{
}
transformation(const transformation &) = default;
transformation(transformation &&) = default;
transformation &operator=(const transformation &) = default;
transformation &operator=(transformation &&) = default;
point operator()(const cell &c, const point &pt) const;
private:
matrix3x3<float> m_rotation;
point m_translation;
};
// --------------------------------------------------------------------
int get_space_group_number(std::string_view spacegroup); // alternative for clipper's parsing code, using space_group_name::full
int get_space_group_number(std::string_view spacegroup, space_group_name type); // alternative for clipper's parsing code
class spacegroup : public std::vector<transformation>
{
public:
spacegroup(std::string_view name)
: spacegroup(get_space_group_number(name))
{
}
spacegroup(std::string_view name, space_group_name type)
: spacegroup(get_space_group_number(name, type))
{
}
spacegroup(int nr);
int get_nr() const { return m_nr; }
std::string get_name() const;
point operator()(const point &pt, const cell &c, sym_op symop);
private:
int m_nr;
size_t m_index;
};
// --------------------------------------------------------------------
int get_space_group_number(std::string spacegroup); // alternative for clipper's parsing code, using space_group_name::full
int get_space_group_number(std::string spacegroup, space_group_name type); // alternative for clipper's parsing code
// class rtop
// {
// public:
// rtop(const spacegroup &sg, const cell &c, int nr);
// friend rtop operator+(rtop rt, cif::point t);
// friend cif::point operator*(cif::point p, rtop rt);
// private:
// cell m_c;
// cif::quaternion m_q;
// cif::point m_t;
// };
// class spacegroup
// {
// public:
// spacegroup(const cif::datablock &db);
// private:
// std::vector<
// };
static_assert(sizeof(sym_op) == 4, "Sym_op should be four bytes");
// --------------------------------------------------------------------
// Symmetry operations on points
template <typename T>
inline point_type<T> operator*(const matrix3x3<T> &m, const point_type<T> &pt)
{
return point_type(m(0, 0) * pt.m_x + m(0, 1) * pt.m_y + m(0, 2) * pt.m_z,
m(1, 0) * pt.m_x + m(1, 1) * pt.m_y + m(1, 2) * pt.m_z,
m(2, 0) * pt.m_x + m(2, 1) * pt.m_y + m(2, 2) * pt.m_z);
}
inline point orthogonal(const point &pt, const cell &c)
{
return c.get_orthogonal_matrix() * pt;
}
inline point fractional(const point &pt, const cell &c)
{
return c.get_fractional_matrix() * pt;
}
} // namespace cif
src/point.cpp
......@@ -25,6 +25,7 @@
*/
#include "cif++/point.hpp"
#include "cif++/matrix.hpp"
#include <cassert>
#include <random>
......@@ -33,245 +34,6 @@ namespace cif
{
// --------------------------------------------------------------------
// We're using expression templates here
template <typename M>
class MatrixExpression
{
public:
uint32_t dim_m() const { return static_cast<const M &>(*this).dim_m(); }
uint32_t dim_n() const { return static_cast<const M &>(*this).dim_n(); }
double &operator()(uint32_t i, uint32_t j)
{
return static_cast<M &>(*this).operator()(i, j);
}
double operator()(uint32_t i, uint32_t j) const
{
return static_cast<const M &>(*this).operator()(i, j);
}
};
// --------------------------------------------------------------------
// matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
// element m i,j is mapped to [i * n + j] and thus storage is row major
class Matrix : public MatrixExpression<Matrix>
{
public:
template <typename M2>
Matrix(const MatrixExpression<M2> &m)
: m_m(m.dim_m())
, m_n(m.dim_n())
, m_data(m_m * m_n)
{
for (uint32_t i = 0; i < m_m; ++i)
{
for (uint32_t j = 0; j < m_n; ++j)
operator()(i, j) = m(i, j);
}
}
Matrix(size_t m, size_t n, double v = 0)
: m_m(m)
, m_n(n)
, m_data(m_m * m_n)
{
std::fill(m_data.begin(), m_data.end(), v);
}
Matrix() = default;
Matrix(Matrix &&m) = default;
Matrix(const Matrix &m) = default;
Matrix &operator=(Matrix &&m) = default;
Matrix &operator=(const Matrix &m) = default;
size_t dim_m() const { return m_m; }
size_t dim_n() const { return m_n; }
double operator()(size_t i, size_t j) const
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
double &operator()(size_t i, size_t j)
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
private:
size_t m_m = 0, m_n = 0;
std::vector<double> m_data;
};
// --------------------------------------------------------------------
class SymmetricMatrix : public MatrixExpression<SymmetricMatrix>
{
public:
SymmetricMatrix(uint32_t n, double v = 0)
: m_n(n)
, m_data((m_n * (m_n + 1)) / 2)
{
std::fill(m_data.begin(), m_data.end(), v);
}
SymmetricMatrix() = default;
SymmetricMatrix(SymmetricMatrix &&m) = default;
SymmetricMatrix(const SymmetricMatrix &m) = default;
SymmetricMatrix &operator=(SymmetricMatrix &&m) = default;
SymmetricMatrix &operator=(const SymmetricMatrix &m) = default;
uint32_t dim_m() const { return m_n; }
uint32_t dim_n() const { return m_n; }
double operator()(uint32_t i, uint32_t j) const
{
return i < j
? m_data[(j * (j + 1)) / 2 + i]
: m_data[(i * (i + 1)) / 2 + j];
}
double &operator()(uint32_t i, uint32_t j)
{
if (i > j)
std::swap(i, j);
assert(j < m_n);
return m_data[(j * (j + 1)) / 2 + i];
}
private:
uint32_t m_n;
std::vector<double> m_data;
};
class IdentityMatrix : public MatrixExpression<IdentityMatrix>
{
public:
IdentityMatrix(uint32_t n)
: m_n(n)
{
}
uint32_t dim_m() const { return m_n; }
uint32_t dim_n() const { return m_n; }
double operator()(uint32_t i, uint32_t j) const
{
return i == j ? 1 : 0;
}
private:
uint32_t m_n;
};
// --------------------------------------------------------------------
// matrix functions, implemented as expression templates
template <typename M1, typename M2>
class MatrixSubtraction : public MatrixExpression<MatrixSubtraction<M1, M2>>
{
public:
MatrixSubtraction(const M1 &m1, const M2 &m2)
: m_m1(m1)
, m_m2(m2)
{
assert(m_m1.dim_m() == m_m2.dim_m());
assert(m_m1.dim_n() == m_m2.dim_n());
}
uint32_t dim_m() const { return m_m1.dim_m(); }
uint32_t dim_n() const { return m_m1.dim_n(); }
double operator()(uint32_t i, uint32_t j) const
{
return m_m1(i, j) - m_m2(i, j);
}
private:
const M1 &m_m1;
const M2 &m_m2;
};
template <typename M1, typename M2>
MatrixSubtraction<M1, M2> operator-(const MatrixExpression<M1> &m1, const MatrixExpression<M2> &m2)
{
return MatrixSubtraction(*static_cast<const M1 *>(&m1), *static_cast<const M2 *>(&m2));
}
template <typename M>
class MatrixMultiplication : public MatrixExpression<MatrixMultiplication<M>>
{
public:
MatrixMultiplication(const M &m, double v)
: m_m(m)
, m_v(v)
{
}
uint32_t dim_m() const { return m_m.dim_m(); }
uint32_t dim_n() const { return m_m.dim_n(); }
double operator()(uint32_t i, uint32_t j) const
{
return m_m(i, j) * m_v;
}
private:
const M &m_m;
double m_v;
};
template <typename M>
MatrixMultiplication<M> operator*(const MatrixExpression<M> &m, double v)
{
return MatrixMultiplication(*static_cast<const M *>(&m), v);
}
// --------------------------------------------------------------------
template <class M1>
Matrix Cofactors(const M1 &m)
{
Matrix cf(m.dim_m(), m.dim_m());
const size_t ixs[4][3] = {
{ 1, 2, 3 },
{ 0, 2, 3 },
{ 0, 1, 3 },
{ 0, 1, 2 }
};
for (size_t x = 0; x < 4; ++x)
{
const size_t *ix = ixs[x];
for (size_t y = 0; y < 4; ++y)
{
const size_t *iy = ixs[y];
cf(x, y) =
m(ix[0], iy[0]) * m(ix[1], iy[1]) * m(ix[2], iy[2]) +
m(ix[0], iy[1]) * m(ix[1], iy[2]) * m(ix[2], iy[0]) +
m(ix[0], iy[2]) * m(ix[1], iy[0]) * m(ix[2], iy[1]) -
m(ix[0], iy[2]) * m(ix[1], iy[1]) * m(ix[2], iy[0]) -
m(ix[0], iy[1]) * m(ix[1], iy[0]) * m(ix[2], iy[2]) -
m(ix[0], iy[0]) * m(ix[1], iy[2]) * m(ix[2], iy[1]);
if ((x + y) % 2 == 1)
cf(x, y) *= -1;
}
}
return cf;
}
// --------------------------------------------------------------------
template<typename T>
quaternion_type<T> normalize(quaternion_type<T> q)
......@@ -443,8 +205,8 @@ double LargestDepressedQuarticSolution(double a, double b, double c)
quaternion align_points(const std::vector<point> &pa, const std::vector<point> &pb)
{
// First calculate M, a 3x3 Matrix containing the sums of products of the coordinates of A and B
Matrix M(3, 3, 0);
// First calculate M, a 3x3 matrix containing the sums of products of the coordinates of A and B
matrix3x3<double> M;
for (uint32_t i = 0; i < pa.size(); ++i)
{
......@@ -462,8 +224,8 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
M(2, 2) += a.m_z * b.m_z;
}
// Now calculate N, a symmetric 4x4 Matrix
SymmetricMatrix N(4);
// Now calculate N, a symmetric 4x4 matrix
symmetric_matrix4x4<double> N(4);
N(0, 0) = M(0, 0) + M(1, 1) + M(2, 2);
N(0, 1) = M(1, 2) - M(2, 1);
......@@ -512,10 +274,10 @@ quaternion align_points(const std::vector<point> &pa, const std::vector<point> &
double lambda = LargestDepressedQuarticSolution(C, D, E);
// calculate t = (N - λI)
Matrix t = N - IdentityMatrix(4) * lambda;
matrix<double> t(N - identity_matrix(4) * lambda);
// calculate a Matrix of cofactors for t
Matrix cf = Cofactors(t);
// calculate a matrix of cofactors for t
auto cf = matrix_cofactors(t);
int maxR = 0;
for (int r = 1; r < 4; ++r)
......
src/symmetry.cpp
/*-
* SPDX-License-Identifier: BSD-2-Clause
*
*
* Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
*
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
*
* 1. Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
*
* 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
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
......@@ -25,6 +25,8 @@
*/
#include "cif++/symmetry.hpp"
#include "cif++/datablock.hpp"
#include "cif++/point.hpp"
#include <stdexcept>
......@@ -34,13 +36,140 @@ namespace cif
{
// --------------------------------------------------------------------
// Unfortunately, clipper has a different numbering scheme than CCP4
cell::cell(float a, float b, float c, float alpha, float beta, float gamma)
: m_a(a)
, m_b(b)
, m_c(c)
, m_alpha(alpha)
, m_beta(beta)
, m_gamma(gamma)
{
init();
}
cell::cell(const datablock &db)
{
auto &_cell = db["cell"];
cif::tie(m_a, m_b, m_c, m_alpha, m_beta, m_gamma) =
_cell.front().get("length_a", "length_b", "length_c", "angle_alpha", "angle_beta", "angle_gamma");
init();
}
void cell::init()
{
auto alpha = (m_alpha * kPI) / 180;
auto beta = (m_beta * kPI) / 180;
auto gamma = (m_gamma * kPI) / 180;
auto alpha_star = std::acos((std::cos(gamma) * std::cos(beta) - std::cos(alpha)) / (std::sin(beta) * std::sin(gamma)));
m_orthogonal = identity_matrix(3);
m_orthogonal(0, 0) = m_a;
m_orthogonal(0, 1) = m_b * std::cos(gamma);
m_orthogonal(0, 2) = m_c * std::cos(beta);
m_orthogonal(1, 1) = m_b * std::sin(gamma);
m_orthogonal(1, 2) = -m_c * std::sin(beta) * std::cos(alpha_star);
m_orthogonal(2, 2) = m_c * std::sin(beta) * std::sin(alpha_star);
m_fractional = inverse(m_orthogonal);
}
// --------------------------------------------------------------------
sym_op::sym_op(std::string_view s)
{
auto b = s.data();
auto e = b + s.length();
int rnri;
auto r = std::from_chars(b, e, rnri);
m_nr = rnri;
m_ta = r.ptr[1] - '0';
m_tb = r.ptr[2] - '0';
m_tc = r.ptr[3] - '0';
if (r.ec != std::errc() or rnri > 192 or r.ptr[0] != '_' or m_ta > 9 or m_tb > 9 or m_tc > 9)
throw std::invalid_argument("Could not convert string into sym_op");
}
// --------------------------------------------------------------------
transformation::transformation(const symop_data &data)
{
const auto &d = data.data();
m_rotation(0, 0) = d[0];
m_rotation(0, 1) = d[1];
m_rotation(0, 2) = d[2];
m_rotation(1, 0) = d[3];
m_rotation(1, 1) = d[4];
m_rotation(1, 2) = d[5];
m_rotation(2, 0) = d[6];
m_rotation(2, 1) = d[7];
m_rotation(2, 2) = d[8];
m_translation.m_x = d[9] == 0 ? 0 : 1.0 * d[9] / d[10];
m_translation.m_y = d[11] == 0 ? 0 : 1.0 * d[11] / d[12];
m_translation.m_y = d[13] == 0 ? 0 : 1.0 * d[13] / d[14];
}
// --------------------------------------------------------------------
spacegroup::spacegroup(int nr)
: m_nr(nr)
{
const size_t N = cif::kSymopNrTableSize;
int32_t L = 0, R = static_cast<int32_t>(N - 1);
while (L <= R)
{
int32_t i = (L + R) / 2;
if (cif::kSymopNrTable[i].spacegroup() < m_nr)
L = i + 1;
else
R = i - 1;
}
m_index = L;
for (size_t i = L; i < N and cif::kSymopNrTable[i].spacegroup() == m_nr; ++i)
emplace_back(cif::kSymopNrTable[i].symop().data());
}
std::string spacegroup::get_name() const
{
for (auto &s : kSpaceGroups)
{
if (s.nr == m_nr)
return s.name;
}
throw std::runtime_error("Spacegroup has an invalid number: " + std::to_string(m_nr));
}
point spacegroup::operator()(const point &pt, const cell &c, sym_op symop)
{
if (symop.m_nr < 1 or symop.m_nr > size())
throw std::out_of_range("symmetry operator number out of range");
transformation t = at(symop.m_nr - 1);
return pt;
}
// --------------------------------------------------------------------
// Unfortunately, clipper has a different numbering scheme than PDB
// for rotation numbers. So we created a table to map those.
// Perhaps a bit over the top, but hey....
// --------------------------------------------------------------------
int get_space_group_number(std::string spacegroup)
int get_space_group_number(std::string_view spacegroup)
{
if (spacegroup == "P 21 21 2 A")
spacegroup = "P 21 21 2 (a)";
......@@ -73,7 +202,7 @@ int get_space_group_number(std::string spacegroup)
{
for (size_t i = 0; i < kNrOfSpaceGroups; ++i)
{
auto& sp = kSpaceGroups[i];
auto &sp = kSpaceGroups[i];
if (sp.xHM == spacegroup)
{
result = sp.nr;
......@@ -83,14 +212,14 @@ int get_space_group_number(std::string spacegroup)
}
if (result == 0)
throw std::runtime_error("Spacegroup name " + spacegroup + " was not found in table");
throw std::runtime_error("Spacegroup name " + std::string(spacegroup) + " was not found in table");
return result;
}
// --------------------------------------------------------------------
int get_space_group_number(std::string spacegroup, space_group_name type)
int get_space_group_number(std::string_view spacegroup, space_group_name type)
{
if (spacegroup == "P 21 21 2 A")
spacegroup = "P 21 21 2 (a)";
......@@ -145,9 +274,9 @@ int get_space_group_number(std::string spacegroup, space_group_name type)
// not found, see if we can find a match based on xHM name
if (result == 0)
throw std::runtime_error("Spacegroup name " + spacegroup + " was not found in table");
throw std::runtime_error("Spacegroup name " + std::string(spacegroup) + " was not found in table");
return result;
}
}
} // namespace cif
src/symop-map-generator.cpp
......@@ -315,7 +315,10 @@ int main(int argc, char* const argv[])
{
case State::skip:
if (line == "begin_spacegroup")
{
state = State::spacegroup;
cur = {};
}
break;
case State::spacegroup:
......
src/symop_table_data.hpp
......@@ -9,6 +9,7 @@ namespace cif
const space_group kSpaceGroups[] =
{
{ "" , "P 2 1 1" , " P 2y (y,z,x)" , 0 },
{ "A 1 2 1" , "A 1 2 1" , " C 2y (z,y,-x)" , 2005 },
{ "A 2" , "A 1 2 1" , " C 2y (z,y,-x)" , 2005 },
{ "A b a 2" , "A b a 2" , " A 2 -2ab" , 41 },
......@@ -20,63 +21,7 @@ const space_group kSpaceGroups[] =
{ "B 1 1 2/m" , "B 1 1 2/m" , "-C 2y (-x,z,y)" , 1012 },
{ "B 1 1 b" , "B 1 1 b" , " C -2yc (-x,z,y)" , 1009 },
{ "B 1 1 m" , "B 1 1 m" , " B -2" , 1008 },
{ "B 2" , "P 1 m 1" , " P -2y" , 6 },
{ "B 2" , "P 1 c 1" , " P -2yc" , 7 },
{ "B 2" , "C 1 m 1" , " C -2y" , 8 },
{ "B 2" , "C 1 c 1" , " C -2yc" , 9 },
{ "B 2" , "P 1 2/m 1" , "-P 2y" , 10 },
{ "B 2" , "P 1 21/m 1" , "-P 2yb" , 11 },
{ "B 2" , "C 1 2/m 1" , "-C 2y" , 12 },
{ "B 2" , "P 1 2/c 1" , "-P 2yc" , 13 },
{ "B 2" , "P 1 21/c 1" , "-P 2ybc" , 14 },
{ "B 2" , "C 1 2/c 1" , "-C 2yc" , 15 },
{ "B 2" , "P 2 2 2" , " P 2 2" , 16 },
{ "B 2" , "P 2 2 21" , " P 2c 2" , 17 },
{ "B 2" , "P 21 21 2" , " P 2 2ab" , 18 },
{ "B 2" , "P 21 21 21" , " P 2ac 2ab" , 19 },
{ "B 2" , "C 2 2 21" , " C 2c 2" , 20 },
{ "B 2" , "C 2 2 2" , " C 2 2" , 21 },
{ "B 2" , "F 2 2 2" , " F 2 2" , 22 },
{ "B 2" , "I 2 2 2" , " I 2 2" , 23 },
{ "B 2" , "I 21 21 21" , " I 2b 2c" , 24 },
{ "B 2" , "P m m 2" , " P 2 -2" , 25 },
{ "B 2" , "P m c 21" , " P 2c -2" , 26 },
{ "B 2" , "P c c 2" , " P 2 -2c" , 27 },
{ "B 2" , "P m a 2" , " P 2 -2a" , 28 },
{ "B 2" , "P c a 21" , " P 2c -2ac" , 29 },
{ "B 2" , "P n c 2" , " P 2 -2bc" , 30 },
{ "B 2" , "P m n 21" , " P 2ac -2" , 31 },
{ "B 2" , "P b a 2" , " P 2 -2ab" , 32 },
{ "B 2" , "P n a 21" , " P 2c -2n" , 33 },
{ "B 2" , "P n n 2" , " P 2 -2n" , 34 },
{ "B 2" , "C m m 2" , " C 2 -2" , 35 },
{ "B 2" , "C m c 21" , " C 2c -2" , 36 },
{ "B 2" , "C c c 2" , " C 2 -2c" , 37 },
{ "B 2" , "A m m 2" , " A 2 -2" , 38 },
{ "B 2" , "A b m 2" , " A 2 -2b" , 39 },
{ "B 2" , "A m a 2" , " A 2 -2a" , 40 },
{ "B 2" , "A b a 2" , " A 2 -2ab" , 41 },
{ "B 2" , "F m m 2" , " F 2 -2" , 42 },
{ "B 2" , "F d d 2" , " F 2 -2d" , 43 },
{ "B 2" , "I m m 2" , " I 2 -2" , 44 },
{ "B 2" , "I b a 2" , " I 2 -2c" , 45 },
{ "B 2" , "I m a 2" , " I 2 -2a" , 46 },
{ "B 2" , "B 1 1 2" , " C 2y (-x,z,y)" , 1005 },
{ "B 2" , "P 1 1 m" , " P -2y (z,x,y)" , 1006 },
{ "B 2" , "P 1 1 b" , " P -2yc (-x,z,y)" , 1007 },
{ "B 2" , "B 1 1 b" , " C -2yc (-x,z,y)" , 1009 },
{ "B 2" , "P 1 1 2/m" , "-P 2y (z,x,y)" , 1010 },
{ "B 2" , "P 1 1 21/m" , "-P 2yb (z,x,y)" , 1011 },
{ "B 2" , "B 1 1 2/m" , "-C 2y (-x,z,y)" , 1012 },
{ "B 2" , "P 1 1 2/b" , "-P 2yc (-x,z,y)" , 1013 },
{ "B 2" , "P 1 1 21/b" , "-P 2ybc (-x,z,y)" , 1014 },
{ "B 2" , "B 1 1 2/b" , "-C 2yc (-x,z,y)" , 1015 },
{ "B 2" , "P 21 2 2" , " P 2c 2 (z,x,y)" , 1017 },
{ "B 2" , "P 1 21/n 1" , "-P 2ybc (x-z,y,z)" , 2014 },
{ "B 2" , "P 2 21 2" , " P 2c 2 (y,z,x)" , 2017 },
{ "B 2" , "P 21 2 21" , " P 2 2ab (y,z,x)" , 2018 },
{ "B 2" , "P 1 21/a 1" , "-P 2ybc (z,y,-x)" , 3014 },
{ "B 2" , "P 2 21 21" , " P 2 2ab (z,x,y)" , 3018 },
{ "C 1 2 1" , "C 1 2 1" , " C 2y" , 5 },
{ "C 1 2/c 1" , "C 1 2/c 1" , "-C 2yc" , 15 },
{ "C 1 2/m 1" , "C 1 2/m 1" , "-C 2y" , 12 },
......@@ -85,8 +30,8 @@ const space_group kSpaceGroups[] =
{ "C 1 m 1" , "C 1 m 1" , " C -2y" , 8 },
{ "C 2 2 2" , "C 2 2 2" , " C 2 2" , 21 },
{ "C 2 2 21" , "C 2 2 21" , " C 2c 2" , 20 },
{ "C 2 2 21a)" , "B 1 1 m" , " C 2c 2 (x+1/4,y,z)" , 1020 },
{ "C 2 2 2a" , "B 1 1 m" , " C 2 2 (x+1/4,y+1/4,z)" , 1021 },
{ "C 2 2 21a)" , "" , " C 2c 2 (x+1/4,y,z)" , 1020 },
{ "C 2 2 2a" , "" , " C 2 2 (x+1/4,y+1/4,z)" , 1021 },
{ "C 2/c 2/c 2/a" , "C c c a :1" , "-C 2a 2ac (x-1/2,y-1/4,z+1/4)" , 68 },
{ "C 2/c 2/c 2/m" , "C c c m" , "-C 2 2c" , 66 },
{ "C 2/m 2/c 21/a" , "C m c a" , "-C 2ac 2" , 64 },
......@@ -105,7 +50,7 @@ const space_group kSpaceGroups[] =
{ "F -4 3 c" , "F -4 3 c" , " F -4a 2 3" , 219 },
{ "F -4 3 m" , "F -4 3 m" , " F -4 2 3" , 216 },
{ "F 2 2 2" , "F 2 2 2" , " F 2 2" , 22 },
{ "F 2 2 2a" , "B 1 1 m" , " F 2 2 (x,y,z+1/4)" , 1022 },
{ "F 2 2 2a" , "" , " F 2 2 (x,y,z+1/4)" , 1022 },
{ "F 2 3" , "F 2 3" , " F 2 2 3" , 196 },
{ "F 2/d -3" , "F d -3 :1" , "-F 2uv 2vw 3 (x+1/8,y+1/8,z+1/8)" , 203 },
{ "F 2/d 2/d 2/d" , "F d d d :1" , "-F 2uv 2vw (x+1/8,y+1/8,z+1/8)" , 70 },
......@@ -147,9 +92,9 @@ const space_group kSpaceGroups[] =
{ "I 1 21 1" , "I 1 21 1" , " C 2y (x+1/4,y+1/4,-x+z-1/4)" , 5005 },
{ "I 2" , "I 1 2 1" , " C 2y (x,y,-x+z)" , 4005 },
{ "I 2 2 2" , "I 2 2 2" , " I 2 2" , 23 },
{ "I 2 2 2a" , "B 1 1 m" , " I 2 2 (x-1/4,y+1/4,z-1/4)" , 1023 },
{ "I 2 2 2a" , "" , " I 2 2 (x-1/4,y+1/4,z-1/4)" , 1023 },
{ "I 2 3" , "I 2 3" , " I 2 2 3" , 197 },
{ "I 2 3a" , "B 1 1 m" , " I 2 2 3 (x+1/4,y+1/4,z+1/4)" , 1197 },
{ "I 2 3a" , "" , " I 2 2 3 (x+1/4,y+1/4,z+1/4)" , 1197 },
{ "I 2/b 2/a 2/m" , "I b a m" , "-I 2 2c" , 72 },
{ "I 2/m -3" , "I m -3" , "-I 2 2 3" , 204 },
{ "I 2/m 2/m 2/m" , "I m m m" , "-I 2 2" , 71 },
......@@ -175,57 +120,9 @@ const space_group kSpaceGroups[] =
{ "I 41/a" , "I 41/a :1" , "-I 4ad (x,y+1/4,z+1/8)" , 88 },
{ "I 41/a -3 2/d" , "I a -3 d" , "-I 4bd 2c 3" , 230 },
{ "I 41/a 2/c 2/d" , "I 41/a c d :1" , "-I 4bd 2c (x-1/2,y+1/4,z-3/8)" , 142 },
{ "I 41/a 2/c 2/d" , "P 3" , " P 3" , 143 },
{ "I 41/a 2/c 2/d" , "P 31" , " P 31" , 144 },
{ "I 41/a 2/c 2/d" , "P 32" , " P 32" , 145 },
{ "I 41/a 2/c 2/d" , "R 3 :H" , " R 3" , 146 },
{ "I 41/a 2/c 2/d" , "P -3" , "-P 3" , 147 },
{ "I 41/a 2/c 2/d" , "R -3 :H" , "-R 3" , 148 },
{ "I 41/a 2/c 2/d" , "P 3 1 2" , " P 3 2" , 149 },
{ "I 41/a 2/c 2/d" , "P 3 2 1" , " P 3 2\"" , 150 },
{ "I 41/a 2/c 2/d" , "P 31 1 2" , " P 31 2 (x,y,z+1/3)" , 151 },
{ "I 41/a 2/c 2/d" , "P 31 2 1" , " P 31 2\"" , 152 },
{ "I 41/a 2/c 2/d" , "P 32 1 2" , " P 32 2 (x,y,z+1/6)" , 153 },
{ "I 41/a 2/c 2/d" , "P 32 2 1" , " P 32 2\"" , 154 },
{ "I 41/a 2/c 2/d" , "R 3 2 :H" , " R 3 2\"" , 155 },
{ "I 41/a 2/c 2/d" , "P 3 m 1" , " P 3 -2\"" , 156 },
{ "I 41/a 2/c 2/d" , "P 3 1 m" , " P 3 -2" , 157 },
{ "I 41/a 2/c 2/d" , "P 3 c 1" , " P 3 -2\"c" , 158 },
{ "I 41/a 2/c 2/d" , "P 3 1 c" , " P 3 -2c" , 159 },
{ "I 41/a 2/c 2/d" , "R 3 m :H" , " R 3 -2\"" , 160 },
{ "I 41/a 2/c 2/d" , "R 3 c :H" , " R 3 -2\"c" , 161 },
{ "I 41/a 2/c 2/d" , "R 3 :R" , " R 3 (-y+z,x+z,-x+y+z)" , 1146 },
{ "I 41/a 2/c 2/d" , "R -3 :R" , "-R 3 (-y+z,x+z,-x+y+z)" , 1148 },
{ "I 41/a 2/c 2/d" , "R 3 2 :R" , " R 3 2\" (-y+z,x+z,-x+y+z)" , 1155 },
{ "I 41/a 2/c 2/d" , "R 3 m :R" , " R 3 -2\" (-y+z,x+z,-x+y+z)" , 1160 },
{ "I 41/a 2/c 2/d" , "R 3 c :R" , " R 3 -2\"c (-y+z,x+z,-x+y+z)" , 1161 },
{ "I 41/a 2/m 2/d" , "I 41/a m d :1" , "-I 4bd 2 (x-1/2,y+1/4,z+1/8)" , 141 },
{ "I a -3" , "I a -3" , "-I 2b 2c 3" , 206 },
{ "I a -3" , "P 4 3 2" , " P 4 2 3" , 207 },
{ "I a -3" , "P 42 3 2" , " P 4n 2 3" , 208 },
{ "I a -3" , "F 4 3 2" , " F 4 2 3" , 209 },
{ "I a -3" , "F 41 3 2" , " F 4d 2 3" , 210 },
{ "I a -3" , "I 4 3 2" , " I 4 2 3" , 211 },
{ "I a -3" , "P 43 3 2" , " P 4acd 2ab 3" , 212 },
{ "I a -3" , "P 41 3 2" , " P 4bd 2ab 3" , 213 },
{ "I a -3" , "I 41 3 2" , " I 4bd 2c 3" , 214 },
{ "I a -3" , "P -4 3 m" , " P -4 2 3" , 215 },
{ "I a -3" , "F -4 3 m" , " F -4 2 3" , 216 },
{ "I a -3" , "I -4 3 m" , " I -4 2 3" , 217 },
{ "I a -3" , "P -4 3 n" , " P -4n 2 3" , 218 },
{ "I a -3" , "F -4 3 c" , " F -4a 2 3" , 219 },
{ "I a -3" , "I -4 3 d" , " I -4bd 2c 3" , 220 },
{ "I a -3 d" , "I a -3 d" , "-I 4bd 2c 3" , 230 },
{ "I a -3 d" , "B 1 1 m" , " B -2" , 1008 },
{ "I a -3 d" , "B 1 1 m" , " P 2 2ab (x+1/4,y+1/4,z)" , 1018 },
{ "I a -3 d" , "B 1 1 m" , " C 2c 2 (x+1/4,y,z)" , 1020 },
{ "I a -3 d" , "B 1 1 m" , " C 2 2 (x+1/4,y+1/4,z)" , 1021 },
{ "I a -3 d" , "B 1 1 m" , " F 2 2 (x,y,z+1/4)" , 1022 },
{ "I a -3 d" , "B 1 1 m" , " I 2 2 (x-1/4,y+1/4,z-1/4)" , 1023 },
{ "I a -3 d" , "B 1 1 m" , " P 4n 2n (x-1/4,y-1/4,z-1/4)" , 1094 },
{ "I a -3 d" , "B 1 1 m" , " I 2 2 3 (x+1/4,y+1/4,z+1/4)" , 1197 },
{ "I a -3 d" , "C 1 21 1" , " C 2y (x+1/4,y+1/4,z)" , 3005 },
{ "I a -3 d" , "I 1 21 1" , " C 2y (x+1/4,y+1/4,-x+z-1/4)" , 5005 },
{ "I b a 2" , "I b a 2" , " I 2 -2c" , 45 },
{ "I b a m" , "I b a m" , "-I 2 2c" , 72 },
{ "I b c a" , "I b c a" , "-I 2b 2c" , 73 },
......@@ -234,54 +131,6 @@ const space_group kSpaceGroups[] =
{ "I m a 2" , "I m a 2" , " I 2 -2a" , 46 },
{ "I m m 2" , "I m m 2" , " I 2 -2" , 44 },
{ "I m m a" , "I m m a" , "-I 2b 2" , 74 },
{ "I m m a" , "P 4" , " P 4" , 75 },
{ "I m m a" , "P 41" , " P 4w" , 76 },
{ "I m m a" , "P 42" , " P 4c" , 77 },
{ "I m m a" , "P 43" , " P 4cw" , 78 },
{ "I m m a" , "I 4" , " I 4" , 79 },
{ "I m m a" , "I 41" , " I 4bw" , 80 },
{ "I m m a" , "P -4" , " P -4" , 81 },
{ "I m m a" , "I -4" , " I -4" , 82 },
{ "I m m a" , "P 4/m" , "-P 4" , 83 },
{ "I m m a" , "P 42/m" , "-P 4c" , 84 },
{ "I m m a" , "P 4/n :1" , "-P 4a (x-1/4,y+1/4,z)" , 85 },
{ "I m m a" , "P 42/n :1" , "-P 4bc (x+1/4,y+1/4,z+1/4)" , 86 },
{ "I m m a" , "I 4/m" , "-I 4" , 87 },
{ "I m m a" , "I 41/a :1" , "-I 4ad (x,y+1/4,z+1/8)" , 88 },
{ "I m m a" , "P 4 2 2" , " P 4 2" , 89 },
{ "I m m a" , "P 4 21 2" , " P 4ab 2ab" , 90 },
{ "I m m a" , "P 41 2 2" , " P 4w 2c" , 91 },
{ "I m m a" , "P 41 21 2" , " P 4abw 2nw" , 92 },
{ "I m m a" , "P 42 2 2" , " P 4c 2" , 93 },
{ "I m m a" , "P 42 21 2" , " P 4n 2n" , 94 },
{ "I m m a" , "P 43 2 2" , " P 4cw 2c" , 95 },
{ "I m m a" , "P 43 21 2" , " P 4nw 2abw" , 96 },
{ "I m m a" , "I 4 2 2" , " I 4 2" , 97 },
{ "I m m a" , "I 41 2 2" , " I 4bw 2bw" , 98 },
{ "I m m a" , "P 4 m m" , " P 4 -2" , 99 },
{ "I m m a" , "P 4 b m" , " P 4 -2ab" , 100 },
{ "I m m a" , "P 42 c m" , " P 4c -2c" , 101 },
{ "I m m a" , "P 42 n m" , " P 4n -2n" , 102 },
{ "I m m a" , "P 4 c c" , " P 4 -2c" , 103 },
{ "I m m a" , "P 4 n c" , " P 4 -2n" , 104 },
{ "I m m a" , "P 42 m c" , " P 4c -2" , 105 },
{ "I m m a" , "P 42 b c" , " P 4c -2ab" , 106 },
{ "I m m a" , "I 4 m m" , " I 4 -2" , 107 },
{ "I m m a" , "I 4 c m" , " I 4 -2c" , 108 },
{ "I m m a" , "I 41 m d" , " I 4bw -2" , 109 },
{ "I m m a" , "I 41 c d" , " I 4bw -2c" , 110 },
{ "I m m a" , "P -4 2 m" , " P -4 2" , 111 },
{ "I m m a" , "P -4 2 c" , " P -4 2c" , 112 },
{ "I m m a" , "P -4 21 m" , " P -4 2ab" , 113 },
{ "I m m a" , "P -4 21 c" , " P -4 2n" , 114 },
{ "I m m a" , "P -4 m 2" , " P -4 -2" , 115 },
{ "I m m a" , "P -4 c 2" , " P -4 -2c" , 116 },
{ "I m m a" , "P -4 b 2" , " P -4 -2ab" , 117 },
{ "I m m a" , "P -4 n 2" , " P -4 -2n" , 118 },
{ "I m m a" , "I -4 m 2" , " I -4 -2" , 119 },
{ "I m m a" , "I -4 c 2" , " I -4 -2c" , 120 },
{ "I m m a" , "I -4 2 m" , " I -4 2" , 121 },
{ "I m m a" , "I -4 2 d" , " I -4 2bw" , 122 },
{ "I m m m" , "I m m m" , "-I 2 2" , 71 },
{ "I4/m c m" , "I 4/m c m" , "-I 4 2c" , 140 },
{ "I4/m m m" , "I 4/m m m" , "-I 4 2" , 139 },
......@@ -349,7 +198,7 @@ const space_group kSpaceGroups[] =
{ "P 21 2 2" , "P 21 2 2" , " P 2c 2 (z,x,y)" , 1017 },
{ "P 21 2 21" , "P 21 2 21" , " P 2 2ab (y,z,x)" , 2018 },
{ "P 21 21 2" , "P 21 21 2" , " P 2 2ab" , 18 },
{ "P 21 21 2 (a)" , "B 1 1 m" , " P 2 2ab (x+1/4,y+1/4,z)" , 1018 },
{ "P 21 21 2 (a)" , "" , " P 2 2ab (x+1/4,y+1/4,z)" , 1018 },
{ "P 21 21 21" , "P 21 21 21" , " P 2ac 2ab" , 19 },
{ "P 21 3" , "P 21 3" , " P 2ac 2ab 3" , 198 },
{ "P 21/a -3" , "P a -3" , "-P 2ac 2ab 3" , 205 },
......@@ -403,7 +252,7 @@ const space_group kSpaceGroups[] =
{ "P 42" , "P 42" , " P 4c" , 77 },
{ "P 42 2 2" , "P 42 2 2" , " P 4c 2" , 93 },
{ "P 42 21 2" , "P 42 21 2" , " P 4n 2n" , 94 },
{ "P 42 21 2a" , "B 1 1 m" , " P 4n 2n (x-1/4,y-1/4,z-1/4)" , 1094 },
{ "P 42 21 2a" , "" , " P 4n 2n (x-1/4,y-1/4,z-1/4)" , 1094 },
{ "P 42 3 2" , "P 42 3 2" , " P 4n 2 3" , 208 },
{ "P 42 b c" , "P 42 b c" , " P 4c -2ab" , 106 },
{ "P 42 c m" , "P 42 c m" , " P 4c -2c" , 101 },
......@@ -447,11 +296,6 @@ const space_group kSpaceGroups[] =
{ "P 63/m 2/m 2/c" , "P 63/m m c" , "-P 6c 2c" , 194 },
{ "P 63/m c m" , "P 63/m c m" , "-P 6c 2" , 193 },
{ "P 63/m m c" , "P 63/m m c" , "-P 6c 2c" , 194 },
{ "P 63/m m c" , "P 2 3" , " P 2 2 3" , 195 },
{ "P 63/m m c" , "F 2 3" , " F 2 2 3" , 196 },
{ "P 63/m m c" , "I 2 3" , " I 2 2 3" , 197 },
{ "P 63/m m c" , "P 21 3" , " P 2ac 2ab 3" , 198 },
{ "P 63/m m c" , "I 21 3" , " I 2b 2c 3" , 199 },
{ "P 64" , "P 64" , " P 64" , 172 },
{ "P 64 2 2" , "P 64 2 2" , " P 64 2 (x,y,z+1/6)" , 181 },
{ "P 65" , "P 65" , " P 65" , 170 },
......@@ -477,7 +321,6 @@ const space_group kSpaceGroups[] =
{ "P m m a" , "P m m a" , "-P 2a 2a" , 51 },
{ "P m m m" , "P m m m" , "-P 2 2" , 47 },
{ "P m m n" , "P m m n :1" , "-P 2ab 2a (x-1/4,y-1/4,z)" , 59 },
{ "P m m n" , "P m m n :2" , "-P 2ab 2a" , 1059 },
{ "P m n 21" , "P m n 21" , " P 2ac -2" , 31 },
{ "P m n a" , "P m n a" , "-P 2ac 2" , 53 },
{ "P n -3" , "P n -3 :1" , "-P 2ab 2bc 3 (x-1/4,y-1/4,z-1/4)" , 201 },
......@@ -509,29 +352,6 @@ const space_group kSpaceGroups[] =
{ "R -3" , "R -3 :R" , "-R 3 (-y+z,x+z,-x+y+z)" , 1148 },
{ "R -3 2/c" , "R -3 c :R" , "-R 3 2\"c (-y+z,x+z,-x+y+z)" , 1167 },
{ "R -3 2/m" , "R -3 m :R" , "-R 3 2\" (-y+z,x+z,-x+y+z)" , 1166 },
{ "R -3 c" , "P 6" , " P 6" , 168 },
{ "R -3 c" , "P 61" , " P 61" , 169 },
{ "R -3 c" , "P 65" , " P 65" , 170 },
{ "R -3 c" , "P 62" , " P 62" , 171 },
{ "R -3 c" , "P 64" , " P 64" , 172 },
{ "R -3 c" , "P 63" , " P 6c" , 173 },
{ "R -3 c" , "P -6" , " P -6" , 174 },
{ "R -3 c" , "P 6/m" , "-P 6" , 175 },
{ "R -3 c" , "P 63/m" , "-P 6c" , 176 },
{ "R -3 c" , "P 6 2 2" , " P 6 2" , 177 },
{ "R -3 c" , "P 61 2 2" , " P 61 2 (x,y,z+5/12)" , 178 },
{ "R -3 c" , "P 65 2 2" , " P 65 2 (x,y,z+1/12)" , 179 },
{ "R -3 c" , "P 62 2 2" , " P 62 2 (x,y,z+1/3)" , 180 },
{ "R -3 c" , "P 64 2 2" , " P 64 2 (x,y,z+1/6)" , 181 },
{ "R -3 c" , "P 63 2 2" , " P 6c 2c" , 182 },
{ "R -3 c" , "P 6 m m" , " P 6 -2" , 183 },
{ "R -3 c" , "P 6 c c" , " P 6 -2c" , 184 },
{ "R -3 c" , "P 63 c m" , " P 6c -2" , 185 },
{ "R -3 c" , "P 63 m c" , " P 6c -2c" , 186 },
{ "R -3 c" , "P -6 m 2" , " P -6 2" , 187 },
{ "R -3 c" , "P -6 c 2" , " P -6c 2" , 188 },
{ "R -3 c" , "P -6 2 m" , " P -6 -2" , 189 },
{ "R -3 c" , "P -6 2 c" , " P -6c -2c" , 190 },
{ "R -3 c" , "R -3 c :R" , "-R 3 2\"c (-y+z,x+z,-x+y+z)" , 1167 },
{ "R -3 m" , "R -3 m :R" , "-R 3 2\" (-y+z,x+z,-x+y+z)" , 1166 },
{ "R 3" , "R 3 :R" , " R 3 (-y+z,x+z,-x+y+z)" , 1146 },
......@@ -5270,12 +5090,12 @@ const symop_datablock kSymopNrTable[] = {
{ 1017, 2, { -1, 0, 0, 0, 1, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, } },
{ 1017, 3, { 1, 0, 0, 0,-1, 0, 0, 0,-1, 1, 2, 0, 0, 0, 0, } },
{ 1017, 4, { -1, 0, 0, 0,-1, 0, 0, 0, 1, 1, 2, 0, 0, 0, 0, } },
// B 1 1 m
//
{ 1018, 1, { 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, } },
{ 1018, 2, { -1, 0, 0, 0,-1, 0, 0, 0, 1, 1, 2, 1, 2, 0, 0, } },
{ 1018, 3, { 1, 0, 0, 0,-1, 0, 0, 0,-1, 1, 2, 0, 0, 0, 0, } },
{ 1018, 4, { -1, 0, 0, 0, 1, 0, 0, 0,-1, 0, 0, 1, 2, 0, 0, } },
// B 1 1 m
//
{ 1020, 1, { 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, } },
{ 1020, 2, { -1, 0, 0, 0,-1, 0, 0, 0, 1, 1, 2, 0, 0, 1, 2, } },
{ 1020, 3, { 1, 0, 0, 0,-1, 0, 0, 0,-1, 1, 2, 1, 2, 0, 0, } },
......@@ -5284,7 +5104,7 @@ const symop_datablock kSymopNrTable[] = {
{ 1020, 6, { -1, 0, 0, 0,-1, 0, 0, 0, 1, 0, 0, 1, 2, 1, 2, } },
{ 1020, 7, { 1, 0, 0, 0,-1, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, } },
{ 1020, 8, { -1, 0, 0, 0, 1, 0, 0, 0,-1, 1, 2, 0, 0, 1, 2, } },
// B 1 1 m
//
{ 1021, 1, { 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, } },
{ 1021, 2, { -1, 0, 0, 0,-1, 0, 0, 0, 1, 1, 2, 1, 2, 0, 0, } },
{ 1021, 3, { 1, 0, 0, 0,-1, 0, 0, 0,-1, 1, 2, 0, 0, 0, 0, } },
......@@ -5293,7 +5113,7 @@ const symop_datablock kSymopNrTable[] = {
{ 1021, 6, { -1, 0, 0, 0,-1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, } },
{ 1021, 7, { 1, 0, 0, 0,-1, 0, 0, 0,-1, 0, 0, 1, 2, 0, 0, } },
{ 1021, 8, { -1, 0, 0, 0, 1, 0, 0, 0,-1, 1, 2, 0, 0, 0, 0, } },
// B 1 1 m
//
{ 1022, 1, { 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, } },
{ 1022, 2, { -1, 0, 0, 0,-1, 0, 0, 0, 1, 1, 2, 1, 2, 0, 0, } },
{ 1022, 3, { 1, 0, 0, 0,-1, 0, 0, 0,-1, 1, 2, 0, 0, 0, 0, } },
......@@ -5310,7 +5130,7 @@ const symop_datablock kSymopNrTable[] = {
{ 1022, 14, { -1, 0, 0, 0,-1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, } },
{ 1022, 15, { 1, 0, 0, 0,-1, 0, 0, 0,-1, 0, 0, 1, 2, 0, 0, } },
{ 1022, 16, { -1, 0, 0, 0, 1, 0, 0, 0,-1, 1, 2, 0, 0, 0, 0, } },
// B 1 1 m
//
{ 1023, 1, { 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, } },
{ 1023, 2, { -1, 0, 0, 0,-1, 0, 0, 0, 1, 1, 2, 1, 2, 0, 0, } },
{ 1023, 3, { 1, 0, 0, 0,-1, 0, 0, 0,-1, 1, 2, 0, 0, 0, 0, } },
......@@ -5328,7 +5148,7 @@ const symop_datablock kSymopNrTable[] = {
{ 1059, 6, { 1, 0, 0, 0, 1, 0, 0, 0,-1, 1, 2, 1, 2, 0, 0, } },
{ 1059, 7, { 1, 0, 0, 0,-1, 0, 0, 0, 1, 0, 0, 1, 2, 0, 0, } },
{ 1059, 8, { -1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 2, 0, 0, 0, 0, } },
// B 1 1 m
//
{ 1094, 1, { 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, } },
{ 1094, 2, { -1, 0, 0, 0,-1, 0, 0, 0, 1, 1, 2, 1, 2, 0, 0, } },
{ 1094, 3, { 0,-1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 2, 1, 2, } },
......@@ -5395,7 +5215,7 @@ const symop_datablock kSymopNrTable[] = {
{ 1167, 10, { 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 2, 1, 2, 1, 2, } },
{ 1167, 11, { 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 2, 1, 2, 1, 2, } },
{ 1167, 12, { 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 2, 1, 2, 1, 2, } },
// B 1 1 m
//
{ 1197, 1, { 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, } },
{ 1197, 2, { -1, 0, 0, 0,-1, 0, 0, 0, 1, 1, 2, 1, 2, 0, 0, } },
{ 1197, 3, { 1, 0, 0, 0,-1, 0, 0, 0,-1, 1, 2, 0, 0, 0, 0, } },
......
test/unit-3d-test.cpp
......@@ -246,4 +246,68 @@ BOOST_AUTO_TEST_CASE(dh_q_1)
auto dh = cif::dihedral_angle(pts[0], pts[1], pts[2], pts[3]);
BOOST_TEST(dh == angle, tt::tolerance(0.1f));
}
}
\ No newline at end of file
}
// --------------------------------------------------------------------
BOOST_AUTO_TEST_CASE(symm_1)
{
cif::cell c(10, 10, 10);
cif::point p{ 1, 1, 1 };
cif::point f = fractional(p, c);
BOOST_TEST(f.m_x == 0.1f, tt::tolerance(0.01));
BOOST_TEST(f.m_y == 0.1f, tt::tolerance(0.01));
BOOST_TEST(f.m_z == 0.1f, tt::tolerance(0.01));
cif::point o = orthogonal(f, c);
BOOST_TEST(o.m_x == 1.f, tt::tolerance(0.01));
BOOST_TEST(o.m_y == 1.f, tt::tolerance(0.01));
BOOST_TEST(o.m_z == 1.f, tt::tolerance(0.01));
}
BOOST_AUTO_TEST_CASE(symm_2)
{
using namespace cif::literals;
auto symop = "1_555"_symop;
BOOST_TEST(symop.is_identity() == true);
}
BOOST_AUTO_TEST_CASE(symm_3)
{
using namespace cif::literals;
cif::spacegroup sg(18);
BOOST_TEST(sg.size() == 4);
BOOST_TEST(sg.get_name() == "P 21 21 2");
}
BOOST_AUTO_TEST_CASE(symm_4)
{
using namespace cif::literals;
// based on 2b8h
auto sg = cif::spacegroup(154); // p 32 2 1
auto c = cif::cell(107.516, 107.516, 338.487, 90.00, 90.00, 120.00);
cif::point a{ -8.688, 79.351, 10.439 }; // O6 NAG A 500
cif::point b{ -35.356, 33.693, -3.236 }; // CG2 THR D 400
cif::point sb( -6.916, 79.34, 3.236); // 4_565 copy of b
BOOST_TEST(distance(a, sg(a, c, "1_455"_symop)) == static_cast<float>(c.get_a()));
BOOST_TEST(distance(a, sg(a, c, "1_545"_symop)) == static_cast<float>(c.get_b()));
BOOST_TEST(distance(a, sg(a, c, "1_554"_symop)) == static_cast<float>(c.get_c()));
auto sb2 = sg(b, c, "4_565"_symop);
BOOST_TEST(sb.m_x == sb2.m_x);
BOOST_TEST(sb.m_y == sb2.m_y);
BOOST_TEST(sb.m_z == sb2.m_z);
BOOST_TEST(distance(a, sb2) == 7.42f);
}
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