takaba_k
c1b57a83e0
Build Packages / build:rpm (ubuntu2404_nocuda) (push) Successful in 12m37s
Build Packages / build:rpm (rocky8_nocuda) (push) Successful in 13m37s
Build Packages / build:rpm (rocky9_nocuda) (push) Successful in 14m14s
Build Packages / build:rpm (rocky8_sls9) (push) Successful in 14m18s
Build Packages / build:rpm (rocky8) (push) Successful in 14m42s
Build Packages / Build documentation (push) Successful in 40s
Build Packages / Generate python client (push) Successful in 1m8s
Build Packages / Create release (push) Has been skipped
Build Packages / build:rpm (rocky9_sls9) (push) Successful in 16m6s
Build Packages / build:rpm (ubuntu2204_nocuda) (push) Successful in 17m3s
Build Packages / build:rpm (rocky9) (push) Successful in 11m4s
Build Packages / build:rpm (ubuntu2404) (push) Successful in 9m48s
Build Packages / build:rpm (ubuntu2204) (push) Successful in 11m21s
Build Packages / Unit tests (push) Successful in 58m23s
- Search only Bravais types tentatively. The second step will be added - Accurate sigma estimation on merging with partiality-weighting - Additional log output: to be more organized
375 lines
12 KiB
C++
375 lines
12 KiB
C++
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "LatticeSearch.h"
|
|
#include <gemmi/cellred.hpp>
|
|
|
|
#include <cmath>
|
|
|
|
struct NiggliClass {
|
|
int number;
|
|
int type;
|
|
bool cond_AB;
|
|
bool cond_BC;
|
|
double cond_D;
|
|
double cond_E;
|
|
double cond_F;
|
|
bool cond_DEF;
|
|
bool cond_2DF;
|
|
gemmi::Mat33 reindex;
|
|
gemmi::CrystalSystem system;
|
|
char centering;
|
|
int lowest_spacegroup_number;
|
|
};
|
|
|
|
LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance, double angle_tolerance, int ngc_id) {
|
|
UnitCell uc = L.GetUnitCell();
|
|
gemmi::UnitCell g_uc(uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
|
|
gemmi::GruberVector g_vec(g_uc, 'P', true);
|
|
g_vec.niggli_reduce();
|
|
|
|
CrystalLattice L_niggli = L;
|
|
if (g_vec.change_of_basis)
|
|
L_niggli = L.Multiply(gemmi::rot_as_mat33(g_vec.change_of_basis->rot).transpose());
|
|
|
|
double A = g_vec.A;
|
|
double B = g_vec.B;
|
|
double C = g_vec.C;
|
|
double D = g_vec.xi / 2;
|
|
double E = g_vec.eta / 2;
|
|
double F = g_vec.zeta / 2;
|
|
|
|
std::vector<NiggliClass> niggli_classes = {
|
|
{
|
|
1, 1,
|
|
true, true, A / 2, A / 2, A / 2, false, false,
|
|
gemmi::Mat33{1, -1, 1, 1, 1, -1, -1, 1, 1},
|
|
gemmi::CrystalSystem::Cubic, 'F', 196
|
|
},
|
|
{
|
|
2, 1,
|
|
true, true, D, D, D, false, false,
|
|
{1, -1, 0, -1, 0, 1, -1, -1, -1},
|
|
gemmi::CrystalSystem::Trigonal, 'R', 146
|
|
},
|
|
{
|
|
3, 2,
|
|
true, true, 0, 0, 0, false, false,
|
|
gemmi::Mat33{1, 0, 0, 0, 1, 0, 0, 0, 1},
|
|
gemmi::CrystalSystem::Cubic, 'P', 195
|
|
},
|
|
{
|
|
5, 2,
|
|
true, true, -A / 3, -A / 3, -A / 3, false, false,
|
|
gemmi::Mat33{1, 0, 1, 1, 1, 0, 0, 1, 1},
|
|
gemmi::CrystalSystem::Cubic, 'I', 197
|
|
},
|
|
{
|
|
4, 2,
|
|
true, true, D, D, D, false, false,
|
|
{1, -1, 0, -1, 0, 1, -1, -1, -1},
|
|
gemmi::CrystalSystem::Trigonal, 'R', 146
|
|
},
|
|
{
|
|
6, 2,
|
|
true, true, D, D, F, true, false,
|
|
{0, 1, 1, 1, 0, 1, 1, 1, 0},
|
|
gemmi::CrystalSystem::Tetragonal, 'I', 79
|
|
},
|
|
{
|
|
7, 2,
|
|
true, true, D, E, E, true, false,
|
|
{1, 0, 1, 1, 1, 0, 0, 1, 1},
|
|
gemmi::CrystalSystem::Tetragonal, 'I', 79
|
|
},
|
|
{
|
|
8, 2,
|
|
true, true, D, E, F, true, false,
|
|
{-1, -1, 0, -1, 0, -1, 0, -1, -1},
|
|
gemmi::CrystalSystem::Orthorhombic, 'I', 23
|
|
},
|
|
{
|
|
9, 1,
|
|
true, false, A / 2, A / 2, A / 2, false, false,
|
|
{1, 0, 0, -1, 1, 0, -1, -1, -3},
|
|
gemmi::CrystalSystem::Trigonal, 'R', 146
|
|
},
|
|
|
|
{
|
|
10, 1,
|
|
true, false, D, D, F, false, false,
|
|
{1, 1, 0, 1, -1, 0, 0, 0, -1},
|
|
gemmi::CrystalSystem::Monoclinic, 'C', 5
|
|
},
|
|
{
|
|
11, 2,
|
|
true, false, 0, 0, 0, false, false,
|
|
{1, 0, 0, 0, 1, 0, 0, 0, 1},
|
|
gemmi::CrystalSystem::Tetragonal, 'P', 75
|
|
},
|
|
{
|
|
12, 2,
|
|
true, false, 0, 0, -A / 2, false, false,
|
|
{1, 0, 0, 0, 1, 0, 0, 0, 1},
|
|
gemmi::CrystalSystem::Hexagonal, 'P', 143
|
|
},
|
|
{
|
|
13, 2,
|
|
true, false, 0, 0, F, false, false,
|
|
{1, 1, 0, -1, 1, 0, 0, 0, 1},
|
|
gemmi::CrystalSystem::Orthorhombic, 'C', 20
|
|
},
|
|
{
|
|
15, 2,
|
|
true, false, -A / 2, -A / 2, 0, false, false,
|
|
{1, 0, 0, 0, 1, 0, 1, 1, 2},
|
|
gemmi::CrystalSystem::Tetragonal, 'I', 79
|
|
},
|
|
{
|
|
16, 2,
|
|
true, false, D, D, F, true, false,
|
|
{-1, -1, 0, 1, -1, 0, 1, 1, 2},
|
|
gemmi::CrystalSystem::Orthorhombic, 'F', 22
|
|
},
|
|
{
|
|
14, 2,
|
|
true, false, D, D, F, false, false,
|
|
{1, 1, 0, -1, 1, 0, 0, 0, 1},
|
|
gemmi::CrystalSystem::Monoclinic, 'C', 5
|
|
},
|
|
{
|
|
17, 2,
|
|
true, false, D, E, F, true, false,
|
|
{1, -1, 0, 1, 1, 0, -1, 0, -1},
|
|
gemmi::CrystalSystem::Monoclinic, 'C', 5
|
|
},
|
|
{
|
|
18, 1,
|
|
false, true, A / 4, A / 2, A / 2, false, false,
|
|
{0, -1, 1, 1, -1, -1, 1, 0, 0},
|
|
gemmi::CrystalSystem::Tetragonal, 'I', 79
|
|
},
|
|
{
|
|
19, 1,
|
|
false, true, D, A / 2, A / 2, false, false,
|
|
{-1, 0, 0, 0, -1, 1, -1, 1, 1},
|
|
gemmi::CrystalSystem::Orthorhombic, 'I', 23
|
|
},
|
|
{
|
|
20, 1,
|
|
false, true, D, E, E, false, false,
|
|
{0, 1, 1, 0, 1, -1, -1, 0, 0},
|
|
gemmi::CrystalSystem::Monoclinic, 'C', 5
|
|
},
|
|
{
|
|
21, 2,
|
|
false, true, 0, 0, 0, false, false,
|
|
{0, 1, 0, 0, 0, 1, 1, 0, 0},
|
|
gemmi::CrystalSystem::Tetragonal, 'P', 75
|
|
},
|
|
{
|
|
22, 2,
|
|
false, true, -B / 2, 0, 0, false, false,
|
|
{0, 1, 0, 0, 0, 1, 1, 0, 0},
|
|
gemmi::CrystalSystem::Hexagonal, 'P', 143
|
|
},
|
|
{
|
|
23, 2,
|
|
false, true, D, 0, 0, false, false,
|
|
{0, 1, 1, 0, -1, 1, 1, 0, 0},
|
|
gemmi::CrystalSystem::Orthorhombic, 'C', 20
|
|
},
|
|
{
|
|
24, 2,
|
|
false, true, D, -A / 3, -A / 3, true, false,
|
|
{1, 2, 1, 0, -1, 1, 1, 0, 0},
|
|
gemmi::CrystalSystem::Trigonal, 'R', 146
|
|
},
|
|
{
|
|
25, 2,
|
|
false, true, D, E, E, false, false,
|
|
{0, 1, 1, 0, -1, 1, 1, 0, 0},
|
|
gemmi::CrystalSystem::Monoclinic, 'C', 5
|
|
},
|
|
{
|
|
26, 1,
|
|
false, false, A / 4, A / 2, A / 2, false, false,
|
|
{1, 0, 0, -1, 2, 0, -1, 0, 2},
|
|
gemmi::CrystalSystem::Orthorhombic, 'F', 22
|
|
},
|
|
{
|
|
27, 1,
|
|
false, false, D, A / 2, A / 2, false, false,
|
|
{-1, 2, 0, -1, 0, 0, 0, -1, 1},
|
|
gemmi::CrystalSystem::Monoclinic, 'C', 5
|
|
},
|
|
{
|
|
28, 1,
|
|
false, false, D, A / 2, 2 * D, false, false,
|
|
{-1, 0, 0, -1, 0, 2, 0, 1, 0},
|
|
gemmi::CrystalSystem::Monoclinic, 'C', 5
|
|
},
|
|
{
|
|
29, 1,
|
|
false, false, D, 2 * D, A / 2, false, false,
|
|
{1, 0, 0, 1, -2, 0, 0, 0, -1},
|
|
gemmi::CrystalSystem::Monoclinic, 'C', 5
|
|
},
|
|
{
|
|
30, 1,
|
|
false, false, B / 2, E, 2 * E, false, false,
|
|
{0, 1, 0, 0, 1, -2, -1, 0, 0},
|
|
gemmi::CrystalSystem::Monoclinic, 'C', 5
|
|
},
|
|
|
|
{
|
|
31, 1,
|
|
false, false, D, E, F, false, false,
|
|
{1, 0, 0, 0, 1, 0, 0, 0, 1},
|
|
gemmi::CrystalSystem::Triclinic, 'P', 75
|
|
},
|
|
|
|
{
|
|
32, 2,
|
|
false, false, 0, 0, 0, false, false,
|
|
{1, 0, 0, 0, 1, 0, 0, 0, 1},
|
|
gemmi::CrystalSystem::Orthorhombic, 'P', 16
|
|
},
|
|
{
|
|
40, 2,
|
|
false, false, -B / 2, 0, 0, false, false,
|
|
{0, -1, 0, -1, 0, 0, 0, 0, -1},
|
|
gemmi::CrystalSystem::Orthorhombic, 'C', 20
|
|
},
|
|
{
|
|
35, 2,
|
|
false, false, D, 0, 0, false, false,
|
|
{0, -1, 0, -1, 0, 0, 0, 0, -1},
|
|
gemmi::CrystalSystem::Monoclinic, 'P', 3
|
|
},
|
|
{
|
|
36, 2,
|
|
false, false, 0, -A / 2, 0, false, false,
|
|
{1, 0, 0, -1, 0, -2, 0, 1, 0},
|
|
gemmi::CrystalSystem::Orthorhombic, 'C', 20
|
|
},
|
|
{
|
|
33, 2,
|
|
false, false, 0, E, 0, false, false,
|
|
{1, 0, 0, 0, 1, 0, 0, 0, 1},
|
|
gemmi::CrystalSystem::Monoclinic, 'P', 3
|
|
},
|
|
{
|
|
38, 2,
|
|
false, false, 0, 0, -A / 2, false, false,
|
|
{-1, 0, 0, 1, 2, 0, 0, 0, -1},
|
|
gemmi::CrystalSystem::Orthorhombic, 'C', 20
|
|
},
|
|
{
|
|
34, 2,
|
|
false, false, 0, 0, F, false, false,
|
|
{-1, 0, 0, 0, 0, -1, 0, -1, 0},
|
|
gemmi::CrystalSystem::Monoclinic, 'P', 3
|
|
},
|
|
{
|
|
42, 2,
|
|
false, false, -B / 2, -A / 2, 0, false, false,
|
|
{-1, 0, 0, 0, -1, 0, 1, 1, 2},
|
|
gemmi::CrystalSystem::Orthorhombic, 'I', 23
|
|
},
|
|
{
|
|
41, 2,
|
|
false, false, -B / 2, E, 0, false, false,
|
|
{0, -1, -2, 0, -1, 0, -1, 0, 0},
|
|
gemmi::CrystalSystem::Monoclinic, 'C', 5
|
|
},
|
|
{
|
|
37, 2,
|
|
false, false, D, -A / 2, 0, false, false,
|
|
{1, 0, 2, 1, 0, 0, 0, 1, 0},
|
|
gemmi::CrystalSystem::Monoclinic, 'C', 5
|
|
},
|
|
{
|
|
39, 2,
|
|
false, false, D, 0, -A / 2, false, false,
|
|
{-1, -2, 0, -1, 0, 0, 0, 0, -1},
|
|
gemmi::CrystalSystem::Monoclinic, 'C', 5
|
|
},
|
|
{
|
|
44, 2,
|
|
false, false, D, E, F, false, false,
|
|
{1, 0, 0, 0, 1, 0, 0, 0, 1},
|
|
gemmi::CrystalSystem::Triclinic, 'P', 1
|
|
}
|
|
};
|
|
|
|
const auto uc_reduced = L_niggli.GetUnitCell();
|
|
|
|
std::vector<NiggliClass> niggli_class_selected;
|
|
if (ngc_id != -1)
|
|
niggli_class_selected = {niggli_classes[ngc_id]};
|
|
else
|
|
niggli_class_selected = niggli_classes;
|
|
|
|
for (const auto &c: niggli_class_selected) {
|
|
if (c.type == 1 && uc_reduced.beta >= 90 - angle_tolerance )
|
|
continue;
|
|
|
|
bool ok = true;
|
|
|
|
if (c.cond_AB && fabs((uc_reduced.a - uc_reduced.b) / (0.5 * (uc_reduced.a + uc_reduced.b))) > dist_tolerance)
|
|
ok = false;
|
|
if (c.cond_BC && fabs((uc_reduced.b - uc_reduced.c) / (0.5 * (uc_reduced.b + uc_reduced.c))) > dist_tolerance)
|
|
ok = false;
|
|
|
|
double expected_alpha = acos(c.cond_D / sqrt(B*C)) * 180 / M_PI;
|
|
double expected_beta = acos(c.cond_E / sqrt(A*C)) * 180 / M_PI;
|
|
double expected_gamma = acos(c.cond_F / sqrt(A*B)) * 180 / M_PI;
|
|
|
|
if (fabs(expected_alpha - uc_reduced.alpha) > angle_tolerance)
|
|
ok = false;
|
|
if (fabs(expected_beta - uc_reduced.beta) > angle_tolerance)
|
|
ok = false;
|
|
if (fabs(expected_gamma - uc_reduced.gamma) > angle_tolerance)
|
|
ok = false;
|
|
|
|
double tmp1 = 2.0 * fabs(D + E + F);
|
|
double tmp2 = A + B;
|
|
|
|
if (c.cond_DEF && fabs((tmp1 - tmp2) / (0.5 * (tmp1 + tmp2))) > dist_tolerance)
|
|
ok = false;
|
|
|
|
if (ok) {
|
|
return LatticeSearchResult{
|
|
.niggli_class = c.number,
|
|
.primitive_reduced = L,
|
|
.conventional = L_niggli.Multiply(c.reindex),
|
|
.system = c.system,
|
|
.centering = c.centering,
|
|
.lowest_spacegroup_number = c.lowest_spacegroup_number,
|
|
};
|
|
}
|
|
}
|
|
|
|
return LatticeSearchResult{
|
|
.niggli_class = 44,
|
|
.primitive_reduced = L,
|
|
.conventional = L,
|
|
.system = gemmi::CrystalSystem::Triclinic,
|
|
.centering = 'P',
|
|
.lowest_spacegroup_number = 1,
|
|
};
|
|
}
|
|
|
|
std::vector<LatticeSearchResult> LatticeCandidates(const CrystalLattice &L, double dist_tolerance, double angle_tolerance) {
|
|
std::vector<LatticeSearchResult> result;
|
|
int n_possible_niggli_classes = 44;
|
|
|
|
for (int id = 1; id <= n_possible_niggli_classes; id++) {
|
|
LatticeSearchResult lattice_found = LatticeSearch(L, dist_tolerance, angle_tolerance, id);
|
|
if (lattice_found.niggli_class != 44)
|
|
result.push_back(lattice_found);
|
|
}
|
|
|
|
return result;
|
|
} |