jfj-process: better spacegroup estimation - intermediate
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
This commit is contained in:
2026-04-16 09:38:49 +02:00
parent de6646b95e
commit c1b57a83e0
6 changed files with 126 additions and 72 deletions
+67 -45
View File
@@ -19,9 +19,10 @@ struct NiggliClass {
gemmi::Mat33 reindex;
gemmi::CrystalSystem system;
char centering;
int lowest_spacegroup_number;
};
LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance, double angle_tolerance) {
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);
@@ -43,268 +44,274 @@ LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
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'
gemmi::CrystalSystem::Triclinic, 'P', 1
}
};
const auto uc_reduced = L_niggli.GetUnitCell();
for (const auto &c: niggli_classes) {
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;
@@ -339,6 +346,7 @@ LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance
.conventional = L_niggli.Multiply(c.reindex),
.system = c.system,
.centering = c.centering,
.lowest_spacegroup_number = c.lowest_spacegroup_number,
};
}
}
@@ -349,5 +357,19 @@ LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance
.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;
}
@@ -16,8 +16,11 @@ struct LatticeSearchResult {
CrystalLattice conventional;
gemmi::CrystalSystem system = gemmi::CrystalSystem::Triclinic;
char centering = 'P'; // 'P','A','B','C','I','F','R'
int lowest_spacegroup_number = 1;
};
LatticeSearchResult LatticeSearch(const CrystalLattice& L, double dist_tolerance = 0.03, double angle_tolerance_deg = 3);
LatticeSearchResult LatticeSearch(const CrystalLattice& L, double dist_tolerance = 0.03, double angle_tolerance_deg = 3, int ngc_id = -1);
std::vector<LatticeSearchResult> LatticeCandidates(const CrystalLattice& L, double dist_tolerance = 0.03, double angle_tolerance_deg = 3);
#endif //JFJOCH_LATTICESEARCH_H
+10 -15
View File
@@ -244,6 +244,7 @@ namespace {
int hkl_slot;
double I_corr;
double sigma_corr;
double weight;
};
void select_reflections_by_quasi_random(const std::vector<ObsRef> &obs,
@@ -516,18 +517,19 @@ namespace {
struct HKLAccum {
double sum_wI = 0.0;
double sum_w = 0.0;
double sum_wI2 = 0.0;
int count = 0;
// double sum_wI2 = 0.0;
double sum_wsigma2 = 0.0;
// int count = 0;
};
std::vector<HKLAccum> accum(nhkl);
for (const auto &co: corr_obs) {
const double w = 1.0 / (co.sigma_corr * co.sigma_corr);
// const double w = 1.0 / (co.sigma_corr * co.sigma_corr);
const double w = co.weight / (co.sigma_corr * co.sigma_corr);
auto &a = accum[co.hkl_slot];
a.sum_wI += w * co.I_corr;
a.sum_w += w;
a.sum_wI2 += w * co.I_corr * co.I_corr;
a.count++;
a.sum_wsigma2 += w * w * co.sigma_corr * co.sigma_corr;
}
for (int h = 0; h < nhkl; ++h) {
@@ -536,15 +538,7 @@ namespace {
continue;
out.merged[h].I = a.sum_wI / a.sum_w;
const double sigma_stat = std::sqrt(1.0 / a.sum_w);
if (a.count >= 2) {
const double var_internal = a.sum_wI2 - (a.sum_wI * a.sum_wI) / a.sum_w;
const double sigma_int = std::sqrt(var_internal / (a.sum_w * (a.count - 1)));
out.merged[h].sigma = std::max(sigma_stat, sigma_int);
} else {
out.merged[h].sigma = sigma_stat;
}
out.merged[h].sigma = std::sqrt(a.sum_wsigma2) / a.sum_w;
}
}
@@ -751,7 +745,8 @@ namespace {
corr_obs.push_back({
o.hkl_slot,
static_cast<double>(r.I) / correction,
o.sigma / correction
o.sigma / correction,
1 / correction,
});
}
}
@@ -8,6 +8,7 @@
#include <cmath>
#include <cctype>
#include <iomanip>
#include <iostream>
#include <limits>
#include <sstream>
#include <tuple>
@@ -394,8 +395,16 @@ SearchSpaceGroupResult SearchSpaceGroup(
if (merged.empty())
return result;
std::vector<gemmi::SpaceGroup> candidates;
const auto lookup = BuildMergedLookup(merged, opt);
const auto candidates = EnumerateCandidateSpaceGroups(opt.crystal_system, opt.centering);
if (opt.candidate_space_groups.has_value())
candidates = opt.candidate_space_groups.value();
else
candidates = EnumerateCandidateSpaceGroups(opt.crystal_system, opt.centering);
std::cout << "No. of SG-search candidates: " << candidates.size() << std::endl;
for (const auto& sg : candidates) {
SpaceGroupCandidateScore score{.space_group = sg};
@@ -91,6 +91,8 @@ struct SearchSpaceGroupOptions {
// Acceptance thresholds for absent/allowed <I/sig> ratios.
double max_absent_to_allowed_i_over_sigma_ratio = 0.20;
double max_absent_to_allowed_i_over_sigma_ratio_in_any_bin = 0.50;
std::optional<std::vector<gemmi::SpaceGroup>> candidate_space_groups;
};
struct SearchSpaceGroupResult {
+33 -10
View File
@@ -28,6 +28,7 @@
#include "../compression/JFJochCompressor.h"
#include "../image_analysis/scale_merge/FrenchWilson.h"
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
#include "../image_analysis/lattice_search/LatticeSearch.h"
#include "../image_analysis/WriteMmcif.h"
void print_usage(Logger &logger) {
@@ -94,6 +95,20 @@ void print_statistics(Logger &logger, const MergeStatistics &stats) {
logger.Info("");
}
void print_unitcell(Logger &logger, const CrystalLattice &lattice, bool print_vector) {
auto vec0 = lattice.Vec0();
auto vec1 = lattice.Vec1();
auto vec2 = lattice.Vec2();
auto uc = lattice.GetUnitCell();
if (print_vector) {
logger.Info("Lattice vec0: {:.3f}, {:.3f}, {:.3f}", vec0.x, vec0.y, vec0.z);
logger.Info("Lattice vec1: {:.3f}, {:.3f}, {:.3f}", vec1.x, vec1.y, vec1.z);
logger.Info("Lattice vec2: {:.3f}, {:.3f}, {:.3f}", vec2.x, vec2.y, vec2.z);
}
logger.Info("Lattice: a={:.2f} b={:.2f} c={:.2f} alpha={:.2f} beta={:.2f} gamma={:.2f}",
uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
}
std::optional<UnitCell> parse_unit_cell_arg(const char* arg) {
if (!arg)
return std::nullopt;
@@ -638,7 +653,9 @@ int main(int argc, char **argv) {
.niggli_class = rotation_indexer_ret->search_result.niggli_class,
.crystal_system = rotation_indexer_ret->search_result.system
};
logger.Info("Rotation Indexing found lattice");
logger.Info("Rotation Indexing found lattice: {} {}",
gemmi::crystal_system_str(rotation_indexer_ret->search_result.system),
rotation_indexer_ret->search_result.centering);
}
// --- Optional: run scaling (mosaicity refinement) on accumulated reflections ---
@@ -652,6 +669,7 @@ int main(int argc, char **argv) {
scale_opts.merge_friedel = !anomalous_mode;
scale_opts.d_min_limit_A = d_min_scale_merge.value_or(0.0);
scale_opts.refine_wedge = refine_wedge;
// scale_opts.min_partiality_for_merge = 0.4;
if (filter_min_per_bin.has_value()) {
scale_opts.selection = true;
scale_opts.filter_min_per_bin = filter_min_per_bin.value();
@@ -681,6 +699,17 @@ int main(int argc, char **argv) {
if (scale_result && !fixed_space_group) {
logger.Info("Searching for space group from P1-merged reflections ...");
std::vector<gemmi::SpaceGroup> candidates_preset;
if (rotation_indexer_ret.has_value()) {
std::vector<LatticeSearchResult> lattice_candidates = LatticeCandidates(rotation_indexer_ret->lattice);
for (const auto &lc: lattice_candidates) {
logger.Info("Possible lattice : {} {} ({})", gemmi::crystal_system_str(lc.system), lc.centering, lc.lowest_spacegroup_number);
print_unitcell(logger, lc.conventional, false);
candidates_preset.push_back(*gemmi::find_spacegroup_by_number(lc.lowest_spacegroup_number));
}
candidates_preset.push_back(*gemmi::find_spacegroup_by_number(1));
}
SearchSpaceGroupOptions sg_opts;
sg_opts.crystal_system.reset();
sg_opts.centering = '\0';
@@ -691,6 +720,8 @@ int main(int argc, char **argv) {
sg_opts.min_pairs_per_operator = 20;
sg_opts.min_total_compared = 100;
sg_opts.test_systematic_absences = true;
if (~candidates_preset.empty())
sg_opts.candidate_space_groups = candidates_preset;
auto sg_search_start = std::chrono::steady_clock::now();
const auto sg_search = SearchSpaceGroup(scale_result->merged, sg_opts);
@@ -860,14 +891,6 @@ int main(int argc, char **argv) {
latt = latt.Multiply(rot);
}
auto vec0 = rotation_indexer_ret->lattice.Vec0();
auto vec1 = rotation_indexer_ret->lattice.Vec1();
auto vec2 = rotation_indexer_ret->lattice.Vec2();
auto uc = rotation_indexer_ret->lattice.GetUnitCell();
logger.Info("Lattice vec0: {:.3f}, {:.3f}, {:.3f}", vec0.x, vec0.y, vec0.z);
logger.Info("Lattice vec1: {:.3f}, {:.3f}, {:.3f}", vec1.x, vec1.y, vec1.z);
logger.Info("Lattice vec2: {:.3f}, {:.3f}, {:.3f}", vec2.x, vec2.y, vec2.z);
logger.Info("Lattice: a={:.2f} b={:.2f} c={:.2f} alpha={:.2f} beta={:.2f} gamma={:.2f}",
uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
print_unitcell(logger, rotation_indexer_ret->lattice, true);
}
}