Gemmi C++ API
Loading...
Searching...
No Matches
twin.hpp
Go to the documentation of this file.
1// Copyright 2022 Global Phasing Ltd.
2//
3// Twinning laws.
4
5#ifndef GEMMI_TWIN_HPP_
6#define GEMMI_TWIN_HPP_
7
8#include <cstdint> // for int8_t
9#include <algorithm> // for sort
10#include <utility> // for pair
11#include "symmetry.hpp" // for Op
12#include "unitcell.hpp" // for UnitCell
13#include "cellred.hpp" // for GruberVector
14
15namespace gemmi {
16
17namespace impl {
18
19// Determination of the lattice symmetry is based on P.H. Zwart et al (2006)
20// http://legacy.ccp4.ac.uk/newsletters/newsletter44/articles/explore_metric_symmetry.html
21// which in turn is based on ideas from
22// Le Page (1982) https://doi.org/10.1107/S0021889882011959
23// Lebedev et al. (2006) https://doi.org/10.1107/S0907444905036759
24
25struct TwoFoldData {
26 std::int8_t matrix[9];
27 std::int8_t ds_axis[3]; // 2-fold axis direction in direct space
28 std::int8_t rs_axis[3]; // 2-fold axis direction in reciprocal space
29};
30
31// The template wrapper here is only to substitute C++17 inline variables
32// https://stackoverflow.com/questions/38043442/how-do-inline-variables-work
33template<class Dummy> struct TwoFold_ {
34 static const TwoFoldData table[81];
35};
36
37// Two-fold twinning operators generated by cctbx/examples/reduced_cell_two_folds.py
38// "Enumeration of the 81 2-fold symmetry operations possible for reduced cells."
39template<class Dummy>
40const TwoFoldData TwoFold_<Dummy>::table[81] = {
41{{-1, -1, -1, 0, 0, 1, 0, 1, 0}, {-1, 1, 1}, {0, 1, 1}},
42{{-1, -1, 0, 0, 1, 0, 0, -1, -1}, {1, -2, 1}, {0, 1, 0}},
43{{-1, -1, 0, 0, 1, 0, 0, 0, -1}, {-1, 2, 0}, {0, 1, 0}},
44{{-1, -1, 0, 0, 1, 0, 0, 1, -1}, {-1, 2, 1}, {0, 1, 0}},
45{{-1, -1, 1, 0, 0, -1, 0, -1, 0}, {1, -1, 1}, {0, -1, 1}},
46{{-1, 0, -1, 0, -1, -1, 0, 0, 1}, {-1, -1, 2}, {0, 0, 1}},
47{{-1, 0, -1, 0, -1, 0, 0, 0, 1}, {-1, 0, 2}, {0, 0, 1}},
48{{-1, 0, -1, 0, -1, 1, 0, 0, 1}, {-1, 1, 2}, {0, 0, 1}},
49{{-1, 0, 0, -1, 0, -1, 1, -1, 0}, {0, -1, 1}, {1, -1, 1}},
50{{-1, 0, 0, -1, 0, 1, -1, 1, 0}, {0, 1, 1}, {-1, 1, 1}},
51{{-1, 0, 0, -1, 1, -1, 0, 0, -1}, {0, 1, 0}, {1, -2, 1}},
52{{-1, 0, 0, -1, 1, 0, 0, 0, -1}, {0, 1, 0}, {-1, 2, 0}},
53{{-1, 0, 0, -1, 1, 1, 0, 0, -1}, {0, 1, 0}, {-1, 2, 1}},
54{{-1, 0, 0, 0, -1, -1, 0, 0, 1}, {0, -1, 2}, {0, 0, 1}},
55{{-1, 0, 0, 0, -1, 0, -1, -1, 1}, {0, 0, 1}, {-1, -1, 2}},
56{{-1, 0, 0, 0, -1, 0, -1, 0, 1}, {0, 0, 1}, {-1, 0, 2}},
57{{-1, 0, 0, 0, -1, 0, -1, 1, 1}, {0, 0, 1}, {-1, 1, 2}},
58{{-1, 0, 0, 0, -1, 0, 0, -1, 1}, {0, 0, 1}, {0, -1, 2}},
59{{-1, 0, 0, 0, -1, 0, 0, 0, 1}, {0, 0, 1}, {0, 0, 1}},
60{{-1, 0, 0, 0, -1, 0, 0, 1, 1}, {0, 0, 1}, {0, 1, 2}},
61{{-1, 0, 0, 0, -1, 0, 1, -1, 1}, {0, 0, 1}, {1, -1, 2}},
62{{-1, 0, 0, 0, -1, 0, 1, 0, 1}, {0, 0, 1}, {1, 0, 2}},
63{{-1, 0, 0, 0, -1, 0, 1, 1, 1}, {0, 0, 1}, {1, 1, 2}},
64{{-1, 0, 0, 0, -1, 1, 0, 0, 1}, {0, 1, 2}, {0, 0, 1}},
65{{-1, 0, 0, 0, 0, -1, 0, -1, 0}, {0, -1, 1}, {0, -1, 1}},
66{{-1, 0, 0, 0, 0, 1, 0, 1, 0}, {0, 1, 1}, {0, 1, 1}},
67{{-1, 0, 0, 0, 1, -1, 0, 0, -1}, {0, 1, 0}, {0, -2, 1}},
68{{-1, 0, 0, 0, 1, 0, 0, -1, -1}, {0, -2, 1}, {0, 1, 0}},
69{{-1, 0, 0, 0, 1, 0, 0, 0, -1}, {0, 1, 0}, {0, 1, 0}},
70{{-1, 0, 0, 0, 1, 0, 0, 1, -1}, {0, 2, 1}, {0, 1, 0}},
71{{-1, 0, 0, 0, 1, 1, 0, 0, -1}, {0, 1, 0}, {0, 2, 1}},
72{{-1, 0, 0, 1, 0, -1, -1, -1, 0}, {0, -1, 1}, {-1, -1, 1}},
73{{-1, 0, 0, 1, 0, 1, 1, 1, 0}, {0, 1, 1}, {1, 1, 1}},
74{{-1, 0, 0, 1, 1, -1, 0, 0, -1}, {0, 1, 0}, {-1, -2, 1}},
75{{-1, 0, 0, 1, 1, 0, 0, 0, -1}, {0, 1, 0}, {1, 2, 0}},
76{{-1, 0, 0, 1, 1, 1, 0, 0, -1}, {0, 1, 0}, {1, 2, 1}},
77{{-1, 0, 1, 0, -1, -1, 0, 0, 1}, {1, -1, 2}, {0, 0, 1}},
78{{-1, 0, 1, 0, -1, 0, 0, 0, 1}, {1, 0, 2}, {0, 0, 1}},
79{{-1, 0, 1, 0, -1, 1, 0, 0, 1}, {1, 1, 2}, {0, 0, 1}},
80{{-1, 1, -1, 0, 0, -1, 0, -1, 0}, {-1, -1, 1}, {0, -1, 1}},
81{{-1, 1, 0, 0, 1, 0, 0, -1, -1}, {-1, -2, 1}, {0, 1, 0}},
82{{-1, 1, 0, 0, 1, 0, 0, 0, -1}, {1, 2, 0}, {0, 1, 0}},
83{{-1, 1, 0, 0, 1, 0, 0, 1, -1}, {1, 2, 1}, {0, 1, 0}},
84{{-1, 1, 1, 0, 0, 1, 0, 1, 0}, {1, 1, 1}, {0, 1, 1}},
85{{0, -1, -1, -1, 0, 1, 0, 0, -1}, {-1, 1, 0}, {-1, 1, 1}},
86{{0, -1, -1, 0, -1, 0, -1, 1, 0}, {-1, 0, 1}, {-1, 1, 1}},
87{{0, -1, 0, -1, 0, 0, -1, 1, -1}, {-1, 1, 1}, {-1, 1, 0}},
88{{0, -1, 0, -1, 0, 0, 0, 0, -1}, {-1, 1, 0}, {-1, 1, 0}},
89{{0, -1, 0, -1, 0, 0, 1, -1, -1}, {1, -1, 1}, {-1, 1, 0}},
90{{0, -1, 1, -1, 0, -1, 0, 0, -1}, {-1, 1, 0}, {1, -1, 1}},
91{{0, -1, 1, 0, -1, 0, 1, -1, 0}, {1, 0, 1}, {1, -1, 1}},
92{{0, 0, -1, -1, -1, 1, -1, 0, 0}, {-1, 1, 1}, {-1, 0, 1}},
93{{0, 0, -1, 0, -1, 0, -1, 0, 0}, {-1, 0, 1}, {-1, 0, 1}},
94{{0, 0, -1, 1, -1, -1, -1, 0, 0}, {-1, -1, 1}, {-1, 0, 1}},
95{{0, 0, 1, -1, -1, -1, 1, 0, 0}, {1, -1, 1}, {1, 0, 1}},
96{{0, 0, 1, 0, -1, 0, 1, 0, 0}, {1, 0, 1}, {1, 0, 1}},
97{{0, 0, 1, 1, -1, 1, 1, 0, 0}, {1, 1, 1}, {1, 0, 1}},
98{{0, 1, -1, 0, -1, 0, -1, -1, 0}, {-1, 0, 1}, {-1, -1, 1}},
99{{0, 1, -1, 1, 0, -1, 0, 0, -1}, {1, 1, 0}, {-1, -1, 1}},
100{{0, 1, 0, 1, 0, 0, -1, -1, -1}, {-1, -1, 1}, {1, 1, 0}},
101{{0, 1, 0, 1, 0, 0, 0, 0, -1}, {1, 1, 0}, {1, 1, 0}},
102{{0, 1, 0, 1, 0, 0, 1, 1, -1}, {1, 1, 1}, {1, 1, 0}},
103{{0, 1, 1, 0, -1, 0, 1, 1, 0}, {1, 0, 1}, {1, 1, 1}},
104{{0, 1, 1, 1, 0, 1, 0, 0, -1}, {1, 1, 0}, {1, 1, 1}},
105{{1, -1, -1, 0, -1, 0, 0, 0, -1}, {1, 0, 0}, {-2, 1, 1}},
106{{1, -1, 0, 0, -1, 0, 0, 0, -1}, {1, 0, 0}, {-2, 1, 0}},
107{{1, -1, 1, 0, -1, 0, 0, 0, -1}, {1, 0, 0}, {2, -1, 1}},
108{{1, 0, -1, 0, -1, 0, 0, 0, -1}, {1, 0, 0}, {-2, 0, 1}},
109{{1, 0, 0, -1, -1, 0, -1, 0, -1}, {-2, 1, 1}, {1, 0, 0}},
110{{1, 0, 0, -1, -1, 0, 0, 0, -1}, {-2, 1, 0}, {1, 0, 0}},
111{{1, 0, 0, -1, -1, 0, 1, 0, -1}, {2, -1, 1}, {1, 0, 0}},
112{{1, 0, 0, 0, -1, 0, -1, 0, -1}, {-2, 0, 1}, {1, 0, 0}},
113{{1, 0, 0, 0, -1, 0, 0, 0, -1}, {1, 0, 0}, {1, 0, 0}},
114{{1, 0, 0, 0, -1, 0, 1, 0, -1}, {2, 0, 1}, {1, 0, 0}},
115{{1, 0, 0, 1, -1, 0, -1, 0, -1}, {-2, -1, 1}, {1, 0, 0}},
116{{1, 0, 0, 1, -1, 0, 0, 0, -1}, {2, 1, 0}, {1, 0, 0}},
117{{1, 0, 0, 1, -1, 0, 1, 0, -1}, {2, 1, 1}, {1, 0, 0}},
118{{1, 0, 1, 0, -1, 0, 0, 0, -1}, {1, 0, 0}, {2, 0, 1}},
119{{1, 1, -1, 0, -1, 0, 0, 0, -1}, {1, 0, 0}, {-2, -1, 1}},
120{{1, 1, 0, 0, -1, 0, 0, 0, -1}, {1, 0, 0}, {2, 1, 0}},
121{{1, 1, 1, 0, -1, 0, 0, 0, -1}, {1, 0, 0}, {2, 1, 1}},
122};
123
124using TwoFold = TwoFold_<void>;
125
126} // namespace impl
127
128using OpObliquity = std::pair<Op, double>;
129
130// Obliquity calculated here is the same as Le Page delta in cctbx.
131// (tested against output of lebedev_2005_perturbation.py from cctbx)
132inline double calculate_cos_obliquity(const UnitCell& reduced_cell,
133 const Vec3& d_axis, const Vec3& r_axis) {
134 // From the Le Page paper:
135 // tan(delta) = |t x tau| / |t . tau| (i.e. delta is angle 0...90deg)
136 // t = L U where U is direct-space axis
137 // tau = (L^T)^-1 h where h is reciprocal-space axis
138 // L is orthogonalization matrix, although in a different order
139 Vec3 t = reduced_cell.orth.mat.multiply(d_axis);
140 Vec3 tau = reduced_cell.frac.mat.left_multiply(r_axis);
141 // it's faster to calculate cos(delta) than tan(delta)
142 return std::min(1.0, std::fabs(t.cos_angle(tau)));
143}
144
145// Reduced cell can be from GruberVector::get_cell() after Niggli reduction.
146// max_obliq is max obliquity (delta) in degrees as defined in Le Page (1982).
147inline std::vector<OpObliquity> find_lattice_2fold_ops(const UnitCell& reduced_cell,
148 double max_obliq) {
149 std::vector<OpObliquity> ret;
150 const double cos_max_obliq = std::cos(rad(max_obliq));
151 for (const impl::TwoFoldData& row : impl::TwoFold::table) {
152 Vec3 d_axis(row.ds_axis[0], row.ds_axis[1], row.ds_axis[2]);
153 Vec3 r_axis(row.rs_axis[0], row.rs_axis[1], row.rs_axis[2]);
154 double cos_delta = calculate_cos_obliquity(reduced_cell, d_axis, r_axis);
155 if (cos_delta > cos_max_obliq) {
156 constexpr int D = Op::DEN;
157 Op op{Op::Rot{{{{row.matrix[0] * D, row.matrix[1] * D, row.matrix[2] * D}},
158 {{row.matrix[3] * D, row.matrix[4] * D, row.matrix[5] * D}},
159 {{row.matrix[6] * D, row.matrix[7] * D, row.matrix[8] * D}}}},
160 Op::Tran{{0, 0, 0}}};
161 ret.emplace_back(op, deg(std::acos(cos_delta)));
162 }
163 }
164 std::sort(ret.begin(), ret.end(),
165 [](const OpObliquity& a, const OpObliquity& b) { return a.second < b.second; });
166 return ret;
167}
168
169// Reduced cell can be from GruberVector::get_cell() after Niggli reduction.
170// max_obliq is max obliquity (delta) in degrees as defined in Le Page (1982).
171// Returns lattice symmetry except inversion.
172inline GroupOps find_lattice_symmetry_r(const UnitCell& reduced_cell, double max_obliq) {
173 std::vector<OpObliquity> gen = find_lattice_2fold_ops(reduced_cell, max_obliq);
174 std::vector<Op> genops;
175 genops.reserve(gen.size());
176 for (const OpObliquity& op_obl : gen)
177 genops.push_back(op_obl.first);
178 GroupOps go;
179 go.sym_ops.push_back(Op::identity());
180 go.cen_ops.push_back({0,0,0});
181 // cf. GroupOps::add_missing_elements()
182 if (!gen.empty())
183 go.sym_ops.push_back(gen[0].first);
184 // no need to try operator^2, we know it must be identity
185 go.add_missing_elements_part2(genops, 24, true);
186 return go;
187}
188
189// Returns lattice symmetry, but without inversion.
190inline GroupOps find_lattice_symmetry(const UnitCell& cell, char centring,
191 double max_obliq) {
192 GruberVector gv(cell, centring, true);
193 gv.niggli_reduce();
194 UnitCell reduced = gv.get_cell();
195 GroupOps gops = find_lattice_symmetry_r(reduced, max_obliq);
196 gops.change_basis_forward(*gv.change_of_basis);
197 return gops;
198}
199
200// Determine potential twinning operators.
201// Returns all operators or only unique ones (coset representatives).
202inline std::vector<Op> find_twin_laws(const UnitCell& cell,
203 const SpaceGroup* sg,
204 double max_obliq,
205 bool all_ops) {
206 if (sg == nullptr)
207 sg = &get_spacegroup_p1();
208 GroupOps go = sg->operations();
209 GroupOps lat_go = find_lattice_symmetry(cell, sg->centring_type(), max_obliq);
210 if (!go.has_same_centring(lat_go))
211 fail("find_twin_laws(): internal error");
212 std::vector<Op> ops;
213 size_t sg_symop_count = go.sym_ops.size();
214 for (const Op& op : lat_go.sym_ops)
215 if (!go.find_by_rotation(op.rot) &&
216 !go.find_by_rotation(op.negated_rot())) {
217 ops.push_back(op);
218 if (!all_ops)
219 for (size_t i = 1; i < sg_symop_count; ++i)
220 go.sym_ops.push_back(op * go.sym_ops[i]);
221 }
222 return ops;
223}
224
225} // namespace gemmi
226#endif
gemmi::get_spacegroup_p1
const SpaceGroup & get_spacegroup_p1()
Definition symmetry.hpp:1994
gemmi::OpObliquity
std::pair< Op, double > OpObliquity
Definition twin.hpp:128
gemmi::find_twin_laws
std::vector< Op > find_twin_laws(const UnitCell &cell, const SpaceGroup *sg, double max_obliq, bool all_ops)
Definition twin.hpp:202
gemmi::calculate_cos_obliquity
double calculate_cos_obliquity(const UnitCell &reduced_cell, const Vec3 &d_axis, const Vec3 &r_axis)
Definition twin.hpp:132
gemmi::deg
constexpr double deg(double angle)
Definition math.hpp:31
gemmi::find_lattice_symmetry
GroupOps find_lattice_symmetry(const UnitCell &cell, char centring, double max_obliq)
Definition twin.hpp:190
gemmi::find_lattice_2fold_ops
std::vector< OpObliquity > find_lattice_2fold_ops(const UnitCell &reduced_cell, double max_obliq)
Definition twin.hpp:147
gemmi::rad
constexpr double rad(double angle)
Definition math.hpp:32
gemmi::Vec3
Vec3_< double > Vec3
Definition math.hpp:101
gemmi::find_lattice_symmetry_r
GroupOps find_lattice_symmetry_r(const UnitCell &reduced_cell, double max_obliq)
Definition twin.hpp:172
gemmi::fail
void fail(const std::string &msg)
Definition fail.hpp:59
gemmi::GroupOps::change_basis_forward
void change_basis_forward(const Op &cob)
Definition symmetry.hpp:573
gemmi::Op::Rot
std::array< std::array< int, 3 >, 3 > Rot
Definition symmetry.hpp:55
gemmi::Op::DEN
static constexpr int DEN
Definition symmetry.hpp:54
gemmi::Op::Tran
std::array< int, 3 > Tran
Definition symmetry.hpp:56
gemmi::Op::identity
static constexpr Op identity()
Definition symmetry.hpp:177
gemmi::UnitCell
Unit cell.
Definition unitcell.hpp:139
Generated by doxygen 1.9.8