diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index d4b81410..f5ca5d0a 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -21,6 +21,7 @@ #include // Mtz (map-coefficient output) #include "../common/Logger.h" +#include "../image_analysis/scale_merge/ReindexAmbiguity.h" // ReindexReflections namespace { @@ -50,7 +51,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector & const UnitCell &cell, const std::string &model_path, const std::string &output_prefix, - Logger &logger) { + Logger &logger, + std::optional data_space_group_number) { ModelValidationResult result; // --- read the atomic model --- @@ -71,9 +73,32 @@ ModelValidationResult ValidateAgainstModel(const std::vector & return result; } + // 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 + // reflections into the model's hand so the two settings agree. This does not change the + // R-factors (which use |F|), but keeps the observed data consistent with the model. + std::vector reindexed; + const std::vector *obs_ptr = &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()); + const gemmi::SpaceGroup *enant = gemmi::find_spacegroup_by_ops(eops); + if (enant && enant->number == sg->number) { + reindexed = ReindexReflections(merged, dsg->change_of_hand_op()); + obs_ptr = &reindexed; + 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); + } + } + } + const std::vector &obs = *obs_ptr; + // Resolution limit from the data (the merged set is already resolution-trimmed). double d_min = 0.0; - for (const MergedReflection &r : merged) + for (const MergedReflection &r : obs) if (r.d > 0 && (d_min == 0.0 || r.d < d_min)) d_min = r.d; if (d_min <= 0.0) { @@ -130,7 +155,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // Observed amplitudes are the French-Wilson |F| already computed at the end of the merge // (MergedReflection.F), so the model R-free / maps use exactly the same amplitudes as the // written reflection file. - for (const MergedReflection &r : merged) { + for (const MergedReflection &r : obs) { if (std::isnan(r.F)) continue; gemmi::Miller h{{r.h, r.k, r.l}}; if (!asu.is_in(h)) h = asu.to_asu(h, gops).first; diff --git a/rugnux/ModelValidation.h b/rugnux/ModelValidation.h index 5a6bc4ed..379753f5 100644 --- a/rugnux/ModelValidation.h +++ b/rugnux/ModelValidation.h @@ -3,6 +3,7 @@ #pragma once +#include #include #include @@ -33,8 +34,12 @@ 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 +// the model's space group (e.g. data P4(1)2(1)2, model P4(3)2(1)2 - indistinguishable from merged +// intensities), the observed reflections are reindexed into the model's hand before comparison. ModelValidationResult ValidateAgainstModel(const std::vector &merged, const UnitCell &cell, const std::string &model_path, const std::string &output_prefix, - Logger &logger); + Logger &logger, + std::optional data_space_group_number = std::nullopt); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 92fada20..89a54a97 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -872,8 +872,10 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { if (result.consensus_cell && write_files && !config_.model_path.empty()) { phase("Validating against model"); + const auto data_sg = experiment_.GetSpaceGroupNumber(); ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path, - config_.output_prefix, logger); + config_.output_prefix, logger, + data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt); } } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 02ac869d..ea166605 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -914,6 +914,19 @@ int main(int argc, char **argv) { : 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); + logger.Info("Fixing space group from reference MTZ: {} ({})", + reference.space_group_name, *space_group_number); + } + if (!fixed_reference_unit_cell.has_value() && reference.cell.has_value()) { + fixed_reference_unit_cell = reference.cell; + logger.Info("Fixing reference unit cell from reference MTZ (indexing may drift within tolerance)"); + } } catch (const std::exception &e) { logger.Error("Error reading reference MTZ {}: {}", ref_mtz, e.what()); exit(EXIT_FAILURE); @@ -1038,9 +1051,12 @@ int main(int argc, char **argv) { error_model_isa > 0 ? fmt::format("{:.2f}", error_model_isa) : "?", twinning, output_prefix); - if (!output_prefix.empty() && !model_pdb.empty()) + if (!output_prefix.empty() && !model_pdb.empty()) { + const auto data_sg = experiment.GetSpaceGroupNumber(); ValidateAgainstModel(merged_reflections, *experiment.GetUnitCell(), model_pdb, - output_prefix, logger); + output_prefix, logger, + data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt); + } return 0; }