From 493be2cf1946b01e859151e3f210f8c0c05f7553 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 30 Aug 2026 03:27:56 +0200 Subject: [PATCH] symmetry: enumerate the settings the cell can host, not only the reference ones The space-group search offered a candidate only if gemmi calls it the reference setting. A setting is a statement about direction, so that restricted the search to the axes the convention chose: a crystal whose 2-fold lies on c had no rung between P1 and 222 and fell to P1, and one whose screws lie on b and c was reported as the group with a single screw on c - the wrong group, not a lower one, because the candidate that predicts a subset of the real absences and nothing else wins on no evidence at all. Both stages now enumerate more, under refusals rather than thresholds. Stage B offers the non-reference settings of the chosen point group. Their rotation set is equal to the chosen one, not merely contained in it, so this cannot raise the symmetry; what it adds is a screw or a centring on the axis the data show it on. A candidate is offered only if the cell's own metric admits the rotations its setting names, and one predicting exactly the absences another candidate already predicts is dropped as the same hypothesis under a second name. A non-reference candidate whose centring class this merge does not contain is refused outright: the reference path may adopt an untested centring because the caller's centred-lattice re-test backs it, and a non-reference setting has no such backing. Stage A offers the rotation sets no reference setting carries - the a-unique and c-unique monoclinic 2-folds, and the two rhombohedral-axes trigonal groups - and only those, so every point group reachable before is still reached by the same group in the same setting. A rung reached only that way may be ADOPTED but does not judge anything else: it is skipped when the reference chi^2 is formed and when a higher promotion's parents are collected. Without that it made higher promotions strictly harder - a promotion answers to the most damning of its parents, and offering two more order-2 subgroups of 222 refused 222 and 432 on crystals that had them. The run report gains SPACE_GROUP_NAME beside SPACE_GROUP_NUMBER, since the number alone does not say which axes a group's symmetry lies on, and a group that is not a reference setting is printed as its extended Hermann-Mauguin name. The two-arm reconciliation compares point groups by that name rather than by number for the same reason. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01T3yNBXk4wKdMZy1ak2NY7f --- docs/CHANGELOG.md | 7 + docs/CPU_DATA_ANALYSIS.md | 2 + docs/RUGNUX.md | 1 + .../scale_merge/SearchSpaceGroup.cpp | 185 ++++++++++++++---- image_analysis/scale_merge/SearchSpaceGroup.h | 24 +++ rugnux/ResultReport.cpp | 6 +- rugnux/Rugnux.cpp | 40 ++-- tests/SearchSpaceGroupTest.cpp | 123 ++++++++++++ 8 files changed, 340 insertions(+), 48 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 4fd9f4b82..eba08924c 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,13 @@ ### 1.0.0-rc.166 * `rugnux` reads a PILATUS miniCBF rotation sweep natively, with no conversion: naming any frame - or the directory holding it - processes the whole sweep that frame's template belongs to. `--mode scale` still needs a `_process.h5`. * A master whose companion `_meta.h5` was not kept still opens: the fields an Eiger master links into that file are treated as absent rather than as a read error, and a run that finds no saturation value anywhere says so and carries on instead of refusing the dataset. +* A miniCBF sweep takes the mounting from the imgCIF axis table its header carries: which laboratory direction the image's columns and rows run along, and which the spindle turns about. Where there is no table, a `+SLOW` on the `Oscillation_axis` line still says the spindle runs along the image's slow direction. All three were previously assumed, and two instruments are not what was assumed. +* Naming a miniCBF frame with no directory - a frame in the working directory - finds its sweep instead of reporting that no images were found. +* A detector swung out on a 2theta arm is placed where it stands, from the depends_on transformation chain of an NXmx master or the `Detector_2theta` line of a miniCBF header; both were previously read and then ignored. +* `rugnux` writes the unmerged MTZ `_unmerged.mtz` by default; `--no-export-unmerged` skips it. +* Every rotation run that determines its own space group also writes `_P1.mtz`, the same observations merged in P1, so a wrong space group can be re-merged, re-solved or re-refined without processing the images again; `--no-p1-crosscheck` declines it, and a run given `-S` writes nothing because its space group's centring absences were never integrated. +* The space-group search names the setting the data support, so a screw or a 2-fold on the a or c axis is reported as such instead of costing the crystal its space group. +* The run report carries `SPACE_GROUP_NAME` beside `SPACE_GROUP_NUMBER`, since the number alone does not say which axes a group's symmetry lies on. * `jfjoch_viewer` opens PILATUS miniCBF sweeps - naming any frame opens the whole sweep - and can run a processing job on one. * A miniCBF sweep takes its rotation axis from the goniometer angles the header states, instead of assuming the axis every such file was previously assumed to have. * A miniCBF sweep takes the mounting from the imgCIF axis table its header carries - which laboratory direction the image's columns and rows run along, and which the spindle turns about - and, where there is no table, from a `+SLOW` on the `Oscillation_axis` line, which says the spindle runs along the image's slow direction. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index c420a2ec7..74df89280 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -1081,6 +1081,8 @@ $$\frac{I}{\sigma} \;\ge\; \frac{T}{\sqrt{1 + (T/\mathrm{ISa})^2}}$$ The converted cut lies strictly below $\mathrm{ISa}$ for every $\mathrm{ISa}$, so it is always reachable by the reflection the ceiling binds hardest, and it is within 1 % of $T$ on any merge with $\mathrm{ISa} \ge 21$ — a healthy merge is left exactly where it was. Multiplicity is taken as 1 deliberately rather than estimated, for the same reason. Where the merge reports no ISa (or $b = 0$) the cut is used as it stands. +The candidates are enumerated in **every setting the refined cell can host**, not only in the settings the International Tables call the reference one. A setting is a statement about direction — $P112_1$ puts its $2_1$ on $c$ where $P12_11$ puts it on $b$, and both are space group 4 — so a search restricted to reference settings can only ever put a screw or a centring on the axis the convention chose, whatever the data say. Two things bound the widened set. A candidate is offered only if the cell's own metric admits the rotations its setting names ($\alpha=\beta=90^\circ$ for a $2$-fold on $c$), and a candidate predicting exactly the absences another candidate already predicts is dropped as the same hypothesis under a second name. A non-reference setting is additionally refused when its centring class holds no reflection in this merge, since there is then nothing to confirm it with. Because the candidates of a point group all share its rotations, this cannot raise the symmetry: it decides *which axes* carry the screws and the centring, never how many operators there are. Where a determined group is not the reference setting, the report and the reflection files name it by its full Hermann–Mauguin symbol, and `SPACE_GROUP_NAME` in the run report carries it beside the number, which alone would be ambiguous. + Several space groups may share an absence pattern exactly. Where they do, the search scores them identically and **all of them are named** in the result rather than one being reported as the answer: some are enantiomorph pairs, which merged intensities cannot distinguish in principle, and others differ only by a screw condition that the centering condition already implies, so the screw has no observable signature at all. The representative reported first is the lowest space-group number, which is a convention and not a measurement. The Lorentz factor $\zeta$ (§8.3) governs how well a reflection can be measured, so when the spindle lies in a plane of the lattice, an operator permuting the two in-plane axes samples a different mixture of measurement qualities than one that only flips signs. The search is therefore run a second time on a merge of only the well-measured observations (`--search-min-zeta`, rotation default 0.85), both answers are reported, and **where they disagree the merge of all the observations decides**. The filter discards 40–80 % of the observations, which can starve an operator correlation the full merge confirms and can equally leave an operator confirmed that the full merge refuses, so the decision — the point group as well as the absences, which live in the weak reflections the filter removes — rests on the arm with every observation behind it. A tie (same order, different symmetry) is reported with both candidates named, for trying in molecular replacement. diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 25d03fa79..155a2ff15 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -223,6 +223,7 @@ it in a script, every value is a `KEY= value` line: ``` grep '^SPACE_GROUP_NUMBER= ' myrun_report.txt +grep '^SPACE_GROUP_NAME= ' myrun_report.txt # the setting, which the number does not name grep '^UNIT_CELL_CONSTANTS= ' myrun_report.txt grep '^INCLUDE_RESOLUTION_RANGE= ' myrun_report.txt grep '^ISA= ' myrun_report.txt diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 39995c9a9..cc6c37e28 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -153,6 +153,57 @@ namespace { return out; } + // Whether a CELL can host a group's rotations. A setting is a statement about DIRECTION - + // P 1 1 2 puts the 2-fold on c and needs alpha = beta = 90, where P 1 2 1 puts it on b and needs + // alpha = gamma = 90 - so a candidate offered in a setting the metric does not have would be + // merged on axes the crystal does not have. Compared on the metric tensor, each element against + // its own scale, which makes the bound a tolerance on a lattice ANGLE and on an axis-length ratio; + // gemmi's own is_compatible_with_groupops takes one absolute eps in A^2 instead, which means + // something different on a 30 A cell and on a 300 A one. + constexpr double CELL_SETTING_TOLERANCE = 2e-3; // ~0.11 deg on an angle, ~0.1% on an axis ratio + + bool CellHostsRotations(const gemmi::UnitCell& cell, const gemmi::GroupOps& gops) { + const auto g = cell.metric_tensor(); + const double G[3][3] = {{g.u11, g.u12, g.u13}, {g.u12, g.u22, g.u23}, {g.u13, g.u23, g.u33}}; + for (const gemmi::Op& op : gops.sym_ops) { + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) { + double v = 0; + for (int p = 0; p < 3; ++p) + for (int q = 0; q < 3; ++q) + v += static_cast(op.rot[p][i]) * G[p][q] * op.rot[q][j]; + v /= static_cast(gemmi::Op::DEN) * gemmi::Op::DEN; + if (std::fabs(v - G[i][j]) > CELL_SETTING_TOLERANCE * std::sqrt(G[i][i] * G[j][j])) + return false; + } + } + return true; + } + + // How a candidate is NAMED in the report. short_name() is what the reference settings have + // always been printed as, and it stays that; it cannot express a setting, though - P 1 2 1 and + // P 1 1 2 are both "P2" - so a non-reference setting is printed as its extended Hermann-Mauguin + // name, which is the only faithful one. + std::string SettingName(const gemmi::SpaceGroup& sg) { + return sg.is_reference_setting() ? sg.short_name() : sg.xhm(); + } + + // The absences a group predicts over a fixed probe of low-index reflections - everything Stage B + // judges a candidate on. Two settings with the same signature are one hypothesis written twice + // (an alternative-origin entry, or a screw whose absences the centering already accounts for), so + // only the first of them is worth scoring. Four indices is enough for every Sohncke setting in + // the table: taking the probe to six or eight admits exactly the same candidates. + std::vector AbsenceSignature(const gemmi::SpaceGroup& sg) { + const gemmi::GroupOps gops = sg.operations(); + std::vector out; + out.reserve(9 * 9 * 9); + for (int h = -4; h <= 4; ++h) + for (int k = -4; k <= 4; ++k) + for (int l = -4; l <= 4; ++l) + out.push_back(gops.is_systematically_absent({{h, k, l}}) ? 1 : 0); + return out; + } + // The rotation part of a space group in the reference setting (identity included), as a // sorted list of matrices - the key that groups space groups into a candidate point group. // It must be the rotation SET, not gemmi's PointGroup enum: P321 and P312 are both "32" yet @@ -190,45 +241,66 @@ namespace { RotationSet rotation_set; std::vector rotations; // non-identity proper rotations const gemmi::SpaceGroup* representative = nullptr; + // This rung exists only because no reference setting names it (the a- and c-unique monoclinic + // 2-folds). It may be ADOPTED like any other, but it must not take part in judging a higher + // promotion - see the two sites that read it below. + bool widened = false; }; // Enumerate candidate point groups. When a holohedry is given (from the lattice metric), keep // only its subgroups - this both skips operators the lattice forbids and avoids accepting a // coincidental higher symmetry; all subgroups down to P1 are still candidates. - std::vector EnumeratePointGroups(const std::optional& holohedry) { + std::vector EnumeratePointGroups(const std::optional& holohedry, + const std::optional& cell, + bool all_rotation_sets) { std::vector out; std::map index; - for (const auto& sg : gemmi::spacegroup_tables::main) { - if (!sg.is_sohncke() || !sg.is_reference_setting()) - continue; + // Two passes. The reference settings first, exactly as before; then, only if asked and only + // for rotation sets NO reference setting carries, a setting that names them. The four such + // sets are the a-unique and c-unique monoclinic 2-folds and the two rhombohedral-axes trigonal + // groups, and without them a crystal whose only 2-fold lies on a or c has no rung to stand on + // between P1 and 222, so it falls to P1. Restricting the second pass to sets the first did not + // reach makes it a pure addition: every point group reachable before is still reached, by the + // same group, in the same setting. + for (int pass = 0; pass < 2; ++pass) { + if (pass == 1 && (!all_rotation_sets || !cell.has_value())) + break; + for (const auto& sg : gemmi::spacegroup_tables::main) { + if (!sg.is_sohncke() || sg.is_reference_setting() != (pass == 0)) + continue; - RotationSet rs = RotationSetOf(sg); - if (holohedry.has_value() && - !std::includes(holohedry->begin(), holohedry->end(), rs.begin(), rs.end())) - continue; + RotationSet rs = RotationSetOf(sg); + if (holohedry.has_value() && + !std::includes(holohedry->begin(), holohedry->end(), rs.begin(), rs.end())) + continue; + if (pass == 1 && + (index.count(rs) > 0 || !CellHostsRotations(*cell, sg.operations()))) + continue; - auto it = index.find(rs); - size_t pos; - if (it == index.end()) { - PointGroupInfo info; - for (const auto& op : sg.operations().derive_symmorphic().sym_ops) { - if (op.rot == gemmi::Op::identity().rot) - continue; - info.rotations.push_back(gemmi::Op{op.rot, {0, 0, 0}, op.notation}); + auto it = index.find(rs); + size_t pos; + if (it == index.end()) { + PointGroupInfo info; + for (const auto& op : sg.operations().derive_symmorphic().sym_ops) { + if (op.rot == gemmi::Op::identity().rot) + continue; + info.rotations.push_back(gemmi::Op{op.rot, {0, 0, 0}, op.notation}); + } + info.rotation_set = rs; + info.widened = (pass == 1); + pos = out.size(); + index[rs] = pos; + out.push_back(std::move(info)); + } else { + pos = it->second; } - info.rotation_set = rs; - pos = out.size(); - index[rs] = pos; - out.push_back(std::move(info)); - } else { - pos = it->second; + // Prefer a symmorphic representative (the plain point-group setting). + auto& info = out[pos]; + if (info.representative == nullptr || + (!info.representative->is_symmorphic() && sg.is_symmorphic())) + info.representative = &sg; } - // Prefer a symmorphic representative (the plain point-group setting). - auto& info = out[pos]; - if (info.representative == nullptr || - (!info.representative->is_symmorphic() && sg.is_symmorphic())) - info.representative = &sg; } return out; } @@ -479,7 +551,8 @@ SearchSpaceGroupResult SearchSpaceGroup( std::optional holohedry; if (opt.lattice_system.has_value()) holohedry = HolohedryRotationSet(opt.lattice_system.value()); - const auto point_groups = EnumeratePointGroups(holohedry); + const auto point_groups = + EnumeratePointGroups(holohedry, opt.cell, opt.enumerate_all_rotation_sets); // Every operator the search can ask about comes from this list, and scoring one is a pass over the // whole merge with a hash lookup per reflection - the most expensive thing in here. They do not @@ -649,7 +722,12 @@ SearchSpaceGroupResult SearchSpaceGroup( const PointGroupInfo& pg = *confirmed[i]; pg_cands.push_back({&pg, static_cast(pg.rotations.size()) + 1, confirmed_cc[i], cand_chi2[i], cand_b[i]}); - if (!pg.rotations.empty() && std::isfinite(cand_chi2[i])) + // A rung the widened enumeration added does not set the reference chi^2. Otherwise adding it + // would make every HIGHER promotion harder to reach - chi2_ref is a minimum - and the search + // could answer LOWER because a candidate was offered, which is the opposite of what offering + // it is for. Measured before this line existed: an F-cubic crystal's 432 and an orthorhombic + // crystal's 222 were both refused once the a- and c-unique 2-folds joined the list. + if (!pg.rotations.empty() && std::isfinite(cand_chi2[i]) && !pg.widened) chi2_ref = std::min(chi2_ref, cand_chi2[i]); } @@ -690,7 +768,11 @@ SearchSpaceGroupResult SearchSpaceGroup( if (!c.pg->rotations.empty()) { int parent_order = 0; for (const auto& s : pg_cands) - if (s.order < c.order && s.order >= parent_order + // Same reason as chi2_ref above: a rung only the widened enumeration offers is not a + // parent. The promotion answers to the most damning of its parents, so admitting two + // more order-2 subgroups of 222 makes the 222 step strictly harder than it was before + // the rung was offered at all. + if (!s.pg->widened && s.order < c.order && s.order >= parent_order && std::includes(c.pg->rotation_set.begin(), c.pg->rotation_set.end(), s.pg->rotation_set.begin(), s.pg->rotation_set.end())) { if (s.order > parent_order) { @@ -874,6 +956,30 @@ SearchSpaceGroupResult SearchSpaceGroup( for (const auto& sg : gemmi::spacegroup_tables::main) if (sg.is_sohncke() && sg.is_reference_setting() && RotationSetOf(sg) == best_pg->rotation_set) sg_cands.push_back(&sg); + // The non-reference SETTINGS of the same point group, when asked for - and unconditionally when + // the point group is one only a non-reference setting carries (Stage A's second pass), since + // otherwise there is no candidate at all and the group would be lost after being found. These are + // alternative namings at the same order, so nothing here can promote the point group; what they + // add is a screw or a centering on the axis the data show it on. Two refusals bound them: a + // setting whose axes the cell does not have is not offered, and one predicting exactly the + // absences a candidate already offered predicts is the same hypothesis under another name. + if (opt.cell.has_value() && + (opt.enumerate_all_settings || (opt.enumerate_all_rotation_sets && sg_cands.empty()))) { + std::vector> signatures; + for (const auto* c : sg_cands) + signatures.push_back(AbsenceSignature(*c)); + for (const auto& sg : gemmi::spacegroup_tables::main) { + if (!sg.is_sohncke() || sg.is_reference_setting() || + RotationSetOf(sg) != best_pg->rotation_set || + !CellHostsRotations(*opt.cell, sg.operations())) + continue; + auto sig = AbsenceSignature(sg); + if (std::find(signatures.begin(), signatures.end(), sig) != signatures.end()) + continue; + signatures.push_back(std::move(sig)); + sg_cands.push_back(&sg); + } + } std::vector sg_scored(sg_cands.size()); ParallelFor(static_cast(sg_cands.size()), std::min(opt.nthreads, sg_cands.size()), [&](int ci) { @@ -1115,8 +1221,15 @@ SearchSpaceGroupResult SearchSpaceGroup( // count measures the sweep's geometry, not the evidence: six axial reflections measured at zero // against a row that averages 1.4x the shell mean settle the question, while twenty uniformly weak // ones settle nothing. So the screw class is judged by its zones' AbsenceEvidence instead. + // A candidate in a non-reference setting must also have had its CENTERING tested here. The + // reference-setting path can adopt an untested centering because something else backs it - the + // caller's centred-lattice re-test reindexes the data into the metric candidate's conventional + // cell and only then commits - and a non-reference setting has no such backing: adopting one on + // an absence class this merge does not contain would name a lattice the data never showed, which + // is worse than the group it displaces. Refused, not warned about. auto eligible = [&](const SpaceGroupCandidateScore& s) { return s.consistent + && !(s.centering_untested && !s.space_group.is_reference_setting()) && (s.centering_absent == 0 || s.centering_absent >= opt.min_absent_observed) && (s.screw_absent == 0 || s.screw_absence_evidence >= opt.min_screw_absence_evidence); }; @@ -1199,7 +1312,7 @@ std::string SearchSpaceGroupResultToText(const SearchSpaceGroupResult& result, for (size_t i = 0; i < count; ++i) { const auto& c = result.candidates[i]; os << (c.selected ? "* " : " ") - << std::setw(10) << std::left << c.space_group.short_name() << std::right + << std::setw(10) << std::left << SettingName(c.space_group) << std::right << std::setw(9) << c.absent_observed << std::setw(7) << c.absent_violations << std::setw(11) << std::fixed << std::setprecision(2) << c.absent_mean_i_over_sigma << std::setw(11) << std::fixed << std::setprecision(2) << c.present_mean_i_over_sigma @@ -1229,7 +1342,7 @@ std::string SearchSpaceGroupResultToText(const SearchSpaceGroupResult& result, for (const auto& c : result.candidates) { if (!c.selected || c.screw_zones.empty()) continue; - os << "Screw conditions of " << c.space_group.short_name() << ", zone by zone:\n"; + os << "Screw conditions of " << SettingName(c.space_group) << ", zone by zone:\n"; for (const auto& z : c.screw_zones) { os << " " << RowLabel(z.row) << " " << std::setw(3) << z.n_absent << " absent "; if (z.n_control == 0) @@ -1242,9 +1355,9 @@ std::string SearchSpaceGroupResultToText(const SearchSpaceGroupResult& result, } if (result.best_space_group.has_value()) { - os << "Best space group: " << result.best_space_group->short_name(); + os << "Best space group: " << SettingName(*result.best_space_group); for (const auto& alt : result.alternatives) - os << " or " << alt.short_name(); + os << " or " << SettingName(alt); if (!result.alternatives.empty()) os << " (indistinguishable from these data)"; os << "\n"; @@ -1260,11 +1373,11 @@ std::string SearchSpaceGroupResultToText(const SearchSpaceGroupResult& result, continue; if (!said_setting) { os << " These are NOT all the same lattice, and one cell is reported for this run -\n" - " " << result.best_space_group->short_name() << "'s (centering " << sel_centring + " " << SettingName(*result.best_space_group) << "'s (centering " << sel_centring << "). These need a cell of their own:\n"; said_setting = true; } - os << " " << std::setw(10) << std::left << alt.short_name() << std::right + os << " " << std::setw(10) << std::left << SettingName(alt) << std::right << "centering " << alt.centring_type() << ", cell volume " << std::fixed << std::setprecision(2) << static_cast(gemmi::centring_vectors(alt.centring_type()).size()) diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index b78b7e42e..45fe7a83b 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -11,6 +11,7 @@ #include "../../common/Reflection.h" #include "gemmi/symmetry.hpp" +#include "gemmi/unitcell.hpp" // Determine the likely space group of a dataset from its P1-merged intensities, in the spirit // of POINTLESS (Evans 2006): @@ -94,6 +95,29 @@ struct SearchSpaceGroupOptions { // in the systematic absences. Stage B therefore tests every centering allowed by the point group // and confirms it from the data (h+k+l etc. absent), rather than trusting a geometric hint. + // The conventional cell the merged reflections are indexed on. Only needed with the two + // enumeration options below: a setting names the symmetry directions by AXIS, so a candidate in a + // setting the cell does not support (a c-unique monoclinic group against a b-unique cell) is + // refused before it is scored. Unset with either option on means no non-reference setting is + // offered at all - the check is a precondition, not an optional extra. + std::optional cell; + + // Enumerate the non-reference SETTINGS of the chosen point group in Stage B - P 2_1 2 2 and + // P 2 2_1 2 beside P 2 2 2_1, A 2 2 2 beside C 2 2 2, I 1 2 1 beside C 1 2 1. These are + // alternative namings at the same group ORDER, so this cannot promote the point group; what it + // adds is the ability to put a screw or a centering on the axis the data actually show it on. + // Without it a crystal whose screws are on b and c is reported as P 2 2 2_1 - the wrong group, + // not a lower one, because the only candidate that predicts the 00l absences and nothing else + // wins on zero extra evidence. + bool enumerate_all_settings = false; + + // Enumerate the rotation sets that only a non-reference setting carries in Stage A - the a-unique + // and c-unique monoclinic 2-folds. Unlike the option above this DOES add promotions: a crystal + // whose only 2-fold lies on a or c currently has no rung to stand on and falls to P1. Requires a + // cell, and requires enumerate_all_settings for those point groups to be nameable in Stage B + // (no reference setting carries them), which this option turns on for them by itself. + bool enumerate_all_rotation_sets = false; + // Friedel mates are treated as equivalent when matching HKLs (i.e. anomalous signal ignored). bool merge_friedel = true; diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index cb2705db3..94a8419d2 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -139,8 +139,12 @@ std::string RenderResultReport(const std::string &output_prefix, Key(os, "LATTICE_FOUND", (result.consensus_cell.has_value() ? "TRUE" : "FALSE")); if (result.consensus_cell.has_value()) Key(os, "UNIT_CELL_CONSTANTS", CellString(*result.consensus_cell)); - if (result.space_group.has_value()) + if (result.space_group.has_value()) { Key(os, "SPACE_GROUP_NUMBER", result.space_group->number); + // The number alone does not name the setting, and the search can now return one that is not + // gemmi's reference setting - P 1 1 2_1 and P 1 2_1 1 are both number 4, on different axes. + Key(os, "SPACE_GROUP_NAME", result.space_group->xhm()); + } if (result.indexing_rate.value_or(0.0f) <= 0.0f) warnings.emplace_back("No image indexed - no crystal lattice was determined from this dataset"); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index b516d7e7d..a772114d5 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -3385,6 +3385,13 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // indexing. Centering is not constrained here - it is determined from the absences. if (end_msg.rotation_lattice_type.has_value()) sg_opts.lattice_system = end_msg.rotation_lattice_type->crystal_system; + // The cell the reflections are indexed on, which is what lets the search offer a candidate + // in a setting other than gemmi's reference one: a setting names the symmetry directions by + // axis, so a candidate the cell cannot host is refused before it is scored. + if (result.consensus_cell.has_value()) + sg_opts.cell = gemmi::UnitCell(*result.consensus_cell); + sg_opts.enumerate_all_settings = true; + sg_opts.enumerate_all_rotation_sets = true; // The error model this merge was scaled with, so the search's "genuinely present" cut can be // read on the counting sigma instead of on the merged one (SearchSpaceGroupOptions::merge_isa). // Set per merge, not once: the two arms below are scaled separately and their ISa differ by up @@ -3432,9 +3439,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // by value and the two become the same result. // Which point group, by its symmorphic representative rather than by name - gemmi calls // both P321 and P312 "32", so two arms that disagree about WHICH 2-folds are real read - // as agreeing. 0 when no point group was chosen. + // as agreeing. By the representative's full symbol rather than its NUMBER, for the same + // reason one level down: P 1 2 1 and P 1 1 2 are both number 3 and put their 2-fold on + // different axes. Empty when no point group was chosen. const auto pg_of = [](const SearchSpaceGroupResult &r) { - return r.point_group_representative.has_value() ? r.point_group_representative->number : 0; + return r.point_group_representative.has_value() ? r.point_group_representative->xhm() + : std::string(); }; if (order_of(alt) > 0 && order_of(alt) == order_of(sg_search) && pg_of(alt) != pg_of(sg_search)) { @@ -3540,6 +3550,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b reindex_hkl(m, reindex); SearchSpaceGroupOptions o2 = sg_opts; o2.lattice_system = cand.system; + o2.cell = gemmi::UnitCell(cand.conventional.GetUnitCell()); const auto s2 = SearchSpaceGroup(merged_c, o2); if (s2.best_space_group.has_value() && s2.best_space_group->number > 1) { // The centring cannot be confirmed from absences here (integrated in the primitive @@ -3575,14 +3586,19 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // reflections it extinguishes were never measured. (The report warns // about it too, from the candidate's centering_untested flag, which // survives into the pass that is written.) - logger.Warning("The {} centering of {} comes from the lattice METRIC, not " - "from the intensities: the crystal is indexed on the " - "primitive sub-cell, so the reflections a {}-centred " - "lattice extinguishes were never predicted or integrated. " - "The point group ({}) IS confirmed, by the operator " - "correlations.", - cand.centering, chosen.short_name(), cand.centering, - s2.point_group_hm); + // Only for a CENTRED candidate: a primitive one extinguishes nothing, so + // there is no untested absence class to warn about. Before the enumeration + // was widened this branch could only ever adopt a centred group, so the + // condition was implicit. + if (cand.centering != 'P') + logger.Warning("The {} centering of {} comes from the lattice METRIC, " + "not from the intensities: the crystal is indexed on " + "the primitive sub-cell, so the reflections a " + "{}-centred lattice extinguishes were never predicted " + "or integrated. The point group ({}) IS confirmed, by " + "the operator correlations.", + cand.centering, chosen.short_name(), cand.centering, + s2.point_group_hm); } } } @@ -3606,7 +3622,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // a later change could tighten without anyone noticing what else was leaning on them. { const gemmi::SpaceGroup sg = sg_search.best_space_group.value_or(gemmi::get_spacegroup_p1()); - logger.Info("Adopting space group {} (number {})", sg.short_name(), sg.number); + // xhm(), not short_name(): the search can now return a non-reference setting, and + // short_name() cannot tell P 1 2 1 from P 1 1 2. + logger.Info("Adopting space group {} (number {})", sg.xhm(), sg.number); // A reindex is committed only when the intensity re-test confirmed a higher symmetry in a // centred setting: bring the integrated reflections + cell into that setting and re-ingest // the rotation merge, so the final scaling/merging folds the equivalents correctly. diff --git a/tests/SearchSpaceGroupTest.cpp b/tests/SearchSpaceGroupTest.cpp index 0426e9649..b9cd6a678 100644 --- a/tests/SearchSpaceGroupTest.cpp +++ b/tests/SearchSpaceGroupTest.cpp @@ -367,3 +367,126 @@ TEST_CASE("SearchSpaceGroup operator correlation reads symmetry, not the resolut CHECK(std::fabs(two_fold_cc(0.0) - two_fold_cc(6.0)) < 0.05); } } + +// The enumeration reaches the settings gemmi does not call the reference one, and only when the cell +// has the axes they name. Both halves are pinned here: `P 1 1 2_1` puts its 2-fold and its screw on +// c, which no reference setting can express (Stage A never offers the rotation and Stage B never +// offers the group), so without the two options the answer is P1; with them, and with a cell whose +// unique axis IS c, it is named; and with a cell whose unique axis is b the same candidate is +// refused rather than adopted on axes the crystal does not have. +TEST_CASE("SearchSpaceGroup names a non-reference setting only on a cell that hosts it") { + const gemmi::SpaceGroup& sg = gemmi::get_spacegroup_by_name("P 1 1 21"); + const auto merged = GenerateMergedReflectionsForSpaceGroup(sg, 12); + + SearchSpaceGroupOptions opt; + opt.merge_friedel = true; + + SECTION("narrow enumeration cannot name it") { + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->number == 1); + } + + SECTION("widened enumeration names it on a c-unique cell") { + opt.cell = gemmi::UnitCell(40.0, 50.0, 60.0, 90.0, 90.0, 100.0); + opt.enumerate_all_settings = true; + opt.enumerate_all_rotation_sets = true; + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->xhm() == "P 1 1 21"); + } + + SECTION("a b-unique cell refuses it") { + opt.cell = gemmi::UnitCell(40.0, 50.0, 60.0, 90.0, 100.0, 90.0); + opt.enumerate_all_settings = true; + opt.enumerate_all_rotation_sets = true; + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->number == 1); + } +} + +// The screw axes of an orthorhombic crystal can lie on any pair of axes, and only one of the three +// namings of #18 is a reference setting. With the narrow enumeration the group that predicts a +// SUBSET of the real absences and nothing else wins on no evidence at all, so the reported group is +// wrong rather than low - the widening is what makes the correct one available. +TEST_CASE("SearchSpaceGroup names an orthorhombic screw pair on the axes it lies on") { + const gemmi::SpaceGroup& sg = gemmi::get_spacegroup_by_name("P 2 21 21"); + const auto merged = GenerateMergedReflectionsForSpaceGroup(sg, 14); + + SearchSpaceGroupOptions opt; + opt.merge_friedel = true; + opt.lattice_system = gemmi::CrystalSystem::Orthorhombic; + + SECTION("narrow enumeration reports the wrong group") { + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->number != 18); + } + + SECTION("widened enumeration reports it") { + opt.cell = gemmi::UnitCell(40.0, 50.0, 60.0, 90.0, 90.0, 90.0); + opt.enumerate_all_settings = true; + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->xhm() == "P 2 21 21"); + } +} + +// The centring half of the same widening. A, B and C centring on one orthorhombic cell are three +// different lattices, and only C is a reference setting, so an A-centred crystal used to have its +// centring refused (its absent class is not the one C predicts) and came out primitive. The +// candidate is now offered, and it has to be adopted from its own absences rather than from the +// metric, which cannot tell A from C at all. +TEST_CASE("SearchSpaceGroup names an A-centred orthorhombic lattice") { + const gemmi::SpaceGroup& sg = gemmi::get_spacegroup_by_name("A 2 2 2"); + const auto merged = GenerateMergedReflectionsForSpaceGroup(sg, 12); + + SearchSpaceGroupOptions opt; + opt.merge_friedel = true; + opt.lattice_system = gemmi::CrystalSystem::Orthorhombic; + + SECTION("narrow enumeration cannot name it") { + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->centring_type() != 'A'); + } + + SECTION("widened enumeration names it") { + opt.cell = gemmi::UnitCell(40.0, 50.0, 60.0, 90.0, 90.0, 90.0); + opt.enumerate_all_settings = true; + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->xhm() == "A 2 2 2"); + } +} + +// The null the widening has to survive. Stage A's second pass offers the a- and c-unique 2-folds on +// any metric that could host them, which is every orthorhombic one - so a genuinely triclinic +// crystal sitting on a pseudo-orthorhombic cell is now offered three promotions where it used to be +// offered one. It must still be refused all three: the added candidates go through the same operator +// correlation as every other, and a rotation the intensities do not have scores nothing. +TEST_CASE("SearchSpaceGroup does not promote triclinic data on a pseudo-orthorhombic cell") { + const gemmi::SpaceGroup& sg = gemmi::get_spacegroup_by_name("P 1"); + const auto merged = GenerateMergedReflectionsForSpaceGroup(sg, 10); + + SearchSpaceGroupOptions opt; + opt.merge_friedel = true; + opt.lattice_system = gemmi::CrystalSystem::Orthorhombic; + opt.cell = gemmi::UnitCell(40.0, 50.0, 60.0, 90.0, 90.0, 90.0); + opt.enumerate_all_settings = true; + opt.enumerate_all_rotation_sets = true; + + const auto result = SearchSpaceGroup(merged, opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.best_space_group.has_value()); + CHECK(result.best_space_group->number == 1); + CHECK(result.point_group_order == 1); +}