diff --git a/image_analysis/scale_merge/ReindexAmbiguity.cpp b/image_analysis/scale_merge/ReindexAmbiguity.cpp index 46753047..e72f6ed2 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.cpp +++ b/image_analysis/scale_merge/ReindexAmbiguity.cpp @@ -5,6 +5,7 @@ #include #include +#include #include "gemmi/twin.hpp" @@ -63,6 +64,34 @@ std::vector ReindexReflections(const std::vector ReindexMergedIntoAsu(const std::vector &merged, + const gemmi::Op &op, int space_group_number, + 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); + + std::vector out = merged; + for (auto &r : out) { + const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}}); + const HKLKey key = key_gen(h[0], h[1], h[2]); + if (!key.plus) { + std::swap(r.I_plus, r.I_minus); + std::swap(r.sigma_plus, r.sigma_minus); + std::swap(r.F_plus, r.F_minus); + std::swap(r.sigmaF_plus, r.sigmaF_minus); + } + const bool at_asu_index = merge_friedel || key.plus; + r.h = at_asu_index ? key.h : -key.h; + r.k = at_asu_index ? key.k : -key.k; + r.l = at_asu_index ? key.l : -key.l; + } + return out; +} + ReindexChoice ChooseReindex(const std::vector &merged, const UnitCell &cell, int space_group_number, const std::function &)> &score, diff --git a/image_analysis/scale_merge/ReindexAmbiguity.h b/image_analysis/scale_merge/ReindexAmbiguity.h index 93f44ce5..354e5598 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.h +++ b/image_analysis/scale_merge/ReindexAmbiguity.h @@ -29,9 +29,19 @@ std::vector ReindexAmbiguityOperators(const UnitCell &cell, int space double max_obliquity_deg = 2.0); // A copy of `merged` with each Miller index transformed by `op` (the intensities are unchanged). +// The indices are left where `op` puts them, which is what the scoring callers want - they reduce to +// the ASU themselves. 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` - +// 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, + bool merge_friedel); + struct ReindexChoice { gemmi::Op op = gemmi::Op::identity(); // the winning operator (identity = keep the current indexing) bool is_identity = true; diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 06b24f45..9fcdd4c3 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -3,9 +3,11 @@ #include "ModelValidation.h" +#include #include #include #include +#include #include #include @@ -72,6 +74,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & logger.Error("Model validation: model {} has no usable space group", model_path); return result; } + result.model_space_group_number = sg->number; // If the data was indexed in the enantiomorph of the model's space group (e.g. data P4(1)2(1)2, // model P4(3)2(1)2 - the merged intensities cannot tell them apart), reindex the observed @@ -88,6 +91,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & if (enant && enant->number == sg->number) { reindexed = ReindexReflections(merged, dsg->change_of_hand_op()); obs_ptr = &reindexed; + result.hand_op = dsg->change_of_hand_op(); logger.Info("Model validation: data space group {} is the enantiomorph of the model {}; " "reindexed the observed reflections into the model's hand", dsg->short_name(), sg->hm); @@ -227,16 +231,22 @@ ModelValidationResult ValidateAgainstModel(const std::vector & Fit best = fit_model(obs); if (probe_indexing_ambiguity) { const auto reindex_ops = ReindexAmbiguityOperators(cell, sg->number); - bool did_reindex = false; + 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)); - if (cand.r_free < best.r_free) { best = std::move(cand); did_reindex = true; } + candidate_r_free.push_back(cand.r_free); + if (cand.r_free < best.r_free) { best = std::move(cand); result.indexing_op = op; } } - if (!reindex_ops.empty()) + if (!reindex_ops.empty()) { + // The runner-up as well as the winner: the margin between them is what says whether the + // choice was made on evidence, and on weak data the two can come out within noise. + std::sort(candidate_r_free.begin(), candidate_r_free.end()); logger.Info("Model validation: no reference - probed {} indexing solution(s) against the model; " - "{} (R-free {:.4f})", reindex_ops.size() + 1, - did_reindex ? "reindexed to the lower-R-free solution" : "kept the current indexing", - best.r_free); + "{} (R-free {:.4f}, runner-up {:.4f})", candidate_r_free.size(), + result.indexing_op == gemmi::Op::identity() ? "kept the current indexing" + : "reindexed to the lower-R-free solution", + candidate_r_free[0], candidate_r_free[1]); + } } gemmi::AsuData> &fmodel = best.fmodel; @@ -322,3 +332,145 @@ ModelValidationResult ValidateAgainstModel(const std::vector & output_prefix, output_prefix, output_prefix); return result; } + +namespace { + +// The reindexing operator as it reads on Miller indices ("k,h,-l" rather than "y,x,-z"). +std::string hkl_triplet(const gemmi::Op &op) { + std::string t = op.triplet(); + std::replace(t.begin(), t.end(), 'x', 'h'); + std::replace(t.begin(), t.end(), 'y', 'k'); + std::replace(t.begin(), t.end(), 'z', 'l'); + return t; +} + +} // 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; + if (!validation.ok) + return space_group_number; + + // The change of hand also changes the space group the reflections are in: P4(1)2(1)2 data written + // in the model's hand are P4(3)2(1)2 data. The merged intensities are the same measurements either + // way - the enantiomorph is exactly what they cannot distinguish - so the model's hand is the only + // evidence there is, and adopting it is what makes the file refinable against that model. + if (!(validation.hand_op == gemmi::Op::identity()) && validation.model_space_group_number > 0) { + merged = ReindexMergedIntoAsu(merged, validation.hand_op, + validation.model_space_group_number, merge_friedel); + space_group_number = validation.model_space_group_number; + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); + logger.Info("Model validation: the written reflections take the model's enantiomorph, {} ({}), " + "reindexed by {}", sg ? sg->short_name() : "?", space_group_number, + hkl_triplet(validation.hand_op)); + } + + // 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); + logger.Info("Model validation: the written reflections take the model's indexing, reindexed by {}", + hkl_triplet(validation.indexing_op)); + } + return space_group_number; +} + +std::vector ModelReferenceIntensities(const std::string &model_path, + const std::optional &cell, + std::optional space_group_number, + double d_min, + Logger &logger) { + std::vector out; + if (!(d_min > 0.0)) { + logger.Warning("Model reference: no resolution limit to compute the model intensities to"); + return out; + } + + gemmi::Structure st; + try { + st = gemmi::read_pdb(gemmi::MaybeGzipped(model_path)); + } catch (const std::exception &e) { + logger.Error("Model reference: cannot read model {}: {}", model_path, e.what()); + return out; + } + if (st.models.empty() || !st.cell.is_crystal()) { + logger.Error("Model reference: model {} has no atoms or no unit cell", model_path); + return out; + } + + // Put the model in the cell and group the run works in, where it knows them, so the reference is + // indexed the way the data are. The correlation that consumes this matches on hkl, so a small cell + // difference costs nothing; the space group is what has to agree. + if (cell.has_value()) { + const gemmi::UnitCell target = *cell; + if (target.is_crystal()) { + const gemmi::UnitCell old = st.cell; + for (gemmi::Model &m : st.models) + for (gemmi::Chain &ch : m.chains) + for (gemmi::Residue &r : ch.residues) + for (gemmi::Atom &a : r.atoms) + a.pos = target.orthogonalize(old.fractionalize(a.pos)); + 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(); + const gemmi::SpaceGroup *sg = st.find_spacegroup(); + if (!sg) { + logger.Error("Model reference: model {} has no usable space group", model_path); + return out; + } + st.setup_cell_images(); + + gemmi::DensityCalculator dc; + dc.d_min = d_min; + dc.rate = 1.5; + dc.set_grid_cell_and_spacegroup(st); + dc.set_refmac_compatible_blur(st.models[0]); + dc.put_model_density_on_grid(st.models[0]); + gemmi::AsuData> fcalc = + gemmi::transform_map_to_f_phi(dc.grid, true).prepare_asu_data(dc.d_min, dc.blur, false, false, false); + + // Flat bulk solvent at the standard constants. Nothing here is fitted - there are no observations + // yet - but without it the few lowest-resolution reflections are the largest and the most wrong, + // and a correlation on raw intensities would be led by them. + constexpr double K_SOL = 0.35; + constexpr double B_SOL = 46.0; + gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac); + gemmi::Grid mask_grid; + mask_grid.unit_cell = dc.grid.unit_cell; + mask_grid.spacegroup = dc.grid.spacegroup; + mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up); + masker.put_mask_on_grid(mask_grid, st.models[0]); + gemmi::AsuData> fmask = + gemmi::transform_map_to_f_phi(mask_grid, true).prepare_asu_data(dc.d_min, 0); + + std::unordered_map> mask_by_hkl; + mask_by_hkl.reserve(fmask.v.size()); + for (const auto &hv : fmask.v) + mask_by_hkl[hkl_key(hv.hkl)] = hv.value; + + const gemmi::UnitCell &ucell = st.cell; + out.reserve(fcalc.v.size()); + for (const auto &hv : fcalc.v) { + const double d = ucell.calculate_d(hv.hkl); + if (!(d > 0.0)) + continue; + std::complex f = hv.value; + const auto it = mask_by_hkl.find(hkl_key(hv.hkl)); + if (it != mask_by_hkl.end()) + f += static_cast(K_SOL * std::exp(-B_SOL / (4.0 * d * d))) * it->second; + const double F = std::abs(f); + out.push_back(MergedReflection{.h = hv.hkl[0], .k = hv.hkl[1], .l = hv.hkl[2], + .I = static_cast(F * F), + .d = static_cast(d)}); + } + + logger.Info("Model reference: {} intensities computed from {} to {:.2f} A, space group {}", + out.size(), model_path, d_min, sg->short_name()); + return out; +} diff --git a/rugnux/ModelValidation.h b/rugnux/ModelValidation.h index 34da9bf8..abd5096c 100644 --- a/rugnux/ModelValidation.h +++ b/rugnux/ModelValidation.h @@ -7,6 +7,8 @@ #include #include +#include "gemmi/symmetry.hpp" + #include "../common/Reflection.h" // MergedReflection #include "../common/UnitCell.h" @@ -24,6 +26,15 @@ struct ModelValidationResult { int n_free = 0; double mean_atom_density_sigma = 0.0; // mean 2Fo-Fc value at atom centres, in sigma std::string maps_prefix; // where the .ccp4 / _maps.mtz were written + + // What had to be done to the observed reflections to bring them into the model's frame, in this + // order: the change of hand, where the data were merged in the model's enantiomorph, and then the + // alternative-indexing operator picked by R-free. Identity where nothing was needed. AdoptModelFrame + // below applies them to the reflections that are written out, so the file, the R-factors and the + // maps all describe one indexing. + gemmi::Op hand_op = gemmi::Op::identity(); + gemmi::Op indexing_op = gemmi::Op::identity(); + int model_space_group_number = 0; }; // Given merged intensities and a PDB atomic model, scale the model structure factors (with an @@ -50,3 +61,25 @@ ModelValidationResult ValidateAgainstModel(const std::vector & Logger &logger, std::optional data_space_group_number = std::nullopt, 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 the change of hand was adopted, and `data_space_group_number` 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); + +// 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 +// 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, + double d_min, + Logger &logger); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index cc33a592..4073553a 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1398,6 +1398,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b std::unique_ptr writer; std::unique_ptr writer_queue; + // Whether anything is in place to break the indexing ambiguity, for the warning further down: a + // reference MTZ, or (serial stills) the reference computed from a model just below. + bool ambiguity_reference = !config_.reference_data.empty(); if (write_files && config_.write_process_h5) { fill_writer_maps(); writer = std::make_unique(start_message, /*check_overwrite_at_start=*/true, @@ -1438,8 +1441,26 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // With no per-image file to write, nothing reads the message's reflection list. if (!writer_queue) indexer->KeepReflectionsInMessage(false); - if (!config_.reference_data.empty()) + // The reference that breaks the per-image indexing ambiguity: a reference MTZ where there is + // one, otherwise intensities computed from the model, which carry the same information. Serial + // stills need it here or not at all - each crystal is indexed in its own hand, so by the time + // there is a merge for --model to fit, the hands have already been averaged together. Rotation + // data have one lattice and one hand, so their ambiguity is settled at the merge instead. + std::vector model_reference; + 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()) { + 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); + if (!model_reference.empty()) { + indexer->ReferenceIntensities(model_reference); + ambiguity_reference = true; + } + } // Second pass of the rotation two-pass: predict with the smoothed mosaicity fitted in the pre-pass. if (!geometry_prepass && !prepass_mosaicity_.empty()) indexer->SetPredictionMosaicityOverride(prepass_mosaicity_); @@ -3059,19 +3080,27 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b std::replace(t.begin(), t.end(), 'z', 'l'); laws += (laws.empty() ? "" : " ; ") + t; } - if (config_.reference_data.empty()) { + // A model resolves it too: on rotation data at the validation step below, by R-free; + // on serial stills only through the per-image reference computed from it at integration + // time, which needs the cell and the space group up front - so a model on its own does + // not always suffice there, and `ambiguity_reference` says whether it did. + const bool will_be_resolved = ambiguity_reference + || (!config_.model_path.empty() && experiment_.IsRotationIndexing()); + if (!will_be_resolved) { const std::string msg = fmt::format( "Indexing ambiguity: this cell / space group admits alternative indexing " "(reindex operator(s): {}). {} are indexed in one hand at random, and rugnux can only " "break this against an external reference. WITHOUT one the merge mixes the hands and " - "CC1/2 is degraded - supply a reference (-z reference.mtz or --model) to resolve it.", - laws, experiment_.IsRotationIndexing() ? "Lattices" : "Serial-stills crystals"); + "CC1/2 is degraded - supply a reference MTZ (-z) or a model ({}) to resolve it.", + laws, experiment_.IsRotationIndexing() ? "Lattices" : "Serial-stills crystals", + experiment_.IsRotationIndexing() ? "--model" : "--model, which needs -C and -S here"); logger.Warning("{}", msg); stats_text << " !! " << msg << "\n\n"; result.warnings.push_back(msg); } else { logger.Info("Indexing ambiguity present (reindex operator(s): {}); resolved against the " - "supplied reference.", laws); + "supplied {}.", laws, + config_.reference_data.empty() ? "model" : "reference"); } } } @@ -3243,6 +3272,46 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } } + // Model validation runs BEFORE the reflection files are written, because it is what settles the + // frame they are written in: the enantiomorph, which merged intensities cannot choose, and - + // with no reference MTZ - the alternative indexing. Both are relabelings of the same + // measurements, so the statistics computed above are unaffected; what changes is that the file, + // 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(); + // 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, + /*probe_indexing_ambiguity=*/config_.reference_data.empty()); + 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)) { + // 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; + } + // 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 + // indexings of one dataset. Same order as the merged set: hand first, then indexing. + for (const gemmi::Op &op : {validation.hand_op, validation.indexing_op}) + if (!(op == gemmi::Op::identity())) + for (auto &io : indexer->GetIntegrationOutcome()) + for (auto &r : io.reflections) { + const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}}); + r.h = h[0]; r.k = h[1]; r.l = h[2]; + } + } + } + // Not on the geometry pre-pass. Pass 1 exists to choose the space group and post-refine the // geometry; pass 2 remakes these files seconds later at the refined geometry, and that is the // answer anyone reads. Writing them twice costs 0.6 s of a 15 s run - a fifth of the merged @@ -3265,17 +3334,6 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // partials, which carry the first-pass per-image scale. ScalingResult(indexer->GetIntegrationOutcome()).SaveToFile(config_.output_prefix); } - - if (result.consensus_cell && write_files && !config_.model_path.empty() && !geometry_prepass) { - phase("Validating against model"); - const auto data_sg = experiment_.GetSpaceGroupNumber(); - // 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. - ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path, - config_.output_prefix, logger, - data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt, - /*probe_indexing_ambiguity=*/config_.reference_data.empty()); - } } // Unmerged observations (--export-unmerged). After the merge, because that is where the space diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 839865d9..0ed8fd15 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1620,6 +1620,35 @@ static int RunRugnux(int argc, char **argv) { std::cout << AnisotropyToText(merged_statistics.anisotropy) << std::endl; } + // Before the reflection files, as in the full pipeline: the model settles the enantiomorph and, + // with no reference MTZ, the alternative indexing, and the files are written in what it settled. + if (!output_prefix.empty() && !model_pdb.empty()) { + const auto data_sg = experiment.GetSpaceGroupNumber(); + // 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, + /*probe_indexing_ambiguity=*/reference_data.empty()); + 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); + // The unmerged export below writes the integrated observations, so it follows the same + // relabeling as the merged set - hand first, then indexing - or the two files would + // describe different indexings of one dataset. + for (const gemmi::Op &op : {validation.hand_op, validation.indexing_op}) + if (!(op == gemmi::Op::identity())) + for (auto &io : reflections) + for (auto &r : io.reflections) { + const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}}); + r.h = h[0]; r.k = h[1]; r.l = h[2]; + } + } + } + // Unmerged observations (--export-unmerged), from the integrated observations rather than the // merged ones: the partiality and the per-image scale are left for the reading program, which // fits a scale model of its own. @@ -1645,16 +1674,6 @@ static int RunRugnux(int argc, char **argv) { em_report, twinning, output_prefix, static_cast(nthreads)); } - if (!output_prefix.empty() && !model_pdb.empty()) { - const auto data_sg = experiment.GetSpaceGroupNumber(); - // 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. - ValidateAgainstModel(merged_reflections, *experiment.GetUnitCell(), model_pdb, - output_prefix, logger, - data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt, - /*probe_indexing_ambiguity=*/reference_data.empty()); - } - // --mode scale re-merges stored reflections, so it determines a merging result and gets the // same report. It has no ProcessResult of its own; fill the parts that exist here. ProcessResult scale_result; @@ -2306,10 +2325,17 @@ static int RunRugnux(int argc, char **argv) { // I-centering - so the choice among them is a convention, not a measurement. The search // reports the representative as the lowest space-group number; saying so here keeps the // summary from claiming a decision the diffraction did not make. + // Not the adopted group itself: where --model has taken the enantiomorph, the group now + // reported is the one the search listed as the alternative, and naming it as its own + // alternative would read as a decision that was not made. + bool named_alternative = false; if (result.space_group_search.has_value()) for (const auto &alt : result.space_group_search->alternatives) - line += fmt::format(" or {} (No. {})", alt.short_name(), alt.number); - if (result.space_group_search.has_value() && !result.space_group_search->alternatives.empty()) + if (alt.number != static_cast(result.space_group_number.value())) { + line += fmt::format(" or {} (No. {})", alt.short_name(), alt.number); + named_alternative = true; + } + if (named_alternative) line += " - indistinguishable from these data"; std::cout << line << std::endl; } diff --git a/tests/ReindexAmbiguityTest.cpp b/tests/ReindexAmbiguityTest.cpp index bfc7953f..c3a5598c 100644 --- a/tests/ReindexAmbiguityTest.cpp +++ b/tests/ReindexAmbiguityTest.cpp @@ -3,6 +3,7 @@ #include +#include #include #include @@ -75,3 +76,69 @@ TEST_CASE("Reindex: a correctly indexed dataset keeps identity", "[reindex]") { CHECK(choice.is_identity); CHECK(choice.score > 0.99); } + +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); + 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); + REQUIRE(fixed.size() == reference.size()); + + const HKLKeyGenerator key(true, *gemmi::find_spacegroup_by_number(sg)); + std::unordered_map ref_by_key; + for (const auto &r : reference) + ref_by_key[key(r).pack()] = r.I; + + std::unordered_set seen; + for (const auto &r : fixed) { + const HKLKey k = key(r); + CHECK(k.h == r.h); // already at its own ASU index + CHECK(k.k == r.k); + CHECK(k.l == r.l); + CHECK(seen.insert(k.pack()).second); // no two reflections land on one label + const auto it = ref_by_key.find(k.pack()); + REQUIRE(it != ref_by_key.end()); + CHECK(it->second == r.I); + } +} + +TEST_CASE("Reindex into the ASU: the change of hand swaps the Bijvoet halves", "[reindex]") { + // The change of hand is the inversion, so it leaves the Laue-ASU label alone and moves the + // anomalous signal instead - which is the whole of what adopting a model's hand does to intensities. + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(96); // P4(3)2(1)2 + REQUIRE(sg != nullptr); + const HKLKeyGenerator key(true, *sg); + const HKLKey asu = key(3, 1, 2); + + MergedReflection r; + r.h = asu.h; r.k = asu.k; r.l = asu.l; + r.I = 100.0f; r.sigma = 1.0f; r.d = 10.0f; + 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); + REQUIRE(merged.size() == 1); + CHECK(merged[0].h == asu.h); + CHECK(merged[0].k == asu.k); + CHECK(merged[0].l == asu.l); + CHECK(merged[0].I == 100.0f); + CHECK(merged[0].I_plus == 90.0f); + CHECK(merged[0].I_minus == 110.0f); + CHECK(merged[0].sigma_plus == 3.0f); + CHECK(merged[0].sigma_minus == 2.0f); + CHECK(merged[0].F_plus == 9.5f); + 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); + REQUIRE(anom.size() == 1); + CHECK(anom[0].h == -asu.h); + CHECK(anom[0].k == -asu.k); + CHECK(anom[0].l == -asu.l); + CHECK(anom[0].I == 100.0f); +}