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); +}