rugnux: fix SG/cell from a reference MTZ, and reindex data into the model's enantiomorph
Two ways the processing now adopts symmetry from an external reference: - Reference MTZ (-z) fixes the space group and unit cell, unless -S / -C override them (the explicit flag always wins - so -S with the "wrong" enantiomorph is allowed). The cell is a soft reference: indexing may still drift within tolerance, as with -C. - rugnux --model: when the data was merged 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), ValidateAgainstModel reindexes the observed reflections into the model's hand (via gemmi change_of_hand_op). This keeps the data setting consistent with the model; it does not change the R-factors, which use |F|. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
@@ -21,6 +21,7 @@
|
||||
#include <gemmi/mtz.hpp> // 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<MergedReflection> &
|
||||
const UnitCell &cell,
|
||||
const std::string &model_path,
|
||||
const std::string &output_prefix,
|
||||
Logger &logger) {
|
||||
Logger &logger,
|
||||
std::optional<int> data_space_group_number) {
|
||||
ModelValidationResult result;
|
||||
|
||||
// --- read the atomic model ---
|
||||
@@ -71,9 +73,32 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
|
||||
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<MergedReflection> reindexed;
|
||||
const std::vector<MergedReflection> *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<MergedReflection> &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<MergedReflection> &
|
||||
// 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;
|
||||
|
||||
@@ -3,6 +3,7 @@
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <optional>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
@@ -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<MergedReflection> &merged,
|
||||
const UnitCell &cell,
|
||||
const std::string &model_path,
|
||||
const std::string &output_prefix,
|
||||
Logger &logger);
|
||||
Logger &logger,
|
||||
std::optional<int> data_space_group_number = std::nullopt);
|
||||
|
||||
+3
-1
@@ -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<int>(static_cast<int>(*data_sg)) : std::nullopt);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
+18
-2
@@ -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<int64_t>(*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<int>(static_cast<int>(*data_sg)) : std::nullopt);
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
Reference in New Issue
Block a user