diff --git a/common/DatasetSettings.cpp b/common/DatasetSettings.cpp index 36b6e2ee4..2bcf4167e 100644 --- a/common/DatasetSettings.cpp +++ b/common/DatasetSettings.cpp @@ -133,7 +133,16 @@ DatasetSettings &DatasetSettings::SpaceGroupNumber(std::optional input) check_min("Space group number", input, 1); check_max("Space group number", input, 230); } - space_group_number = input; + // A bare number can only name the reference setting. + if (input) + space_group = *gemmi::find_spacegroup_by_number(static_cast(*input)); + else + space_group.reset(); + return *this; +} + +DatasetSettings &DatasetSettings::SetSpaceGroup(const std::optional &input) { + space_group = input; return *this; } @@ -205,7 +214,13 @@ std::optional DatasetSettings::GetUnitCell() const { } std::optional DatasetSettings::GetSpaceGroupNumber() const { - return space_group_number; + if (!space_group) + return std::nullopt; + return space_group->number; +} + +const std::optional &DatasetSettings::GetSpaceGroup() const { + return space_group; } std::string DatasetSettings::GetSampleName() const { diff --git a/common/DatasetSettings.h b/common/DatasetSettings.h index 09164076e..0b5980507 100644 --- a/common/DatasetSettings.h +++ b/common/DatasetSettings.h @@ -13,6 +13,8 @@ #include "JFJochMessages.h" #include "XrayFluorescenceSpectrum.h" +#include "gemmi/symmetry.hpp" + class DatasetSettings { int64_t images_per_trigger; int64_t ntrigger; @@ -29,7 +31,9 @@ class DatasetSettings { std::string sample_name; std::optional unit_cell; - std::optional space_group_number; + // The whole group, not its number: a number can only name the reference setting, so + // "P 1 1 2" and "R 3:R" are lost the moment a number is all that is kept. + std::optional space_group; std::optional total_flux; std::optional attenuator_transmission; @@ -82,6 +86,7 @@ public: DatasetSettings& Compression(CompressionAlgorithm input); DatasetSettings& SetUnitCell(const std::optional &cell); DatasetSettings& SpaceGroupNumber(std::optional input); + DatasetSettings& SetSpaceGroup(const std::optional &input); DatasetSettings& SampleName(std::string input); DatasetSettings& AttenuatorTransmission(const std::optional &input); DatasetSettings& TotalFlux(const std::optional &input); @@ -126,6 +131,7 @@ public: const nlohmann::json& GetImageAppendix() const; std::optional GetUnitCell() const; std::optional GetSpaceGroupNumber() const; + const std::optional &GetSpaceGroup() const; std::string GetSampleName() const; float GetPhotonEnergy_keV() const; float GetBeamX_pxl() const; diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index 719f28e31..fa535919b 100644 --- a/common/DiffractionExperiment.cpp +++ b/common/DiffractionExperiment.cpp @@ -198,6 +198,11 @@ DiffractionExperiment &DiffractionExperiment::SpaceGroupNumber(std::optional &input) { + dataset.SetSpaceGroup(input); + return *this; +} + DiffractionExperiment &DiffractionExperiment::StorageCells(int64_t input) { detector_settings.StorageCells(input); return *this; @@ -1792,14 +1797,13 @@ std::string DiffractionExperiment::GetSpaceGroupName() const { return sg->short_name(); } -std::optional DiffractionExperiment::GetGemmiSpaceGroup() const { - auto sg = GetSpaceGroupNumber(); - if (!sg) - return std::nullopt; - const gemmi::SpaceGroup *g = gemmi::find_spacegroup_by_number(sg.value()); - if (g == nullptr) - return std::nullopt; - return *g; +const std::optional &DiffractionExperiment::GetGemmiSpaceGroup() const { + return dataset.GetSpaceGroup(); +} + +const gemmi::SpaceGroup &DiffractionExperiment::GetSpaceGroupOrP1() const { + const auto &sg = dataset.GetSpaceGroup(); + return sg ? *sg : *gemmi::find_spacegroup_by_number(1); } char DiffractionExperiment::GetCentering() const { diff --git a/common/DiffractionExperiment.h b/common/DiffractionExperiment.h index 84d35205f..83f7cc350 100644 --- a/common/DiffractionExperiment.h +++ b/common/DiffractionExperiment.h @@ -140,6 +140,7 @@ public: DiffractionExperiment& Compression(CompressionAlgorithm input); DiffractionExperiment& SetUnitCell(const std::optional &cell); DiffractionExperiment& SpaceGroupNumber(std::optional input); + DiffractionExperiment& SetSpaceGroup(const std::optional &input); DiffractionExperiment& SampleName(const std::string &input); DiffractionExperiment& AttenuatorTransmission(const std::optional &input); DiffractionExperiment& TotalFlux(const std::optional &input); @@ -434,7 +435,10 @@ public: float GetPoniRot2_rad() const; float GetPoniRot3_rad() const; - std::optional GetGemmiSpaceGroup() const; + const std::optional &GetGemmiSpaceGroup() const; + // The adopted space group, or P1 when none has been determined - the convention every + // scaling and merging step uses. + const gemmi::SpaceGroup &GetSpaceGroupOrP1() const; gemmi::CrystalSystem GetCrystalSystem() const; std::string GetSpaceGroupName() const; char GetCentering() const; diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index b8e0fd0d8..5e3360708 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -373,7 +373,11 @@ struct EndMessage { std::vector rotation_extra_lattices; std::optional unit_cell; // Space group determined by the offline analysis (overrides the start message when writing the - // master, since it is only known after merging). + // master, since it is only known after merging). Both spellings are carried and a reader prefers + // the name: only the name keeps the setting - "P 1 1 2" and "R 3:R" both come back from a number + // as the reference setting - while the number is what a reader written before the name existed + // understands. Set them together; a name without its number is a half-filled message. + std::optional space_group_name; std::optional space_group_number; // Change of basis (3x3 integers, row major) from the setting the PER-IMAGE reflections and // lattices were written in to the setting of unit_cell / rotation_lattice above: diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 17b5c07d8..4fd9f4b82 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -55,6 +55,7 @@ * `rugnux --mode calibration --no-refine-tilt`, and the "Refine detector tilt" tick box in the viewer's calibration panel, fit the beam centre and the distance with the detector tilt held fixed, for a calibration meant for a program that cannot express a tilted detector. * In `jfjoch_viewer`, Alt and the mouse wheel step through the dataset one image at a time. * The API description of `fft_high_resolution_A` says what the setting does: it sizes the FFT's projection histogram and does not filter spots by resolution. +* The space group is carried as its full name rather than its number - through the analysis, in `-S`, and in `/entry/sample/space_group` - so a setting a number cannot name (`P 1 1 2`, `R 3:R`) is no longer lost. ### 1.0.0-rc.165 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. diff --git a/docs/HDF5.md b/docs/HDF5.md index db3a8f6b7..c4cee0d33 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -242,6 +242,8 @@ beam); for a detector this system assembled itself they are the McStas form of t | `temperature` | NXmx | K | | `transformations/` (NXtransformations) | NXmx | the sample chain, written in mounting order; hard-linked as `/entry/sample/goniometer` | | `unit_cell` | base | `[a, b, c, α, β, γ]` | +| `space_group_number` | base | International Tables number | +| `space_group` | base | extended Hermann-Mauguin name, e.g. `P 43 21 2`, `R 3:H` — this is the field that carries the **setting**, and the one the reader takes the group from; the number alone always reads back as the reference setting | | `ub_matrix` | base | `[1, 3, 3]`, Angstrom⁻¹ | The chain is written from the base outwards, so the innermost axis — the one `depends_on` names — is diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index 5b1605b9f..d6e8a753e 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -1519,6 +1519,10 @@ namespace { GetCBORUInt8Array(value, message.niggli_class); else if (key == "unit_cell") message.unit_cell = ProcessUnitCellElement(value); + else if (key == "space_group_name") + message.space_group_name = GetCBORString(value); + else if (key == "space_group_number") + message.space_group_number = GetCBORUInt(value); else if (key == "pixel_sum") GetCBORInt64Array(value, message.pixel_sum); else if (key == "rotation_lattice") { diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index a80c97937..a37c135df 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -833,6 +833,10 @@ void CBORStream2Serializer::SerializeSequenceEnd(const EndMessage& message) { CBOR_ENC(mapEncoder, "niggli_class", message.niggli_class); CBOR_ENC(mapEncoder, "pixel_sum", message.pixel_sum); CBOR_ENC(mapEncoder, "unit_cell", message.unit_cell); + // Both spellings of the determined group. The name carries the setting and is what a reader + // should take; the number is beside it for a reader written before the name existed. + CBOR_ENC(mapEncoder, "space_group_name", message.space_group_name); + CBOR_ENC(mapEncoder, "space_group_number", message.space_group_number); cborErr(cbor_encoder_close_container(&encoder, &mapEncoder)); diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index df82fc39f..cb0a75a7f 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -101,7 +101,7 @@ std::vector BuildMergedRows(const std::vector &r // Anomalous: group the two mates by their (shared) Friedel-merged ASU representative, then form // IMEAN / F as their inverse-variance Friedel mean. A single generator gives both the group key // (its hkl, identical for +hkl and -hkl) and which mate this row is (.plus). - const HKLKeyGenerator key_gen(false, experiment.GetSpaceGroupNumber().value_or(1)); + const HKLKeyGenerator key_gen(false, experiment.GetSpaceGroupOrP1()); struct AnomRow { int h = 0, k = 0, l = 0; float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN; @@ -186,13 +186,11 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_cell.angle_beta " << Fmt(unitCell.beta, 2) << "\n"; out << "_cell.angle_gamma " << Fmt(unitCell.gamma, 2) << "\n"; - auto *sg = gemmi::find_spacegroup_by_number(experiment.GetSpaceGroupNumber().value_or(1)); - if (sg == nullptr) - throw std::runtime_error("WriteMmcifReflections: invalid space group number"); + const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1(); // ---------- _symmetry ---------- - out << "_symmetry.space_group_name_H-M " << CifStr(sg->hm) << "\n"; - out << "_symmetry.Int_Tables_number " << sg->number << "\n"; + out << "_symmetry.space_group_name_H-M " << CifStr(sg.hm) << "\n"; + out << "_symmetry.Int_Tables_number " << sg.number << "\n"; out << "#\n"; // ---------- _diffrn_source / _diffrn_detector ---------- @@ -439,8 +437,7 @@ void WriteMtzReflections(const std::vector &reflections, gemmi::Mtz mtz; // Optional but recommended metadata - mtz.spacegroup = gemmi::find_spacegroup_by_number( - experiment.GetSpaceGroupNumber().value_or(1)); + mtz.spacegroup = &experiment.GetSpaceGroupOrP1(); mtz.set_cell_for_all(unitCell); // Producing-software provenance in the MTZ header (title + HISTORY, the CCP4 convention). @@ -698,7 +695,7 @@ void WriteUnmergedMtzReflections(const std::vector &outcomes bool sum_partials, const std::string &filename) { gemmi::Mtz mtz; - mtz.spacegroup = gemmi::find_spacegroup_by_number(experiment.GetSpaceGroupNumber().value_or(1)); + mtz.spacegroup = &experiment.GetSpaceGroupOrP1(); mtz.set_cell_for_all(unitCell); mtz.title = "Rugnux unmerged observations"; mtz.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction"); diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 0b89a67f5..c2ac868b8 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -522,20 +522,11 @@ bool RotationIndexer::AccumulationFull() const { void RotationIndexer::ForceLattice(const CrystalLattice &lattice) { indexed_lattice = lattice; - auto sg_num = experiment.GetSpaceGroupNumber().value_or(1); - auto sg = gemmi::find_spacegroup_by_number(sg_num); - if (sg != nullptr) { - search_result_ = LatticeSearchResult{ - .niggli_class = 0, // Since Niggli class was not searched for, we don't know which one - .conventional = lattice, // If lattice provided, it is for now primitive == conventional - .system = sg->crystal_system(), - .centering = sg->centring_type(), - }; - } else - search_result_ = LatticeSearchResult{ - .niggli_class = 0, // Since Niggli class was not searched for, we don't know which one - .conventional = lattice, // If lattice provided, it is for now primitive == conventional - .system = gemmi::CrystalSystem::Triclinic, - .centering = 'P', - }; + const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1(); + search_result_ = LatticeSearchResult{ + .niggli_class = 0, // Since Niggli class was not searched for, we don't know which one + .conventional = lattice, // If lattice provided, it is for now primitive == conventional + .system = sg.crystal_system(), + .centering = sg.centring_type(), + }; } diff --git a/image_analysis/scale_merge/FrenchWilson.cpp b/image_analysis/scale_merge/FrenchWilson.cpp index 35d222640..8b32ff824 100644 --- a/image_analysis/scale_merge/FrenchWilson.cpp +++ b/image_analysis/scale_merge/FrenchWilson.cpp @@ -68,7 +68,7 @@ Posterior integrate_posterior(double I, double sigma, double sigma_wilson, bool } // namespace -void ApplyFrenchWilson(std::vector &merged, int32_t space_group_number, +void ApplyFrenchWilson(std::vector &merged, const gemmi::SpaceGroup &space_group, const FrenchWilsonOptions &opts) { // Naive amplitude sqrt(max(I,0)) for a missing / strong / untrusted intensity; NaN in -> NaN out // (a missing Bijvoet hand stays missing). Fills one (F, sigmaF) pair. @@ -85,12 +85,9 @@ void ApplyFrenchWilson(std::vector &merged, int32_t space_grou naive_one(r.I_minus, r.sigma_minus, r.F_minus, r.sigmaF_minus); }; - const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); - if (sg == nullptr || merged.empty()) { - for (auto &r : merged) naive_all(r); + if (merged.empty()) return; - } - const gemmi::GroupOps gops = sg->operations(); + const gemmi::GroupOps gops = space_group.operations(); float d_min = std::numeric_limits::max(), d_max = 0.0f; for (const auto &r : merged) diff --git a/image_analysis/scale_merge/FrenchWilson.h b/image_analysis/scale_merge/FrenchWilson.h index fd41e020b..f8c1e390a 100644 --- a/image_analysis/scale_merge/FrenchWilson.h +++ b/image_analysis/scale_merge/FrenchWilson.h @@ -7,6 +7,7 @@ #include #include "../../common/Reflection.h" +#include "gemmi/symmetry.hpp" struct FrenchWilsonOptions { int num_shells = 20; // resolution shells for the Wilson mean @@ -21,5 +22,5 @@ struct FrenchWilsonOptions { // prior uses the resolution-shell mean intensity, the correct centric/acentric form, and the // reflection's epsilon (symmetry-enhancement) multiplicity. Strong reflections reduce to sqrt(I); // reflections with an unusable I/sigma fall back to sqrt(max(I,0)) with propagated sigma. -void ApplyFrenchWilson(std::vector &merged, int32_t space_group_number, +void ApplyFrenchWilson(std::vector &merged, const gemmi::SpaceGroup &space_group, const FrenchWilsonOptions &opts = {}); diff --git a/image_analysis/scale_merge/HKLKey.cpp b/image_analysis/scale_merge/HKLKey.cpp index ba60ab77c..ef088cf04 100644 --- a/image_analysis/scale_merge/HKLKey.cpp +++ b/image_analysis/scale_merge/HKLKey.cpp @@ -26,10 +26,6 @@ uint64_t HKLKey::pack() const { return (hh << 1) | (kk << (bits + 1)) | (ll << (2 * bits + 1)) | (plus ? 1ULL : 0ULL); } -HKLKeyGenerator::HKLKeyGenerator(bool merge_friedel, int32_t space_group_number) - : HKLKeyGenerator(merge_friedel, *gemmi::find_spacegroup_by_number(space_group_number)) { -} - HKLKeyGenerator::HKLKeyGenerator(bool merge_friedel, const gemmi::SpaceGroup &sg) : merge_friedel(merge_friedel), sg(sg), diff --git a/image_analysis/scale_merge/HKLKey.h b/image_analysis/scale_merge/HKLKey.h index 26fc8bad3..2fb74b4ba 100644 --- a/image_analysis/scale_merge/HKLKey.h +++ b/image_analysis/scale_merge/HKLKey.h @@ -34,7 +34,6 @@ class HKLKeyGenerator { gemmi::ReciprocalAsu asu; public: - HKLKeyGenerator(bool merge_friedel, int32_t space_group_number); HKLKeyGenerator(bool merge_friedel, const gemmi::SpaceGroup &sg); HKLKey operator()(const Reflection &r) const; HKLKey operator()(const MergedReflection &r) const; diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index d865b8843..b55e72f24 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -33,7 +33,7 @@ namespace { } MergeOnTheFly::MergeOnTheFly(const DiffractionExperiment &x) - : space_group_number(x.GetSpaceGroupNumber().value_or(1)), + : space_group(x.GetSpaceGroupOrP1()), scaling_settings(x.GetScalingSettings()), indexing_settings(x.GetIndexingSettings()), high_resolution_limit(scaling_settings.GetHighResolutionLimit_A()), @@ -45,7 +45,7 @@ MergeOnTheFly::MergeOnTheFly(const DiffractionExperiment &x) ? std::optional(scaling_settings.GetMinCCForImage()) : std::nullopt), min_partiality(scaling_settings.GetMinPartiality()), - generator(scaling_settings.GetMergeFriedel(), space_group_number), + generator(scaling_settings.GetMergeFriedel(), space_group), reject_outliers(scaling_settings.GetOutlierRejectNsigma() > 0.0), reject_nsigma(scaling_settings.GetOutlierRejectNsigma()) { } @@ -377,8 +377,8 @@ std::vector MergeOnTheFly::ExportReflections() { out.emplace_back(mr); } - AssignRfreeFlags(out, space_group_number, scaling_settings.GetRfreeFraction()); - ApplyFrenchWilson(out, space_group_number); + AssignRfreeFlags(out, space_group, scaling_settings.GetRfreeFraction()); + ApplyFrenchWilson(out, space_group); return out; } @@ -461,7 +461,7 @@ std::pair ImageReferenceCC(const std::vector &reflec return {cov / std::sqrt(var_x * var_y), n}; } -void CalcPossibleReflections(int space_group_number , +void CalcPossibleReflections(const gemmi::SpaceGroup &sg, const UnitCell &cell, double d_min, double d_max, @@ -469,14 +469,10 @@ void CalcPossibleReflections(int space_group_number , std::vector &acc, bool merge_friedel) { gemmi::UnitCell gemmi_cell = cell; - const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); - if (sg == nullptr) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Invalid space group number " + std::to_string(space_group_number)); // Generate unique reflections - std::vector possible_hkls = gemmi::make_miller_vector(gemmi_cell, sg, d_min, d_max, true); - const gemmi::GroupOps gops = sg->operations(); + std::vector possible_hkls = gemmi::make_miller_vector(gemmi_cell, &sg, d_min, d_max, true); + const gemmi::GroupOps gops = sg.operations(); CrystalLattice lattice(cell); const auto astar = lattice.Astar(); const auto bstar = lattice.Bstar(); @@ -553,7 +549,7 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m std::vector acc(n_shells); if (reference_cell.has_value()) - CalcPossibleReflections(space_group_number, reference_cell.value(), + CalcPossibleReflections(space_group, reference_cell.value(), d_min_pad, d_max_pad, shells, acc, scaling_settings.GetMergeFriedel()); CorrelationCoefficient cc_half_overall; diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 0d85450eb..48f6290b1 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -153,7 +153,7 @@ inline XdsErrorModel ToXdsErrorModel(double a, double b) { class MergeOnTheFly { mutable std::mutex merged_mutex; - const int space_group_number = 1; + const gemmi::SpaceGroup space_group; ScalingSettings scaling_settings; IndexingSettings indexing_settings; diff --git a/image_analysis/scale_merge/ReindexAmbiguity.cpp b/image_analysis/scale_merge/ReindexAmbiguity.cpp index 55336cf77..73a26cbf1 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.cpp +++ b/image_analysis/scale_merge/ReindexAmbiguity.cpp @@ -41,12 +41,9 @@ namespace { } } -std::vector ReindexAmbiguityOperators(const UnitCell &cell, int space_group_number, +std::vector ReindexAmbiguityOperators(const UnitCell &cell, const gemmi::SpaceGroup &space_group, double max_obliquity_deg) { - const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); - if (sg == nullptr) - return {}; - return gemmi::find_twin_laws(static_cast(cell), sg, max_obliquity_deg, + return gemmi::find_twin_laws(static_cast(cell), &space_group, max_obliquity_deg, /*all_ops=*/false); } @@ -65,14 +62,11 @@ std::vector ReindexReflections(const std::vector ReindexMergedIntoAsu(const std::vector &merged, - const gemmi::Op &op, int space_group_number, + const gemmi::Op &op, const gemmi::SpaceGroup &space_group, bool merge_friedel) { - const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); - if (sg == nullptr) - return merged; // merge_friedel=false so the generator reports which side of the Friedel pair the ASU index was // reached from; that sign is what the two hands have to follow. - const HKLKeyGenerator key_gen(/*merge_friedel=*/false, *sg); + const HKLKeyGenerator key_gen(/*merge_friedel=*/false, space_group); std::vector out = merged; for (auto &r : out) { @@ -102,10 +96,10 @@ std::vector ReindexMergedIntoAsu(const std::vector &merged, - const UnitCell &cell, int space_group_number, + const UnitCell &cell, const gemmi::SpaceGroup &space_group, const std::function &)> &score, double max_obliquity_deg) { - const auto ops = ReindexAmbiguityOperators(cell, space_group_number, max_obliquity_deg); + const auto ops = ReindexAmbiguityOperators(cell, space_group, max_obliquity_deg); const BestReindex best = PickBestReindex( score(merged), ops, [&](const gemmi::Op &op) { return score(ReindexReflections(merged, op)); }); @@ -121,11 +115,8 @@ ReindexChoice ChooseReindex(const std::vector &merged, double ReferenceIntensityCC(const std::vector &merged, const std::vector &reference, - int space_group_number) { - const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); - if (sg == nullptr) - return 0.0; - const HKLKeyGenerator key(/*merge_friedel=*/true, *sg); + const gemmi::SpaceGroup &space_group) { + const HKLKeyGenerator key(/*merge_friedel=*/true, space_group); std::unordered_map ref; ref.reserve(reference.size()); @@ -156,11 +147,11 @@ double ReferenceIntensityCC(const std::vector &merged, ReindexAmbiguityResolver::ReindexAmbiguityResolver(const DiffractionExperiment &x, const std::vector &reference) : s(x.GetScalingSettings()), - hkl_key_generator(s.GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)) { + hkl_key_generator(s.GetMergeFriedel(), x.GetSpaceGroupOrP1()) { for (const auto &r : reference) reference_data[hkl_key_generator(r)] = r.I; - if (x.GetUnitCell().has_value() && x.GetSpaceGroupNumber().has_value()) - ops = ReindexAmbiguityOperators(*x.GetUnitCell(), static_cast(*x.GetSpaceGroupNumber())); + if (x.GetUnitCell().has_value() && x.GetGemmiSpaceGroup().has_value()) + ops = ReindexAmbiguityOperators(*x.GetUnitCell(), x.GetSpaceGroupOrP1()); } bool ReindexAmbiguityResolver::Accept(const Reflection &r) const { diff --git a/image_analysis/scale_merge/ReindexAmbiguity.h b/image_analysis/scale_merge/ReindexAmbiguity.h index 354e55982..1fb274387 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.h +++ b/image_analysis/scale_merge/ReindexAmbiguity.h @@ -25,7 +25,7 @@ // The candidate reindexing operators = twin laws = cosets of the crystal point group in the lattice // symmetry. Empty when the lattice symmetry equals the crystal Laue symmetry (holohedral crystal), // i.e. there is no indexing ambiguity. Identity is NOT included in the returned list. -std::vector ReindexAmbiguityOperators(const UnitCell &cell, int space_group_number, +std::vector ReindexAmbiguityOperators(const UnitCell &cell, const gemmi::SpaceGroup &space_group, double max_obliquity_deg = 2.0); // A copy of `merged` with each Miller index transformed by `op` (the intensities are unchanged). @@ -34,12 +34,12 @@ std::vector ReindexAmbiguityOperators(const UnitCell &cell, int space std::vector ReindexReflections(const std::vector &merged, const gemmi::Op &op); -// As above, but the reindexed reflections are brought back into the ASU of `space_group_number` - +// As above, but the reindexed reflections are brought back into the ASU of `space_group` - // what a reflection file has to be written in. Where the reduction reaches a reflection through its // Friedel mate, the Bijvoet halves travel with it. `merge_friedel` is the merge's own setting, since // it fixes where the minus hand is stored (at -hkl when the mates are kept apart). std::vector ReindexMergedIntoAsu(const std::vector &merged, - const gemmi::Op &op, int space_group_number, + const gemmi::Op &op, const gemmi::SpaceGroup &space_group, bool merge_friedel); struct ReindexChoice { @@ -54,15 +54,15 @@ struct ReindexChoice { // `score` is higher-is-better (e.g. correlation against a reference, or minus the R-free against a // model). Returns identity when there is no ambiguity or identity scores best. ReindexChoice ChooseReindex(const std::vector &merged, - const UnitCell &cell, int space_group_number, + const UnitCell &cell, const gemmi::SpaceGroup &space_group, const std::function &)> &score, double max_obliquity_deg = 2.0); // Ready-made scorer: Pearson correlation of merged intensities against a reference set, matched by the -// Friedel-merged (Laue) ASU index in `space_group_number`. Returns 0 when too few reflections match. +// Friedel-merged (Laue) ASU index in `space_group`. Returns 0 when too few reflections match. double ReferenceIntensityCC(const std::vector &merged, const std::vector &reference, - int space_group_number); + const gemmi::SpaceGroup &space_group); // Per-image resolver for serial stills. Each crystal is indexed independently in one of the // merohedrally-equivalent hands at random, so the ambiguity has to be broken per image (unlike rotation diff --git a/image_analysis/scale_merge/RfreeFlags.cpp b/image_analysis/scale_merge/RfreeFlags.cpp index fc5b7c41d..6b7b17cd2 100644 --- a/image_analysis/scale_merge/RfreeFlags.cpp +++ b/image_analysis/scale_merge/RfreeFlags.cpp @@ -22,7 +22,7 @@ namespace { } } -void AssignRfreeFlags(std::vector &merged, int32_t space_group_number, +void AssignRfreeFlags(std::vector &merged, const gemmi::SpaceGroup &space_group, double rfree_fraction, int min_free_reflections) { for (auto &r : merged) r.rfree_flag = false; @@ -36,7 +36,7 @@ void AssignRfreeFlags(std::vector &merged, int32_t space_group // dataset of one crystal form gets the SAME free set, which is what a multi-dataset campaign // (ensemble refinement, PanDDA) needs. A uniform hash draws ~rfree_fraction of the distinct // reflections free; a stratified per-shell draw would be tied to the dataset and break that. - const HKLKeyGenerator laue_key(/*merge_friedel=*/true, space_group_number); + const HKLKeyGenerator laue_key(/*merge_friedel=*/true, space_group); // Count the distinct test-eligible reflections (distinct Laue-ASU keys; mates collapse to one) so // the fraction can be floored to a usable test-set size on small data. @@ -57,11 +57,11 @@ void AssignRfreeFlags(std::vector &merged, int32_t space_group r.rfree_flag = UniformFromKey(laue_key(r).pack()) < eff_fraction; } -size_t ApplyReferenceFreeFlags(std::vector &merged, int32_t space_group_number, +size_t ApplyReferenceFreeFlags(std::vector &merged, const gemmi::SpaceGroup &space_group, const std::vector &reference) { // Reference free/work partition keyed by the Friedel-merged (Laue) ASU index, so it transfers // regardless of which Bijvoet mate / symmetry equivalent each dataset happens to have measured. - const HKLKeyGenerator laue_key(/*merge_friedel=*/true, space_group_number); + const HKLKeyGenerator laue_key(/*merge_friedel=*/true, space_group); std::unordered_map ref_flag; ref_flag.reserve(reference.size()); for (const auto &r : reference) diff --git a/image_analysis/scale_merge/RfreeFlags.h b/image_analysis/scale_merge/RfreeFlags.h index 8b10e1e18..fd0402577 100644 --- a/image_analysis/scale_merge/RfreeFlags.h +++ b/image_analysis/scale_merge/RfreeFlags.h @@ -8,6 +8,7 @@ #include #include "../../common/Reflection.h" +#include "gemmi/symmetry.hpp" // Assign R-free (test-set) flags to merged reflections. The flag is a pure function of the // Friedel-merged (Laue) ASU key, which gives three properties: @@ -26,7 +27,7 @@ // the floor only lifts the fraction on genuinely small datasets, where per-dataset R-free stability // matters more than cross-dataset identity (and a shared reference/ApplyReferenceFreeFlags is the way // to keep exact identity there). Pass 0 to disable the floor. -void AssignRfreeFlags(std::vector &merged, int32_t space_group_number, +void AssignRfreeFlags(std::vector &merged, const gemmi::SpaceGroup &space_group, double rfree_fraction, int min_free_reflections = 500); // Override the R-free flags on `merged` from a reference set (e.g. a reference MTZ's FreeR_flag @@ -34,5 +35,5 @@ void AssignRfreeFlags(std::vector &merged, int32_t space_group // (Laue) ASU key matches a reference reflection takes the reference's flag; reflections with no // match keep the flag they already carry (the deterministic per-hkl hash from AssignRfreeFlags). // Returns the number of merged reflections matched to the reference (for logging). -size_t ApplyReferenceFreeFlags(std::vector &merged, int32_t space_group_number, +size_t ApplyReferenceFreeFlags(std::vector &merged, const gemmi::SpaceGroup &space_group, const std::vector &reference); diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 5ad43d00d..003842e8e 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -2631,14 +2631,12 @@ namespace { // a Miller vector of the survivors first is that many more megabytes for a column of ten integers. // Inlined, the h planes are independent, so one thread takes each and they are summed back in h // order - integer addition, so the counts are the serial ones whatever the split. - void PossiblePerShell(int space_group_number, const UnitCell &cell, double d_min, double d_max, + void PossiblePerShell(const gemmi::SpaceGroup &sg, const UnitCell &cell, double d_min, double d_max, const ResolutionShells &shells, bool merge_friedel, std::vector &possible, size_t nthreads) { const gemmi::UnitCell gemmi_cell = cell; - const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); - if (sg == nullptr) return; - const gemmi::GroupOps gops = sg->operations(); - const gemmi::ReciprocalAsu asu(sg); + const gemmi::GroupOps gops = sg.operations(); + const gemmi::ReciprocalAsu asu(&sg); const gemmi::Miller lim = gemmi_cell.get_hkl_limits(d_min); const double inv_dmin2 = 1.0 / (d_min * d_min); const double inv_dmax2 = d_max > 0 ? 1.0 / (d_max * d_max) : 0.0; @@ -3210,7 +3208,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // is not written needs neither. French-Wilson is deferred until after the anomalous accumulator // below has attached I(+)/I(-), so the two hands get their amplitudes in one pass. if (full_stats) - AssignRfreeFlags(result.merged, x.GetSpaceGroupNumber().value_or(1), rfree_fraction); + AssignRfreeFlags(result.merged, x.GetSpaceGroupOrP1(), rfree_fraction); if (reject_count > 0) logger.Info("Merge outlier rejection: dropped {} observations", reject_count); @@ -3243,7 +3241,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // than an asymmetric unit of it - so enumerating it there is the largest single piece of work in // the merge that nothing goes on to read. if (reference_cell && !for_search) - PossiblePerShell(x.GetSpaceGroupNumber().value_or(1), *reference_cell, d_min_pad, d_max_pad, + PossiblePerShell(x.GetSpaceGroupOrP1(), *reference_cell, d_min_pad, d_max_pad, shells, merge_friedel, possible, nthreads); for (int s = 0; s < n_shells; ++s) sa[s].possible = possible[s]; @@ -3341,9 +3339,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool double sig_num_all = 0.0, sig_den_all = 0.0; size_t sig_n = 0; if (!for_search && full_stats) { - const int sg_num = x.GetSpaceGroupNumber().value_or(1); - const HKLKeyGenerator anom_keygen(/*merge_friedel=*/false, sg_num); - const gemmi::GroupOps gops = gemmi::find_spacegroup_by_number(sg_num)->operations(); + const HKLKeyGenerator anom_keygen(/*merge_friedel=*/false, x.GetSpaceGroupOrP1()); + const gemmi::GroupOps gops = x.GetSpaceGroupOrP1().operations(); struct AnomAcc { double swI[2] = {}; double sw[2] = {}; int32_t h = 0, k = 0, l = 0; float d = NAN; }; // [hand] 0=I(+) 1=I(-) // The Friedel-merged key an observation lands on, and which hand it is, are properties of its // ASU GROUP, not of the observation - group_h/k/l is that group's SIGNED representative, so the @@ -3444,7 +3441,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // Attach the per-reflection anomalous split so the writer can emit I(+)/I(-) by default (each merged // reflection maps to its Friedel-ASU key; in an anomalous merge both mates map to the same key). if (!anom_export.empty()) { - const HKLKeyGenerator anom_keygen(/*merge_friedel=*/false, x.GetSpaceGroupNumber().value_or(1)); + const HKLKeyGenerator anom_keygen(/*merge_friedel=*/false, x.GetSpaceGroupOrP1()); for (auto &r : result.merged) { const HKLKey ak = anom_keygen(r.h, r.k, r.l); const auto it = anom_export.find(HKLKey{ak.h, ak.k, ak.l, true}.pack()); @@ -3458,7 +3455,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool if (full_stats) { FrenchWilsonOptions fw_opts; fw_opts.num_threads = static_cast(nthreads); - ApplyFrenchWilson(result.merged, x.GetSpaceGroupNumber().value_or(1), fw_opts); + ApplyFrenchWilson(result.merged, x.GetSpaceGroupOrP1(), fw_opts); } logger.Info("Merge complete ({} unique reflections)", result.merged.size()); @@ -3467,8 +3464,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_stats, bool measure_cc_before_corrections) { - const int sg_number = x.GetSpaceGroupNumber().value_or(1); - HKLKeyGenerator keygen(merge_friedel, sg_number); + HKLKeyGenerator keygen(merge_friedel, x.GetSpaceGroupOrP1()); // Start from the corr Ingest built, so this pass runs the scaling_iter iterations it was asked for // rather than continuing the previous pass's (see the header). diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index b1dc18872..6801b818a 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -77,7 +77,7 @@ namespace { ScaleOnTheFly::ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref) : s(x.GetScalingSettings()), - hkl_key_generator(s.GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)) { + hkl_key_generator(s.GetMergeFriedel(), x.GetSpaceGroupOrP1()) { for (const auto &r: ref) { const auto key = hkl_key_generator(r); reference_data[key] = r.I; diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp index 6396d35eb..a7166317a 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.cpp +++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp @@ -141,7 +141,7 @@ namespace { StillsPartialityRefine::StillsPartialityRefine(const DiffractionExperiment &x) : experiment_(x), - hkl_key_generator_(x.GetScalingSettings().GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)), + hkl_key_generator_(x.GetScalingSettings().GetMergeFriedel(), x.GetSpaceGroupOrP1()), d_min_limit_(x.GetScalingSettings().GetHighResolutionLimit_A()), d_max_limit_(x.GetScalingSettings().GetLowResolutionLimit_A()), min_partiality_(x.GetScalingSettings().GetMinPartiality()), diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 65a0bbeb6..a7a531f12 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -902,7 +902,12 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen .beta = tmp[4], .gamma = tmp[5] }); - dataset->experiment.SpaceGroupNumber(master_file->GetOptInt("/entry/sample/space_group_number")); + // The name carries the setting, the number cannot ("R 3:R" reads back as "R 3:H"), so the + // name is preferred; the number is the fallback for a file written before it was recorded. + if (const auto *sg = gemmi::find_spacegroup_by_name(master_file->GetString("/entry/sample/space_group"))) + dataset->experiment.SetSpaceGroup(*sg); + else + dataset->experiment.SpaceGroupNumber(master_file->GetOptInt("/entry/sample/space_group_number")); // The setting the cell and space group just read are in, relative to the setting the per-image // reflections and lattices were written in. Absent on every file written before the offline // analysis started recording it, and on every run that never re-seated its lattice - both mean diff --git a/reader/JFJochReader.cpp b/reader/JFJochReader.cpp index fb49e3588..11da4bd96 100644 --- a/reader/JFJochReader.cpp +++ b/reader/JFJochReader.cpp @@ -83,7 +83,7 @@ void JFJochReader::UpdateGeomMetadata(const DiffractionExperiment &experiment) { new_dataset->experiment.PoniRot2_rad(experiment.GetDatasetSettings().GetPoniRot2_rad()); new_dataset->experiment.PoniRot3_rad(experiment.GetDatasetSettings().GetPoniRot3_rad()); new_dataset->experiment.SetUnitCell(experiment.GetUnitCell()); - new_dataset->experiment.SpaceGroupNumber(experiment.GetSpaceGroupNumber()); + new_dataset->experiment.SetSpaceGroup(experiment.GetGemmiSpaceGroup()); new_dataset->experiment.PolarizationFactor(experiment.GetPolarizationFactor()); new_dataset->experiment.Goniometer(experiment.GetGoniometer()); new_dataset->experiment.GridScan(experiment.GetGridScan()); diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 6c2e02526..e00576174 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -70,7 +70,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & const std::string &model_path, const std::string &output_prefix, Logger &logger, - std::optional data_space_group_number, + const gemmi::SpaceGroup *data_space_group, bool probe_indexing_ambiguity) { ModelValidationResult result; result.model_path = model_path; @@ -116,18 +116,17 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // disagreement but manufactures agreement. The anomalous difference map below is the only honest // arbiter, and it is used to report the disagreement rather than to bury it. const std::vector &obs = merged; - if (data_space_group_number && *data_space_group_number != sg->number) { - const gemmi::SpaceGroup *dsg = gemmi::find_spacegroup_by_number(*data_space_group_number); - if (dsg && dsg->is_enantiomorphic() && sg->is_enantiomorphic()) { - gemmi::GroupOps eops = dsg->operations(); - eops.change_basis_forward(dsg->change_of_hand_op()); + if (data_space_group && data_space_group->number != sg->number) { + if (data_space_group->is_enantiomorphic() && sg->is_enantiomorphic()) { + gemmi::GroupOps eops = data_space_group->operations(); + eops.change_basis_forward(data_space_group->change_of_hand_op()); const gemmi::SpaceGroup *enant = gemmi::find_spacegroup_by_ops(eops); if (enant && enant->number == sg->number) { result.adopted_model_enantiomorph = true; logger.Info("Model validation: data space group {} is the enantiomorph of the model {}; " "taking the model's group as the label, with no reindexing - the two groups " "index identically, and reindexing would flip the anomalous differences", - dsg->short_name(), sg->hm); + data_space_group->short_name(), sg->hm); } } } @@ -263,7 +262,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // it cannot distinguish them - that is taken from the model hand above. Fit best = fit_model(obs); if (probe_indexing_ambiguity) { - const auto reindex_ops = ReindexAmbiguityOperators(cell, sg->number); + const auto reindex_ops = ReindexAmbiguityOperators(cell, *sg); std::vector candidate_r_free{best.r_free}; // identity first, then the twin laws for (const auto &op : reindex_ops) { Fit cand = fit_model(ReindexReflections(obs, op)); @@ -320,10 +319,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // Friedel ASU of that frame, so an anomalous merge - which keeps each mate as a row of its own, // both carrying the same pair - is read on its + rows only. Taking the - rows as well would // give one reflection both signs of its difference, and the last row written would decide. - const gemmi::SpaceGroup *data_sg = data_space_group_number - ? gemmi::find_spacegroup_by_number(*data_space_group_number) : nullptr; - if (data_sg == nullptr) - data_sg = sg; + const gemmi::SpaceGroup *data_sg = data_space_group != nullptr ? data_space_group : sg; const gemmi::ReciprocalAsu data_asu(data_sg); const gemmi::GroupOps data_gops = data_sg->operations(); @@ -471,14 +467,14 @@ std::string hkl_triplet(const gemmi::Op &op) { } // namespace -int AdoptModelFrame(const ModelValidationResult &validation, - std::vector &merged, - int data_space_group_number, - bool merge_friedel, - Logger &logger) { - int space_group_number = data_space_group_number; +const gemmi::SpaceGroup *AdoptModelFrame(const ModelValidationResult &validation, + std::vector &merged, + const gemmi::SpaceGroup &data_space_group, + bool merge_friedel, + Logger &logger) { + const gemmi::SpaceGroup *space_group = &data_space_group; if (!validation.ok) - return space_group_number; + return space_group; // Adopting the model's enantiomorph is a change of the space-group LABEL and nothing else. The // two groups have the same rotation operations, so the same reflections, indexed the way they @@ -486,24 +482,24 @@ int AdoptModelFrame(const ModelValidationResult &validation, // group that agrees with the model it will be refined against. Reindexing here would swap the // Bijvoet mates and so change the data - see the note in ValidateAgainstModel. if (validation.adopted_model_enantiomorph && validation.model_space_group_number > 0) { - space_group_number = validation.model_space_group_number; - const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); + space_group = gemmi::find_spacegroup_by_number(validation.model_space_group_number); logger.Info("Model validation: the written reflections take the model's enantiomorph, {} ({}), " - "as a label - no reflection moved", sg ? sg->short_name() : "?", space_group_number); + "as a label - no reflection moved", + space_group ? space_group->short_name() : "?", validation.model_space_group_number); } // The alternative indexing, by contrast, is metric- and group-preserving: only the labels move. if (!(validation.indexing_op == gemmi::Op::identity())) { - merged = ReindexMergedIntoAsu(merged, validation.indexing_op, space_group_number, merge_friedel); + merged = ReindexMergedIntoAsu(merged, validation.indexing_op, *space_group, merge_friedel); logger.Info("Model validation: the written reflections take the model's indexing, reindexed by {}", hkl_triplet(validation.indexing_op)); } - return space_group_number; + return space_group; } std::vector ModelReferenceIntensities(const std::string &model_path, const std::optional &cell, - std::optional space_group_number, + const gemmi::SpaceGroup *space_group, double d_min, Logger &logger) { std::vector out; @@ -542,9 +538,8 @@ std::vector ModelReferenceIntensities(const std::string &model st.cell = target; } } - if (space_group_number.has_value()) - if (const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(*space_group_number)) - st.spacegroup_hm = sg->xhm(); + if (space_group != nullptr) + st.spacegroup_hm = space_group->xhm(); const gemmi::SpaceGroup *sg = st.find_spacegroup(); if (!sg) { logger.Error("Model reference: model {} has no usable space group", model_path); diff --git a/rugnux/ModelValidation.h b/rugnux/ModelValidation.h index ce1e38104..24f5073ec 100644 --- a/rugnux/ModelValidation.h +++ b/rugnux/ModelValidation.h @@ -68,7 +68,7 @@ struct ModelValidationResult { // No refinement of the structure itself is done. The model is only re-fractionalized into the // data unit cell (a cheap rigid cell adjustment) so a deposited model with a slightly different // cell still lines up with the processed data. Returns ok=false (and logs) on any failure. -// data_space_group_number is the space group the data was merged in. If it is the enantiomorph of +// data_space_group is the space group the data was merged in. If it is the enantiomorph of // the model's space group (e.g. data P4(1)2(1)2, model P4(3)2(1)2 - indistinguishable from merged // intensities), the model's group is adopted as a label, without touching the reflections. // @@ -82,27 +82,27 @@ ModelValidationResult ValidateAgainstModel(const std::vector & const std::string &model_path, const std::string &output_prefix, Logger &logger, - std::optional data_space_group_number = std::nullopt, + const gemmi::SpaceGroup *data_space_group = nullptr, bool probe_indexing_ambiguity = true); // Reindex `merged` into the frame ValidateAgainstModel reported, so the reflection files that are // written describe the same indexing as the R-factors and the maps. Returns the space group they are -// now in - the model's, where its enantiomorph was adopted, and `data_space_group_number` otherwise +// now in - the model's, where its enantiomorph was adopted, and `data_space_group` otherwise // (a change of indexing is metric- and group-preserving). A no-op where the model needed neither. -int AdoptModelFrame(const ModelValidationResult &validation, - std::vector &merged, - int data_space_group_number, - bool merge_friedel, - Logger &logger); +const gemmi::SpaceGroup *AdoptModelFrame(const ModelValidationResult &validation, + std::vector &merged, + const gemmi::SpaceGroup &data_space_group, + bool merge_friedel, + Logger &logger); // Reference intensities computed from an atomic model: |F_model|^2, with a flat bulk-solvent // contribution at the standard constants, for every reflection to `d_min`. This is the same // information a reference MTZ carries, and it is what lets a model break the indexing ambiguity of a // serial-stills run: there the ambiguity has to be broken per image, at integration time, long before -// there is a merge to fit a model to. `cell` and `space_group_number` override the model's own, where +// there is a merge to fit a model to. `cell` and `space_group` override the model's own, where // the run already knows them. Empty on any failure (which is logged). std::vector ModelReferenceIntensities(const std::string &model_path, const std::optional &cell, - std::optional space_group_number, + const gemmi::SpaceGroup *space_group, double d_min, Logger &logger); diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 07822fe5a..cb2705db3 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -139,8 +139,8 @@ 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_number.has_value()) - Key(os, "SPACE_GROUP_NUMBER", result.space_group_number.value()); + if (result.space_group.has_value()) + Key(os, "SPACE_GROUP_NUMBER", result.space_group->number); if (result.indexing_rate.value_or(0.0f) <= 0.0f) warnings.emplace_back("No image indexed - no crystal lattice was determined from this dataset"); @@ -196,7 +196,7 @@ std::string RenderResultReport(const std::string &output_prefix, search.best_space_group->centring_type(), search.best_space_group->short_name(), search.best_space_group->centring_type())); - } else if (result.space_group_number.has_value()) { + } else if (result.space_group.has_value()) { Key(os, "SPACE_GROUP_SEARCH", "FIXED"); os << "\n The space group was given, not determined here.\n"; } else { diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 7c09187da..b516d7e7d 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -398,7 +398,7 @@ Rugnux::Rugnux(JFJochReader &reader, DiffractionExperiment experiment, PixelMask pixel_mask, ProcessConfig config) : reader_(reader), experiment_(std::move(experiment)), pixel_mask_(std::move(pixel_mask)), config_(std::move(config)), - user_fixed_sg_(experiment_.GetSpaceGroupNumber()) { + user_fixed_sg_(experiment_.GetGemmiSpaceGroup()) { // Bit 9 describes where THIS run found the beam stop, so a mask read back from a file that // already carries one starts clear; the user mask (bit 8) is left as it was loaded. pixel_mask_.ClearBeamStopMask(experiment_); @@ -1526,7 +1526,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // ASKED for. A pass that determined one leaves it on experiment_, and the pass after it - the // re-runs in Run(), or a second Run() on the same object - would otherwise inherit that answer and // skip the search that is the whole point of running again. - experiment_.SpaceGroupNumber(user_fixed_sg_); + experiment_.SetSpaceGroup(user_fixed_sg_); const auto dataset = reader_.GetDataset(); if (!dataset) @@ -1735,12 +1735,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (!config_.reference_data.empty()) { indexer->ReferenceIntensities(config_.reference_data); } else if (!config_.model_path.empty() && !experiment_.IsRotationIndexing() - && experiment_.GetUnitCell().has_value() && experiment_.GetSpaceGroupNumber().has_value()) { + && experiment_.GetUnitCell().has_value() && experiment_.GetGemmiSpaceGroup().has_value()) { const auto d_min = experiment_.GetScalingSettings().GetHighResolutionLimit_A() .value_or(experiment_.GetDetectorMaxResolution_A()); model_reference = ModelReferenceIntensities( config_.model_path, experiment_.GetUnitCell(), - static_cast(*experiment_.GetSpaceGroupNumber()), d_min, logger); + &experiment_.GetSpaceGroupOrP1(), d_min, logger); if (!model_reference.empty()) { indexer->ReferenceIntensities(model_reference); ambiguity_reference = true; @@ -3059,13 +3059,13 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // A group the RUN determined is left alone - the arm above and the centring check below are its // handling, and a de-novo run must behave exactly as it did. if (user_fixed_sg_ && end_msg.rotation_lattice.has_value() && end_msg.rotation_lattice_type.has_value()) { - const auto *fixed_sg = gemmi::find_spacegroup_by_number(static_cast(*user_fixed_sg_)); + const gemmi::SpaceGroup &fixed_sg = *user_fixed_sg_; const char indexed_centering = end_msg.rotation_lattice_type->centering; - if (fixed_sg && indexed_centering != fixed_sg->centring_type()) { + if (indexed_centering != fixed_sg.centring_type()) { // A trigonal-P group sits on a hexagonal-P lattice; every other system names its own. - const auto want_system = fixed_sg->crystal_system() == gemmi::CrystalSystem::Trigonal - && fixed_sg->centring_type() == 'P' - ? gemmi::CrystalSystem::Hexagonal : fixed_sg->crystal_system(); + const auto want_system = fixed_sg.crystal_system() == gemmi::CrystalSystem::Trigonal + && fixed_sg.centring_type() == 'P' + ? gemmi::CrystalSystem::Hexagonal : fixed_sg.crystal_system(); // Search the PRIMITIVE cell. rotation_lattice is the conventional one and its centring // is carried beside it, but the character table is reached through a Niggli reduction // that takes whatever it is handed to be primitive - so a centred conventional cell @@ -3077,13 +3077,13 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // The other lattice arithmetic in this file goes through ToPrimitive for the same reason. const auto cand = LatticeSearchForClass( end_msg.rotation_lattice->ToPrimitive(indexed_centering), want_system, - fixed_sg->centring_type()); + fixed_sg.centring_type()); if (cand && reindex_into(*cand)) { const auto &uc = *result.consensus_cell; logger.Info("Reindexed the {}-centred indexed lattice into the {}-centred setting the " "fixed space group {} needs: a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} " "beta={:.2f} gamma={:.2f}", indexed_centering, cand->centering, - fixed_sg->xhm(), uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma); + fixed_sg.xhm(), uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma); } } } @@ -3102,15 +3102,15 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // the group is the user's assertion, and quietly determining a different one would answer a // question that was not asked. Name the cell that WAS indexed so the user can act on it. if (user_fixed_sg_ && end_msg.rotation_lattice_type.has_value() && end_msg.unit_cell.has_value()) { - const auto *fixed_sg = gemmi::find_spacegroup_by_number(static_cast(*user_fixed_sg_)); - if (fixed_sg && end_msg.rotation_lattice_type->centering != fixed_sg->centring_type()) { + const gemmi::SpaceGroup &fixed_sg = *user_fixed_sg_; + if (end_msg.rotation_lattice_type->centering != fixed_sg.centring_type()) { const auto &uc = *end_msg.unit_cell; if (prepass_result_) { logger.Warning("Two-pass: the second pass indexed a {}-centred lattice, but the fixed " "space group {} has a {}-centred lattice - the group cannot describe " "this lattice, so this pass will not be adopted", - end_msg.rotation_lattice_type->centering, fixed_sg->xhm(), - fixed_sg->centring_type()); + end_msg.rotation_lattice_type->centering, fixed_sg.xhm(), + fixed_sg.centring_type()); result.lattice_conflicts_with_fixed_sg = true; } else { throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, fmt::format( @@ -3121,9 +3121,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b "its frame, so nothing it reported would describe this crystal. Re-run without a fixed " "space group to have it determined from the data, or fix one whose lattice this " "crystal has.", - fixed_sg->xhm(), fixed_sg->centring_type(), end_msg.rotation_lattice_type->centering, + fixed_sg.xhm(), fixed_sg.centring_type(), end_msg.rotation_lattice_type->centering, gemmi::crystal_system_str(end_msg.rotation_lattice_type->crystal_system), - uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, fixed_sg->xhm())); + uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, fixed_sg.xhm())); } } } @@ -3630,20 +3630,23 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b rsm->Ingest(); } } - experiment_.SpaceGroupNumber(sg.number); - end_msg.space_group_number = sg.number; - result.space_group_number = sg.number; + experiment_.SetSpaceGroup(sg); + end_msg.space_group_name = sg.xhm(); + end_msg.space_group_number = static_cast(sg.number); + result.space_group = sg; phase("Re-scaling in space group " + sg.short_name()); sm = scale_and_merge(sg.short_name(), false); } result.space_group_search = sg_search; } else { // A space group was fixed by the user; surface it so the viewer/CLI can still show it. - result.space_group_number = experiment_.GetSpaceGroupNumber(); + result.space_group = experiment_.GetGemmiSpaceGroup(); // ...and record it, exactly as the search arm above does with the group it determined, or // the _process.h5 comes out with a cell but no group and re-merging it defaults to P1. - if (const auto sg = experiment_.GetSpaceGroupNumber(); sg.has_value()) - end_msg.space_group_number = static_cast(*sg); + if (const auto &sg = experiment_.GetGemmiSpaceGroup(); sg.has_value()) { + end_msg.space_group_name = sg->xhm(); + end_msg.space_group_number = static_cast(sg->number); + } } // Reference-based indexing-ambiguity resolution (rotation). When a reference MTZ is supplied and @@ -3653,11 +3656,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // Stills resolve the ambiguity per image with ReindexAmbiguityResolver, not here. if (rsm && !config_.reference_data.empty() && result.consensus_cell && experiment_.GetSpaceGroupNumber().has_value()) { - const int sg_num = static_cast(*experiment_.GetSpaceGroupNumber()); + const gemmi::SpaceGroup &sg = experiment_.GetSpaceGroupOrP1(); const auto choice = ChooseReindex( - sm.merged, *result.consensus_cell, sg_num, + sm.merged, *result.consensus_cell, sg, [&](const std::vector &m) { - return ReferenceIntensityCC(m, config_.reference_data, sg_num); + return ReferenceIntensityCC(m, config_.reference_data, sg); }); if (!choice.is_identity) { logger.Info("Reference: resolved indexing ambiguity by reindexing to match the reference " @@ -3676,9 +3679,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } } - const auto twin_sg_number = experiment_.GetSpaceGroupNumber(); - const gemmi::SpaceGroup *twin_sg = twin_sg_number - ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr; + const auto &twin_sg_opt = experiment_.GetGemmiSpaceGroup(); + const gemmi::SpaceGroup *twin_sg = twin_sg_opt ? &*twin_sg_opt : nullptr; // Not on the geometry pre-pass: the analysis goes into that pass's statistics text and its // written reflections, and neither survives the run. The promotion flag below is a different // thing - it is what the SEARCH did, the second pass reads it, and it is set either way. @@ -3729,9 +3731,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // operators. For an obvious merohedral case (P3/P4/P6...) users expect this; but a PSEUDO-merohedral // metric (e.g. a C2 crystal whose beta makes it pseudo-F-orthorhombic) is easy to miss, so surface it. The // reindex operator is applied directly to (h,k,l), so its triplet reads as an h,k,l transform. - if (result.consensus_cell && twin_sg_number) { - const auto twin_ops = ReindexAmbiguityOperators(*result.consensus_cell, - static_cast(*twin_sg_number), 2.0); + if (result.consensus_cell && twin_sg != nullptr) { + const auto twin_ops = ReindexAmbiguityOperators(*result.consensus_cell, *twin_sg, 2.0); if (!twin_ops.empty()) { std::string laws; for (const auto &op : twin_ops) { @@ -3840,8 +3841,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // Inherit the campaign's shared R-free test set from the reference MTZ (overriding the // per-hkl hash that the merge assigned), so every dataset flags the same free reflections. if (config_.reference_has_free_flags && !config_.reference_data.empty() && !sm.merged.empty()) { - const auto sg = experiment_.GetSpaceGroupNumber().value_or(1); - const size_t matched = ApplyReferenceFreeFlags(sm.merged, static_cast(sg), + const size_t matched = ApplyReferenceFreeFlags(sm.merged, experiment_.GetSpaceGroupOrP1(), config_.reference_data); logger.Info("R-free flags: inherited the reference test set ({} of {} merged reflections matched)", matched, sm.merged.size()); @@ -3945,13 +3945,13 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // the R-factors and the maps then all describe one indexing instead of two. if (result.consensus_cell && write_files && !config_.model_path.empty() && !geometry_prepass) { phase("Validating against model"); - const auto data_sg = experiment_.GetSpaceGroupNumber(); + const auto data_sg = experiment_.GetGemmiSpaceGroup(); // With a reference MTZ the merohedral indexing was already resolved against it (rotation // merge / stills scaling), so trust that; only probe indexing by R-free when model-only. const auto validation = ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path, config_.output_prefix, logger, - data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt, + data_sg ? &*data_sg : nullptr, /*probe_indexing_ambiguity=*/config_.reference_data.empty()); // A model that was asked for and could not be used has to say so where anyone will see // it. Without this the run ends successfully with no R-free, no maps and nothing in the @@ -3960,16 +3960,17 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (!validation.failure_reason.empty()) result.warnings.push_back("Model validation did not run: " + validation.failure_reason); if (data_sg.has_value()) { - const int adopted = AdoptModelFrame(validation, sm.merged, static_cast(*data_sg), - experiment_.GetScalingSettings().GetMergeFriedel(), logger); - if (adopted != static_cast(*data_sg)) { + const gemmi::SpaceGroup *adopted = + AdoptModelFrame(validation, sm.merged, *data_sg, + experiment_.GetScalingSettings().GetMergeFriedel(), logger); + if (adopted->number != data_sg->number) { // The written reflections and the report that describes them, but deliberately not // end_msg: the _process.h5 already holds the per-image reflections as they were // integrated, and a file whose recorded group did not match its own reflections // would mis-merge on a later --mode scale. An enantiomorphic pair merges identically // anyway - same Laue class, same absences - so nothing is lost by leaving it. - experiment_.SpaceGroupNumber(adopted); - result.space_group_number = adopted; + experiment_.SetSpaceGroup(*adopted); + result.space_group = *adopted; } // The unmerged export below comes from the integrated observations, not from the merge, // so it follows the same relabeling - otherwise the two files would describe different diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 064b7f21c..3733b5d57 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -251,8 +251,9 @@ struct ProcessResult { std::optional model_validation; // Space group used for the final re-scale/merge and written to the master file: determined by the - // search when the user did not fix one, otherwise the fixed group. - std::optional space_group_number; + // search when the user did not fix one, otherwise the fixed group. The whole group, so it keeps + // its setting - a number only ever names the reference one. + std::optional space_group; // Calibration mode: the geometry fitted to the powder rings of the whole run. std::optional calibration; @@ -340,7 +341,7 @@ class Rugnux { // determined. A group the user asserted is treated differently from one the run determined: it is // re-seated onto its own Bravais lattice if the metric has that setting, and the run refuses rather // than merging in it if the metric does not. Each pass is also restored to it before it starts. - const std::optional user_fixed_sg_; + const std::optional user_fixed_sg_; // Whether the beam centre has already been placed by a measurement the pre-scan cannot better: diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 12b18d939..9c3466732 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -699,7 +699,7 @@ static int RunRugnux(int argc, char **argv) { bool no_scaling_corrections = false; // --no-scaling-corrections: disable rot3d decay+absorption+modulation surfaces bool no_expected_variance_merge = false; // --no-expected-variance-merge: restore observed-sigma stills merge weighting bool anomalous_mode = false; - std::optional space_group_number; + const gemmi::SpaceGroup *space_group = nullptr; std::optional fixed_reference_unit_cell; std::optional max_spot_count_override; float sigma_spot_finding = 4.0; @@ -979,13 +979,14 @@ static int RunRugnux(int argc, char **argv) { break; case 'S': { // Accept a space-group number ("92") or a Hermann-Mauguin symbol ("P43212", "P 43 21 2"). + // The symbol is kept as the group it names, not reduced to its number: a number can only + // ever name the reference setting, so -S P1121 would silently become P 1 21 1. char *end = nullptr; const long as_number = strtol(optarg, &end, 10); - if (end != optarg && *end == '\0') { - space_group_number = as_number; - } else if (const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(optarg)) { - space_group_number = sg->number; - } else { + space_group = (end != optarg && *end == '\0') + ? gemmi::find_spacegroup_by_number(static_cast(as_number)) + : gemmi::find_spacegroup_by_name(optarg); + if (space_group == nullptr) { logger.Error("Unknown space group '{}' (use a number like 92 or a symbol like P43212)", optarg); exit(EXIT_FAILURE); } @@ -1339,16 +1340,8 @@ static int RunRugnux(int argc, char **argv) { nthreads = hw > 0 ? static_cast(hw) : 1; } - // Validate space group number early - const gemmi::SpaceGroup *space_group = nullptr; - if (space_group_number.has_value()) { - space_group = gemmi::find_spacegroup_by_number(space_group_number.value()); - if (!space_group) { - logger.Error("Unknown space group number {}", space_group_number.value()); - exit(EXIT_FAILURE); - } - logger.Info("Using space group {} (number {})", space_group->hm, space_group_number.value()); - } + if (space_group != nullptr) + logger.Info("Using space group {} (number {})", space_group->hm, space_group->number); // Off the critical path, so cuFFT's one-time initialisation overlaps the file open, the mask read // and the beam-stop pre-scan instead of landing inside the first pass. @@ -1424,18 +1417,17 @@ static int RunRugnux(int argc, char **argv) { mode == RugnuxMode::Scale ? dataset->experiment.GetUnitCell() : fixed_reference_unit_cell; const auto warning = ReferenceConsistencyWarning( reference, data_cell, - space_group_number.has_value() ? std::optional(static_cast(*space_group_number)) - : std::nullopt); + space_group != nullptr ? std::optional(space_group->number) : std::nullopt); if (!warning.empty()) logger.Warning("{}", warning); // A reference MTZ fixes the space group and unit cell, unless -S / -C override them. // (-S with the wrong enantiomorph, or -C with a different cell, is allowed - the explicit // flag always wins.) The cell is a soft reference: indexing may drift within tolerance. - if (!space_group_number.has_value() && reference.space_group_number.has_value()) { - space_group_number = static_cast(*reference.space_group_number); + if (space_group == nullptr && reference.space_group_number.has_value()) { + space_group = gemmi::find_spacegroup_by_number(static_cast(*reference.space_group_number)); logger.Info("Fixing space group from reference MTZ: {} ({})", - reference.space_group_name, *space_group_number); + reference.space_group_name, *reference.space_group_number); } if (!fixed_reference_unit_cell.has_value() && reference.cell.has_value()) { fixed_reference_unit_cell = reference.cell; @@ -1464,8 +1456,8 @@ static int RunRugnux(int argc, char **argv) { configure_offline_output(experiment, output_prefix); // The reflections in the file are already indexed, so the cell and space group they were // integrated in are the file's to supply here - but an explicit -S / -C still wins. - if (space_group_number.has_value()) - experiment.SpaceGroupNumber(space_group_number); + if (space_group != nullptr) + experiment.SetSpaceGroup(*space_group); if (fixed_reference_unit_cell.has_value()) experiment.SetUnitCell(fixed_reference_unit_cell); // A rotation (goniometer) dataset uses RotationScaleMerge unless --force-still asks for stills scaling. @@ -1708,8 +1700,7 @@ static int RunRugnux(int argc, char **argv) { // Inherit the campaign's shared R-free test set from the reference MTZ (overriding the // per-hkl hash the merge assigned), so every dataset flags the same free reflections. if (reference_has_free_flags && !reference_data.empty() && !merged_reflections.empty()) { - const auto sg = experiment.GetSpaceGroupNumber().value_or(1); - const size_t matched = ApplyReferenceFreeFlags(merged_reflections, static_cast(sg), + const size_t matched = ApplyReferenceFreeFlags(merged_reflections, experiment.GetSpaceGroupOrP1(), reference_data); logger.Info("R-free flags: inherited the reference test set ({} of {} merged reflections matched)", matched, merged_reflections.size()); @@ -1737,9 +1728,8 @@ static int RunRugnux(int argc, char **argv) { "Re-run rugnux (which determines and stores the space group) or pass " "-S to scale and merge in the correct symmetry."); - const auto twin_sg_number = experiment.GetSpaceGroupNumber(); - const gemmi::SpaceGroup *twin_sg = twin_sg_number - ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr; + const auto &twin_sg_opt = experiment.GetGemmiSpaceGroup(); + const gemmi::SpaceGroup *twin_sg = twin_sg_opt ? &*twin_sg_opt : nullptr; const auto twinning = AnalyzeTwinning(merged_reflections, twin_sg); std::cout << std::endl << TwinningAnalysisToText(twinning) << std::endl; @@ -1765,22 +1755,23 @@ static int RunRugnux(int argc, char **argv) { std::string model_validation_failure; std::optional model_validation; if (!output_prefix.empty() && !model_pdb.empty()) { - const auto data_sg = experiment.GetSpaceGroupNumber(); + const auto data_sg = experiment.GetGemmiSpaceGroup(); // With a reference MTZ the merohedral indexing was already resolved (stills per-image // scaling); only probe indexing by R-free when model-only, with no reference. const auto validation = ValidateAgainstModel(merged_reflections, *experiment.GetUnitCell(), model_pdb, output_prefix, logger, - data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt, + data_sg ? &*data_sg : nullptr, /*probe_indexing_ambiguity=*/reference_data.empty()); model_validation = validation; if (!validation.failure_reason.empty()) model_validation_failure = validation.failure_reason; if (data_sg.has_value()) { - const int adopted = AdoptModelFrame(validation, merged_reflections, static_cast(*data_sg), - experiment.GetScalingSettings().GetMergeFriedel(), logger); - if (adopted != static_cast(*data_sg)) - experiment.SpaceGroupNumber(adopted); + const gemmi::SpaceGroup *adopted = + AdoptModelFrame(validation, merged_reflections, *data_sg, + experiment.GetScalingSettings().GetMergeFriedel(), logger); + if (adopted->number != data_sg->number) + experiment.SetSpaceGroup(*adopted); // The unmerged export below writes the integrated observations, so it follows the same // relabeling as the merged set - only the indexing moves a reflection; the enantiomorph // is a label - or the two files would describe different indexings of one dataset. @@ -1825,7 +1816,7 @@ static int RunRugnux(int argc, char **argv) { ProcessResult scale_result; scale_result.images_processed = reader.GetNumberOfImages(); scale_result.consensus_cell = experiment.GetUnitCell(); - scale_result.space_group_number = experiment.GetSpaceGroupNumber(); + scale_result.space_group = experiment.GetGemmiSpaceGroup(); scale_result.used_beam_x_pxl = experiment.GetBeamX_pxl(); scale_result.used_beam_y_pxl = experiment.GetBeamY_pxl(); scale_result.used_distance_mm = experiment.GetDetectorDistance_mm(); @@ -2076,7 +2067,8 @@ static int RunRugnux(int argc, char **argv) { // factor, so an explicit --polarization has to land after it or it is silently overwritten. if (polarization_factor) experiment.PolarizationFactor(polarization_factor.value()); ClearStoredCrystal(experiment); // shared with the viewer; -S / -C below override it - experiment.SpaceGroupNumber(space_group_number); + if (space_group != nullptr) + experiment.SetSpaceGroup(*space_group); experiment.ImagesPerTrigger(images_to_process); // Re-determine the unit cell from scratch: discard any cell stored in the input file so @@ -2489,11 +2481,9 @@ static int RunRugnux(int argc, char **argv) { // search happened to return and no reflection was measured on it, so printing it as the run's answer // states a result the data do not support. const bool anything_indexed = result.indexing_rate.value_or(0.0f) > 0.0f; - if (result.space_group_number.has_value() && anything_indexed) { - const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number( - static_cast(result.space_group_number.value())); - std::string line = fmt::format("Space group: {} (No. {})", sg ? sg->short_name() : "?", - result.space_group_number.value()); + if (result.space_group.has_value() && anything_indexed) { + std::string line = fmt::format("Space group: {} (No. {})", result.space_group->short_name(), + result.space_group->number); // Name every group the data cannot separate, not just the representative. Some pairs share // their whole absence pattern - an enantiomorphic pair (P4_1 vs P4_3), or I23 vs I2_13 and // I222 vs I2_12_12_1, where the screw condition h00: h=2n is already implied by the @@ -2506,7 +2496,7 @@ static int RunRugnux(int argc, char **argv) { bool named_alternative = false; if (result.space_group_search.has_value()) for (const auto &alt : result.space_group_search->alternatives) - if (alt.number != static_cast(result.space_group_number.value())) { + if (alt.number != result.space_group->number) { line += fmt::format(" or {} (No. {})", alt.short_name(), alt.number); named_alternative = true; } diff --git a/tests/CBORTest.cpp b/tests/CBORTest.cpp index 4b55f6a53..33a677382 100644 --- a/tests/CBORTest.cpp +++ b/tests/CBORTest.cpp @@ -495,6 +495,52 @@ TEST_CASE("CBORSerialize_End", "[CBOR]") { CHECK(output_message.rotation_lattice->GetUnitCell().c == Catch::Approx(60.0)); } +TEST_CASE("CBORSerialize_End_SpaceGroup", "[CBOR]") { + // The determined group travels as both spellings. Only the name keeps the setting, so a + // non-reference setting has to survive the round trip as itself and not as its number's group. + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 1 1 2"); + REQUIRE(sg != nullptr); + REQUIRE(gemmi::find_spacegroup_by_number(sg->number)->xhm() != sg->xhm()); + + std::vector buffer(1024 * 1024); + CBORStream2Serializer serializer(buffer.data(), buffer.size()); + + EndMessage message{}; + message.max_image_number = 0; + message.space_group_name = sg->xhm(); + message.space_group_number = static_cast(sg->number); + REQUIRE_NOTHROW(serializer.SerializeSequenceEnd(message)); + + auto deserialized = CBORStream2Deserialize(buffer.data(), serializer.GetBufferSize()); + REQUIRE(deserialized); + REQUIRE(deserialized->end_message); + const EndMessage &out = *deserialized->end_message; + REQUIRE(out.space_group_name.has_value()); + CHECK(*out.space_group_name == sg->xhm()); + REQUIRE(out.space_group_number.has_value()); + CHECK(*out.space_group_number == static_cast(sg->number)); +} + +TEST_CASE("CBORSerialize_End_SpaceGroupNumberOnly", "[CBOR]") { + // A sender written before the name existed fills in only the number; the message still arrives + // and the missing name is absent rather than empty, so a reader can tell it apart. + std::vector buffer(1024 * 1024); + CBORStream2Serializer serializer(buffer.data(), buffer.size()); + + EndMessage message{}; + message.max_image_number = 0; + message.space_group_number = 96; + REQUIRE_NOTHROW(serializer.SerializeSequenceEnd(message)); + + auto deserialized = CBORStream2Deserialize(buffer.data(), serializer.GetBufferSize()); + REQUIRE(deserialized); + REQUIRE(deserialized->end_message); + const EndMessage &out = *deserialized->end_message; + CHECK_FALSE(out.space_group_name.has_value()); + REQUIRE(out.space_group_number.has_value()); + CHECK(*out.space_group_number == 96); +} + TEST_CASE("CBORSerialize_End_RadIntResult", "[CBOR]") { std::vector buffer(8 * 1024 * 1024); CBORStream2Serializer serializer(buffer.data(), buffer.size()); diff --git a/tests/DiffractionExperimentTest.cpp b/tests/DiffractionExperimentTest.cpp index 24f329001..b38b3c395 100644 --- a/tests/DiffractionExperimentTest.cpp +++ b/tests/DiffractionExperimentTest.cpp @@ -357,6 +357,31 @@ TEST_CASE("DiffractionExperiment_SpaceGroup", "[DiffractionExperiment]") { REQUIRE(x.GetCrystalSystem() == gemmi::CrystalSystem::Tetragonal); REQUIRE(x.GetGemmiSpaceGroup()); REQUIRE(x.GetSpaceGroupName() == "P43212"); + REQUIRE(x.GetSpaceGroupOrP1().number == 96); + + x.SpaceGroupNumber(std::nullopt); + REQUIRE(!x.GetGemmiSpaceGroup()); + REQUIRE(x.GetSpaceGroupOrP1().number == 1); +} + +TEST_CASE("DiffractionExperiment_SpaceGroupSetting", "[DiffractionExperiment]") { + // The experiment carries the group itself, so a non-reference setting survives it. Set by number + // the same groups come back as their reference settings, which is what the number can name. + DiffractionExperiment x; + for (const char *setting : {"P 1 1 2", "I 1 1 2", "R 3:R", "P 2 21 2"}) { + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(setting); + REQUIRE(sg != nullptr); + x.SetSpaceGroup(*sg); + CHECK(x.GetGemmiSpaceGroup()->xhm() == setting); + CHECK(x.GetSpaceGroupNumber() == sg->number); + CHECK(x.GetCentering() == sg->centring_type()); + + x.SpaceGroupNumber(sg->number); + CHECK(x.GetGemmiSpaceGroup()->xhm() != setting); + } + + x.SetSpaceGroup(std::nullopt); + CHECK(!x.GetGemmiSpaceGroup()); } TEST_CASE("DiffractionExperiment_SampleName", "[DiffractionExperiment]") { diff --git a/tests/FrenchWilsonTest.cpp b/tests/FrenchWilsonTest.cpp index dab69ea9b..c2599220f 100644 --- a/tests/FrenchWilsonTest.cpp +++ b/tests/FrenchWilsonTest.cpp @@ -9,6 +9,8 @@ #include "../image_analysis/scale_merge/FrenchWilson.h" namespace { + const gemmi::SpaceGroup &SG(int number) { return *gemmi::find_spacegroup_by_number(number); } + MergedReflection Refl(int h, int k, int l, float d, float I, float sigma) { MergedReflection r; r.h = h; r.k = k; r.l = l; r.d = d; r.I = I; r.sigma = sigma; @@ -28,7 +30,7 @@ namespace { TEST_CASE("French-Wilson: strong reflections reduce to sqrt(I)", "[french_wilson]") { auto v = Background(); v.push_back(Refl(1, 0, 0, 25.0f, 40000.0f, 50.0f)); // I/sigma = 800, clearly strong - ApplyFrenchWilson(v, 1); + ApplyFrenchWilson(v, SG(1)); CHECK(v.back().F == Catch::Approx(std::sqrt(40000.0)).epsilon(0.02)); // ~200 CHECK(v.back().sigmaF >= 0.0f); CHECK(std::isfinite(v.back().sigmaF)); @@ -38,7 +40,7 @@ TEST_CASE("French-Wilson: weak and negative intensities get a positive amplitude auto v = Background(); v.push_back(Refl(2, 0, 0, 20.0f, -40.0f, 50.0f)); // negative measured intensity v.push_back(Refl(3, 0, 0, 15.0f, 10.0f, 50.0f)); // weak, I < sigma - ApplyFrenchWilson(v, 1); + ApplyFrenchWilson(v, SG(1)); const auto& neg = v[v.size() - 2]; const auto& weak = v.back(); CHECK(std::isfinite(neg.F)); @@ -55,7 +57,7 @@ TEST_CASE("French-Wilson: amplitudes are always finite and non-negative", "[fren const float I = (i % 7 == 0) ? -30.0f : static_cast((i % 50) * 40); v.push_back(Refl(1 + i, 2, 3, d, I, 25.0f)); } - ApplyFrenchWilson(v, 96); // P4(3)2(1)2 (has centric reflections + epsilon>1 axes) + ApplyFrenchWilson(v, SG(96)); // P4(3)2(1)2 (has centric reflections + epsilon>1 axes) for (const auto& r : v) { CHECK(std::isfinite(r.F)); CHECK(r.F >= 0.0f); @@ -79,7 +81,7 @@ TEST_CASE("French-Wilson: centric weak reflection gets a smaller amplitude than // acentric. d is set directly, so both share a shell (hence Sigma) regardless of the cell. v.push_back(Refl(3, 1, 0, 5.0f, 0.0f, 10.0f)); v.push_back(Refl(3, 1, 4, 5.0f, 0.0f, 10.0f)); - ApplyFrenchWilson(v, 75); // P4 + ApplyFrenchWilson(v, SG(75)); // P4 const auto& centric = v[v.size() - 2]; const auto& acentric = v.back(); CHECK(centric.F > 0.0f); @@ -91,7 +93,7 @@ TEST_CASE("French-Wilson: unusable sigma falls back to sqrt(max(I,0))", "[french auto v = Background(); v.push_back(Refl(4, 0, 0, 12.0f, 144.0f, NAN)); // no sigma v.push_back(Refl(5, 0, 0, 11.0f, -5.0f, NAN)); // no sigma, negative I - ApplyFrenchWilson(v, 1); + ApplyFrenchWilson(v, SG(1)); CHECK(v[v.size() - 2].F == Catch::Approx(12.0f)); // sqrt(144) CHECK(v.back().F == Catch::Approx(0.0f)); // sqrt(max(-5,0)) } diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index ce35635c3..505105085 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -44,6 +44,76 @@ TEST_CASE("HDF5DataType_ElemType","[HDF5]") { CHECK(type_u32.IsInteger()); } +TEST_CASE("JFJochReader_SpaceGroupSetting", "[HDF5][Full]") { + // A space group is carried through the master file as its name, not its number, because a number + // only ever names the reference setting. Both groups here are non-reference settings that a + // number destroys: "P 1 1 2" comes back from 3 as "P 1 2 1", and "R 3:R" from 146 as "R 3:H" + // (short_name() loses that one too - only xhm() is faithful). + const auto setting = GENERATE(std::string("P 1 1 2"), std::string("R 3:R")); + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(setting); + REQUIRE(sg != nullptr); + + // The number is not a carrier for it, which is the whole reason the name is written. + DiffractionExperiment by_number(DetJF(1)); + by_number.SpaceGroupNumber(sg->number); + CHECK(by_number.GetGemmiSpaceGroup()->xhm() != setting); + + DiffractionExperiment x(DetJF(1)); + x.FilePrefix("test_sg_setting").ImagesPerTrigger(1).OverwriteExistingFiles(true); + x.SetSpaceGroup(*sg); + CHECK(x.GetGemmiSpaceGroup()->xhm() == setting); + CHECK(x.GetSpaceGroupNumber() == sg->number); + + RegisterHDF5Filter(); + { + StartMessage start_message; + x.FillMessage(start_message); + + EndMessage end_message; + end_message.max_image_number = 0; + end_message.space_group_name = sg->xhm(); + std::unique_ptr master = std::make_unique(start_message); + master->Finalize(end_message); + } + { + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("test_sg_setting_master.h5")); + auto dataset = reader.GetDataset(); + REQUIRE(dataset->experiment.GetGemmiSpaceGroup().has_value()); + CHECK(dataset->experiment.GetGemmiSpaceGroup()->xhm() == setting); + } + remove("test_sg_setting_master.h5"); + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} + +TEST_CASE("JFJochReader_SpaceGroupNumberOnly", "[HDF5][Full]") { + // A file written before the name was recorded carries only the number; it still reads back, as + // the reference setting the number names. + DiffractionExperiment x(DetJF(1)); + x.FilePrefix("test_sg_number").ImagesPerTrigger(1).OverwriteExistingFiles(true); + x.SpaceGroupNumber(96); + + RegisterHDF5Filter(); + { + StartMessage start_message; + x.FillMessage(start_message); + + EndMessage end_message; + end_message.max_image_number = 0; + std::unique_ptr master = std::make_unique(start_message); + master->Finalize(end_message); + } + { + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("test_sg_number_master.h5")); + auto dataset = reader.GetDataset(); + REQUIRE(dataset->experiment.GetGemmiSpaceGroup().has_value()); + CHECK(dataset->experiment.GetGemmiSpaceGroup()->xhm() == "P 43 21 2"); + } + remove("test_sg_number_master.h5"); + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} + TEST_CASE("JFJochReader_MasterFile", "[HDF5][Full]") { DiffractionExperiment x(DetJF(1)); diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 45e336243..381491ec5 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -5,14 +5,14 @@ #include "../image_analysis/scale_merge/HKLKey.h" TEST_CASE("HKLKey_NoSG_noMergeFriedel") { - HKLKeyGenerator hkl_key_gen(false, 1); + HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(1)); CHECK(hkl_key_gen(-1, -2, -3) != hkl_key_gen(1,2,3)); CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3)); CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3)); } TEST_CASE("HKLKey_NoSG_MergeFriedel") { - HKLKeyGenerator hkl_key_gen(true, 1); + HKLKeyGenerator hkl_key_gen(true, *gemmi::find_spacegroup_by_number(1)); CHECK(hkl_key_gen(-1, -2, -3) == hkl_key_gen(1,2,3)); CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3)); CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3)); diff --git a/tests/ReindexAmbiguityTest.cpp b/tests/ReindexAmbiguityTest.cpp index f018ceebc..8eaa76522 100644 --- a/tests/ReindexAmbiguityTest.cpp +++ b/tests/ReindexAmbiguityTest.cpp @@ -11,6 +11,8 @@ #include "../image_analysis/scale_merge/ReindexAmbiguity.h" namespace { + const gemmi::SpaceGroup &SG(int number) { return *gemmi::find_spacegroup_by_number(number); } + UnitCell Tetragonal() { return UnitCell{78.0f, 78.0f, 37.0f, 90.0f, 90.0f, 90.0f}; } // A reference set: one distinct intensity per Laue-ASU reflection, so that reflections related by @@ -38,28 +40,28 @@ namespace { TEST_CASE("Reindex: holohedral crystal has no indexing ambiguity", "[reindex]") { // P4(3)2(1)2 (422, holohedral for the tetragonal lattice) -> no twin laws. - CHECK(ReindexAmbiguityOperators(Tetragonal(), 96).empty()); + CHECK(ReindexAmbiguityOperators(Tetragonal(), SG(96)).empty()); // P422 likewise. - CHECK(ReindexAmbiguityOperators(Tetragonal(), 89).empty()); + CHECK(ReindexAmbiguityOperators(Tetragonal(), SG(89)).empty()); } TEST_CASE("Reindex: merohedral crystal exposes the ambiguity operators", "[reindex]") { // P4 (point group 4) in a tetragonal lattice (422) -> a non-trivial reindexing coset. - CHECK_FALSE(ReindexAmbiguityOperators(Tetragonal(), 75).empty()); + CHECK_FALSE(ReindexAmbiguityOperators(Tetragonal(), SG(75)).empty()); } TEST_CASE("Reindex: reference agreement recovers a misindexed dataset", "[reindex]") { const int sg = 75; // P4 const auto reference = DistinctReference(sg); - const auto laws = ReindexAmbiguityOperators(Tetragonal(), sg); + const auto laws = ReindexAmbiguityOperators(Tetragonal(), SG(sg)); REQUIRE_FALSE(laws.empty()); // Deliberately mis-index the data by one twin law. const auto data = ReindexReflections(reference, laws.front()); const auto choice = ChooseReindex( - data, Tetragonal(), sg, - [&](const std::vector &m) { return ReferenceIntensityCC(m, reference, sg); }); + data, Tetragonal(), SG(sg), + [&](const std::vector &m) { return ReferenceIntensityCC(m, reference, SG(sg)); }); CHECK_FALSE(choice.is_identity); // a reindex was needed CHECK(choice.score > 0.99); // the winner realigns with the reference @@ -71,8 +73,8 @@ TEST_CASE("Reindex: a correctly indexed dataset keeps identity", "[reindex]") { const int sg = 75; const auto reference = DistinctReference(sg); const auto choice = ChooseReindex( - reference, Tetragonal(), sg, - [&](const std::vector &m) { return ReferenceIntensityCC(m, reference, sg); }); + reference, Tetragonal(), SG(sg), + [&](const std::vector &m) { return ReferenceIntensityCC(m, reference, SG(sg)); }); CHECK(choice.is_identity); CHECK(choice.score > 0.99); } @@ -80,13 +82,13 @@ TEST_CASE("Reindex: a correctly indexed dataset keeps identity", "[reindex]") { TEST_CASE("Reindex into the ASU: a twin law permutes the reflections without losing any", "[reindex]") { const int sg = 75; // P4 const auto reference = DistinctReference(sg); - const auto laws = ReindexAmbiguityOperators(Tetragonal(), sg); + const auto laws = ReindexAmbiguityOperators(Tetragonal(), SG(sg)); REQUIRE_FALSE(laws.empty()); // Mis-index by a twin law, then reindex back into the ASU: the labels must land where an export // needs them, and the reflection the label carries must be the one the reference has there. const auto misindexed = ReindexReflections(reference, laws.front()); - const auto fixed = ReindexMergedIntoAsu(misindexed, laws.front(), sg, /*merge_friedel=*/true); + const auto fixed = ReindexMergedIntoAsu(misindexed, laws.front(), SG(sg), /*merge_friedel=*/true); REQUIRE(fixed.size() == reference.size()); const HKLKeyGenerator key(true, *gemmi::find_spacegroup_by_number(sg)); @@ -121,7 +123,7 @@ TEST_CASE("Reindex into the ASU: the change of hand swaps the Bijvoet halves", " r.I_plus = 110.0f; r.sigma_plus = 2.0f; r.I_minus = 90.0f; r.sigma_minus = 3.0f; r.F_plus = 10.5f; r.F_minus = 9.5f; - const auto merged = ReindexMergedIntoAsu({r}, sg->change_of_hand_op(), 96, /*merge_friedel=*/true); + const auto merged = ReindexMergedIntoAsu({r}, sg->change_of_hand_op(), SG(96), /*merge_friedel=*/true); REQUIRE(merged.size() == 1); CHECK(merged[0].h == asu.h); CHECK(merged[0].k == asu.k); @@ -135,7 +137,7 @@ TEST_CASE("Reindex into the ASU: the change of hand swaps the Bijvoet halves", " CHECK(merged[0].F_minus == 10.5f); // With the mates kept apart the merge stores the minus hand at -hkl, so the row moves there. - const auto anom = ReindexMergedIntoAsu({r}, sg->change_of_hand_op(), 96, /*merge_friedel=*/false); + const auto anom = ReindexMergedIntoAsu({r}, sg->change_of_hand_op(), SG(96), /*merge_friedel=*/false); REQUIRE(anom.size() == 1); CHECK(anom[0].h == -asu.h); CHECK(anom[0].k == -asu.k); @@ -168,7 +170,7 @@ TEST_CASE("Reindex into the ASU: both mates of an anomalous pair follow the same MergedReflection minus = plus; // same anomalous split, stored at -hkl_asu minus.h = -asu.h; minus.k = -asu.k; minus.l = -asu.l; - const auto out = ReindexMergedIntoAsu({plus, minus}, op, 75, /*merge_friedel=*/false); + const auto out = ReindexMergedIntoAsu({plus, minus}, op, SG(75), /*merge_friedel=*/false); REQUIRE(out.size() == 2); // The two rows are still one pair: same Laue label up to the Friedel sign, opposite hands. diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index 9271413ac..bc27855c1 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -59,7 +59,7 @@ TEST_CASE("ResultReport_Render", "[Diagnostics]") { result.indexing_rate = 0.87f; result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f, .alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f}; - result.space_group_number = 96; + result.space_group = *gemmi::find_spacegroup_by_number(96); result.used_beam_x_pxl = 766.62f; result.used_beam_y_pxl = 846.87f; result.used_distance_mm = 243.53f; diff --git a/tests/RfreeFlagsTest.cpp b/tests/RfreeFlagsTest.cpp index 320b6abb9..bb78a8f8b 100644 --- a/tests/RfreeFlagsTest.cpp +++ b/tests/RfreeFlagsTest.cpp @@ -10,6 +10,8 @@ #include "../image_analysis/scale_merge/RfreeFlags.h" namespace { + const gemmi::SpaceGroup &SG(int number) { return *gemmi::find_spacegroup_by_number(number); } + MergedReflection Refl(int h, int k, int l, float d) { MergedReflection r; r.h = h; r.k = k; r.l = l; r.d = d; @@ -37,8 +39,8 @@ TEST_CASE("R-free flags are deterministic and hit the requested fraction", "[rfr auto a = Grid(-15, 15); auto b = a; // Floor off (min_free = 0) so this isolates the pure-hash fraction on this modest grid. - AssignRfreeFlags(a, 1, 0.05, /*min_free=*/0); - AssignRfreeFlags(b, 1, 0.05, /*min_free=*/0); + AssignRfreeFlags(a, SG(1), 0.05, /*min_free=*/0); + AssignRfreeFlags(b, SG(1), 0.05, /*min_free=*/0); REQUIRE(a.size() == b.size()); for (size_t i = 0; i < a.size(); ++i) CHECK(a[i].rfree_flag == b[i].rfree_flag); // pure function of the reflection @@ -58,7 +60,7 @@ TEST_CASE("R-free flags never split a Friedel/Bijvoet pair", "[rfree]") { v.push_back(Refl(h, k, l, d)); v.push_back(Refl(-h, -k, -l, d)); } - AssignRfreeFlags(v, 1, 0.10); // P1 -> only Friedel relates the mates + AssignRfreeFlags(v, SG(1), 0.10); // P1 -> only Friedel relates the mates std::map, bool> flag; for (const auto& r : v) flag[{r.h, r.k, r.l}] = r.rfree_flag; @@ -80,7 +82,7 @@ TEST_CASE("R-free flags are shared across symmetry equivalents", "[rfree]") { if (h == 0 && k == 0 && l == 0) continue; v.push_back(Refl(h, k, l, 60.0f / (1 + h * h + k * k + l * l))); } - AssignRfreeFlags(v, 75, 0.10); // P4 + AssignRfreeFlags(v, SG(75), 0.10); // P4 std::map, bool> flag; for (const auto& r : v) flag[{r.h, r.k, r.l}] = r.rfree_flag; @@ -101,7 +103,7 @@ TEST_CASE("R-free flags spread across resolution", "[rfree]") { v.push_back(Refl(2, 1 + i, 3, 4.0f)); // mid res v.push_back(Refl(2, 3, 1 + i, 2.0f)); // high res } - AssignRfreeFlags(v, 1, 0.10); + AssignRfreeFlags(v, SG(1), 0.10); int lo = 0, mid = 0, hi = 0; for (const auto& r : v) { if (!r.rfree_flag) continue; @@ -124,8 +126,8 @@ TEST_CASE("R-free flags are identical across datasets of one crystal form", "[rf narrow.push_back(r); // Floor off so the two different-sized sets share one effective fraction (the floor is the only // thing that ties the fraction to the dataset; with it off this is the pure per-hkl guarantee). - AssignRfreeFlags(wide, 96, 0.05, /*min_free=*/0); - AssignRfreeFlags(narrow, 96, 0.05, /*min_free=*/0); + AssignRfreeFlags(wide, SG(96), 0.05, /*min_free=*/0); + AssignRfreeFlags(narrow, SG(96), 0.05, /*min_free=*/0); std::map, bool> flag; for (const auto& r : wide) flag[{r.h, r.k, r.l}] = r.rfree_flag; @@ -143,13 +145,13 @@ TEST_CASE("ApplyReferenceFreeFlags imports the reference test set", "[rfree]") { // A reference with its own free set; a dataset that starts from the per-hkl hash must, after the // import, carry exactly the reference's flags on every reflection they share. auto reference = Grid(-12, 12); - AssignRfreeFlags(reference, 96, 0.07); + AssignRfreeFlags(reference, SG(96), 0.07); std::map, bool> ref_flag; for (const auto& r : reference) ref_flag[{r.h, r.k, r.l}] = r.rfree_flag; auto data = Grid(-12, 12); - AssignRfreeFlags(data, 96, 0.30); // deliberately a different fraction/hash split - const size_t matched = ApplyReferenceFreeFlags(data, 96, reference); + AssignRfreeFlags(data, SG(96), 0.30); // deliberately a different fraction/hash split + const size_t matched = ApplyReferenceFreeFlags(data, SG(96), reference); CHECK(matched == data.size()); for (const auto& r : data) CHECK(r.rfree_flag == ref_flag[{r.h, r.k, r.l}]); @@ -166,9 +168,9 @@ TEST_CASE("R-free flags floor the test set size on small data", "[rfree]") { const double n = static_cast(v.size()); auto plain = v, floored = v, capped = v; - AssignRfreeFlags(plain, 1, 0.05, /*min_free=*/0); // pure 5% ~= 37 free - AssignRfreeFlags(floored, 1, 0.05, /*min_free=*/60); // 60/738 = 8.1% (above 5%, under the 10% cap) - AssignRfreeFlags(capped, 1, 0.05, /*min_free=*/100000); // floor wants ~all; capped near 10% + AssignRfreeFlags(plain, SG(1), 0.05, /*min_free=*/0); // pure 5% ~= 37 free + AssignRfreeFlags(floored, SG(1), 0.05, /*min_free=*/60); // 60/738 = 8.1% (above 5%, under the 10% cap) + AssignRfreeFlags(capped, SG(1), 0.05, /*min_free=*/100000); // floor wants ~all; capped near 10% CHECK(FreeFraction(floored) > FreeFraction(plain)); // the floor lifted the test set CHECK(FreeFraction(floored) > 0.06); // ~8%, clearly above the bare 5% @@ -177,6 +179,6 @@ TEST_CASE("R-free flags floor the test set size on small data", "[rfree]") { TEST_CASE("R-free fraction of zero flags nothing", "[rfree]") { auto v = Grid(1, 6); - AssignRfreeFlags(v, 1, 0.0); + AssignRfreeFlags(v, SG(1), 0.0); for (const auto& r : v) CHECK(!r.rfree_flag); } diff --git a/viewer/windows/JFJochMergeStatsWindow.cpp b/viewer/windows/JFJochMergeStatsWindow.cpp index 7dd4f7efb..520d6b964 100644 --- a/viewer/windows/JFJochMergeStatsWindow.cpp +++ b/viewer/windows/JFJochMergeStatsWindow.cpp @@ -24,8 +24,6 @@ #include #include -#include "gemmi/symmetry.hpp" // find_spacegroup_by_number - #include "../charts/JFJochSimpleChartView.h" #include "../widgets/ToolbarIcons.h" @@ -71,7 +69,7 @@ namespace { JFJochMergeStatsWindow::JFJochMergeStatsWindow(const QString &title, const MergeStatistics &stats, double isa, bool has_reference, const TwinningAnalysisResult &twinning, - std::optional space_group_number, + const std::optional &space_group, const std::optional &space_group_search, const std::vector> &merged_i_sigma, QWidget *parent) @@ -111,12 +109,8 @@ JFJochMergeStatsWindow::JFJochMergeStatsWindow(const QString &title, const Merge } if (has_reference_) hero->addWidget(MakeCard(pct(o.cc_ref * 100.0), "CCref", this)); - if (space_group_number.has_value()) { - const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(static_cast(*space_group_number)); - const QString sg_name = sg ? QString::fromStdString(sg->short_name()) - : QString::number(*space_group_number); - hero->addWidget(MakeCard(sg_name, "Space group", this)); - } + if (space_group.has_value()) + hero->addWidget(MakeCard(QString::fromStdString(space_group->short_name()), "Space group", this)); layout->addLayout(hero); // De-novo space-group search: the point group + the ranked candidates it scored (only present when diff --git a/viewer/windows/JFJochMergeStatsWindow.h b/viewer/windows/JFJochMergeStatsWindow.h index c79d75494..e8402121b 100644 --- a/viewer/windows/JFJochMergeStatsWindow.h +++ b/viewer/windows/JFJochMergeStatsWindow.h @@ -31,7 +31,7 @@ public: JFJochMergeStatsWindow(const QString &title, const MergeStatistics &stats, double isa, bool has_reference, const TwinningAnalysisResult &twinning, - std::optional space_group_number, + const std::optional &space_group, const std::optional &space_group_search, const std::vector> &merged_i_sigma, QWidget *parent = nullptr); diff --git a/viewer/windows/JFJochProcessingJobsWindow.cpp b/viewer/windows/JFJochProcessingJobsWindow.cpp index 57ba45a2f..2c7e0cb44 100644 --- a/viewer/windows/JFJochProcessingJobsWindow.cpp +++ b/viewer/windows/JFJochProcessingJobsWindow.cpp @@ -591,7 +591,7 @@ void JFJochProcessingJobsWindow::showStats(const QString &id) { } if (j.id == id && j.has_merge_stats) { auto *win = new JFJochMergeStatsWindow(j.label, j.merge_stats, j.isa, j.merge_has_reference, - j.twinning, j.space_group_number, j.space_group_search, + j.twinning, j.space_group, j.space_group_search, j.merged_i_sigma, window()); win->show(); return; @@ -762,7 +762,7 @@ void JFJochProcessingJobsWindow::onFinished(ProcessResult result) { jobs_[row].isa = result.error_model_isa; jobs_[row].merge_has_reference = result.has_reference; jobs_[row].twinning = result.twinning; - jobs_[row].space_group_number = result.space_group_number; + jobs_[row].space_group = result.space_group; jobs_[row].space_group_search = result.space_group_search; jobs_[row].merged_i_sigma = result.merged_i_sigma; if (jobs_[row].graph_btn) diff --git a/viewer/windows/JFJochProcessingJobsWindow.h b/viewer/windows/JFJochProcessingJobsWindow.h index 64331c0c5..5d33fa5d3 100644 --- a/viewer/windows/JFJochProcessingJobsWindow.h +++ b/viewer/windows/JFJochProcessingJobsWindow.h @@ -69,7 +69,7 @@ private: double isa = 0.0; bool merge_has_reference = false; TwinningAnalysisResult twinning; // twinning test of the merged intensities - std::optional space_group_number; // final space group (searched or fixed) + std::optional space_group; // final space group (searched or fixed) std::optional space_group_search; // ranked candidates, when a search ran std::vector> merged_i_sigma; // ISa diagnostic scatter // Detector calibration result (Calibration mode): the fit, the header geometry it is compared diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 89a5be7ac..89db926c1 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -742,14 +742,20 @@ void NXmx::Sample(const StartMessage &start, const EndMessage &end) { group.SaveScalar("name", start.sample_name); // The offline analysis determines the space group only after merging, so it arrives on the end - // message; prefer it over the (usually empty) start-message value. - const auto space_group_number = end.space_group_number ? end.space_group_number - : start.space_group_number; - if (space_group_number) { - group.SaveScalar("space_group_number", space_group_number.value()); - auto *sg = gemmi::find_spacegroup_by_number(space_group_number.value()); - if (sg != nullptr) - group.SaveScalar("space_group", sg->short_name()); + // message; prefer it over the (usually empty) start-message value the user supplied. Within the + // end message the NAME is preferred and the number is the fallback, because only the name keeps + // the setting - a number always reads back as the reference setting - and a sender written + // before the name existed fills in only the number. + const gemmi::SpaceGroup *sg = nullptr; + if (end.space_group_name) + sg = gemmi::find_spacegroup_by_name(*end.space_group_name); + if (sg == nullptr && end.space_group_number) + sg = gemmi::find_spacegroup_by_number(static_cast(*end.space_group_number)); + if (sg == nullptr && start.space_group_number) + sg = gemmi::find_spacegroup_by_number(static_cast(*start.space_group_number)); + if (sg != nullptr) { + group.SaveScalar("space_group_number", static_cast(sg->number)); + group.SaveScalar("space_group", sg->xhm()); } std::optional unit_cell;