From 7dc066c0f34a0a84103f6ce7441a51613cf11e2d Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 16 Sep 2026 11:52:24 +0200 Subject: [PATCH] model: the placed coordinates carry the group they obey, not the data's WritePlacedModel stamped the reflection files' space group on the placed model unconditionally. That is right whenever the coordinates obey that group - the ordinary case, and the subgroup case too, where the file just holds more than one asymmetric unit of it - but wrong when the data were merged in a supergroup of the model's group, or in the other hand. The extra operations are then a claim the contents contradict: a refinement program expands them, refines a structure that is not the model, and reports nothing. Measured on a public deposition whose data this run over-merged across a pseudo-symmetry twofold: refining the written pair gave R-free 0.495, where the same coordinates and the same reflections under the model's own group give 0.285. The label alone was worth 0.21 of R-free. The group is now taken from the argument only where every one of its operations is one the model already has, and the model keeps its own otherwise. The two files then disagree, which is the truth about them and is warned about, rather than agreeing on something false. Co-Authored-By: Claude Opus 5 (1M context) --- rugnux/WriteModel.cpp | 54 ++++++++++++++++++++++++----- rugnux/WriteModel.h | 7 ++-- tests/ModelValidationTest.cpp | 64 +++++++++++++++++++++++++++++++++-- 3 files changed, 112 insertions(+), 13 deletions(-) diff --git a/rugnux/WriteModel.cpp b/rugnux/WriteModel.cpp index 9b90faef7..1d0885f21 100644 --- a/rugnux/WriteModel.cpp +++ b/rugnux/WriteModel.cpp @@ -3,6 +3,7 @@ #include "WriteModel.h" +#include #include #include // make_mmcif_document @@ -12,6 +13,22 @@ #include "../common/Logger.h" +namespace { + +// Do the coordinates of a model built in `coordinates_group` obey `label` as well? They do when every +// operation of `label` is one of theirs, i.e. when `label` is a subgroup: the file then holds more +// than one asymmetric unit of it, which is allowed and is what a model always does in P1. They do not +// when `label` has an operation of its own, and then the label is a lie a reader acts on - it +// generates atoms from that operation, and the crystal it builds is not the one in the file. +bool coordinates_obey(const gemmi::SpaceGroup &coordinates_group, const gemmi::SpaceGroup &label) { + const std::vector coordinate_ops = coordinates_group.operations().all_ops_sorted(); + const std::vector label_ops = label.operations().all_ops_sorted(); + return std::includes(coordinate_ops.begin(), coordinate_ops.end(), + label_ops.begin(), label_ops.end()); +} + +} // namespace + void WritePlacedModel(const gemmi::Structure &placed, const UnitCell &cell, const gemmi::SpaceGroup &space_group, @@ -25,13 +42,23 @@ void WritePlacedModel(const gemmi::Structure &placed, const std::string pdb_path = output_prefix + "_model.pdb"; gemmi::Structure st = placed; - // Not the model's own cell and group but the data's, taken from the same two values the reflection - // files are written from. The coordinates already sit in this cell; what is set here is the label - // the file carries, which must match the .mtz beside it - the enantiomorph in particular, since - // --model can adopt the model's and that is neither the data's original label nor, necessarily, - // the one the input model arrived with. + // The cell is the data's, taken from the same value the reflection files are written from: the + // coordinates were re-fractionalized into it and genuinely sit there. + // + // The group is the data's too wherever the coordinates obey it, which is the ordinary case and the + // one that matters - the label then matches the .mtz beside it, the enantiomorph in particular, + // since --model can adopt the model's and that is neither the data's original label nor, + // necessarily, the one the input model arrived with. Where they do not obey it the model's own + // group is written instead. Data merged in a supergroup of the model's - a run that over-merged + // across a pseudo-symmetry twofold, or a hand the model's fit did not earn - would otherwise put + // operations on the file that its contents contradict, and a refinement program expands them: the + // structure it refines is not the model, and it says so nowhere. A pair whose two files disagree + // is worth a warning; a coordinate file that is quietly wrong on its own is not worth writing. + const gemmi::SpaceGroup *model_group = st.find_spacegroup(); + const bool takes_model_group = model_group != nullptr && !coordinates_obey(*model_group, space_group); + const gemmi::SpaceGroup &label = takes_model_group ? *model_group : space_group; st.cell = cell; - st.spacegroup_hm = space_group.xhm(); + st.spacegroup_hm = label.xhm(); st.setup_cell_images(); // Fills in entity types and label_asym_id for a model read from a PDB, which carries neither. // Both are no-ops where the input already had them, i.e. for an mmCIF input. @@ -60,7 +87,16 @@ void WritePlacedModel(const gemmi::Structure &placed, logger.Warning("Model validation: the placed model could not also be written as PDB ({}); " "{} has it", e.what(), cif_path); } - logger.Info("Model validation: the model as placed against these data written to {}{} " - "(cell and space group {} as in the reflection files)", - cif_path, pdb_written ? " and " + pdb_path : "", space_group.short_name()); + if (!takes_model_group) { + logger.Info("Model validation: the model as placed against these data written to {}{} " + "(cell and space group {} as in the reflection files)", + cif_path, pdb_written ? " and " + pdb_path : "", space_group.short_name()); + } else { + logger.Warning("Model validation: the model as placed written to {}{} in ITS OWN space group " + "{}, not the {} the reflection files carry - its coordinates do not obey {}, so " + "the two files describe the same crystal in different symmetry and cannot be " + "refined against each other as they stand", + cif_path, pdb_written ? " and " + pdb_path : "", + label.short_name(), space_group.short_name(), space_group.short_name()); + } } diff --git a/rugnux/WriteModel.h b/rugnux/WriteModel.h index 356a47952..1c5bf58d1 100644 --- a/rugnux/WriteModel.h +++ b/rugnux/WriteModel.h @@ -18,8 +18,11 @@ class Logger; // // `cell` and `space_group` must be the ones the reflection files beside it are written in, i.e. the // unit cell and DiffractionExperiment::GetSpaceGroupOrP1() that WriteReflections is given after -// AdoptModelFrame has settled the enantiomorph. A coordinate file in a different frame from the .mtz -// next to it is worse than no coordinate file at all. +// AdoptModelFrame has settled the enantiomorph. The cell is always taken; the group is taken wherever +// the model's coordinates obey it, which is the ordinary case. Where they do not - data merged in a +// supergroup of the model's group, or in the other hand - the model keeps its own group and a warning +// says the two files disagree, because a label whose operations the coordinates contradict makes a +// refinement program build a different structure and say nothing. void WritePlacedModel(const gemmi::Structure &placed, const UnitCell &cell, const gemmi::SpaceGroup &space_group, diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index 5d98a6318..b7e4ec822 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -398,9 +398,12 @@ TEST_CASE("WriteModel_KeepsTheContentAndTakesTheGivenFrame", "[ModelValidation]" const auto input = WriteTemp("write_model_test_input.pdb", kPdbRich); gemmi::Structure st = gemmi::read_structure_gz(input, gemmi::CoorFormat::Detect); + // The model in the data's own group, which is the ordinary case: the frame it is written in is + // then wholly the caller's. The cell it arrives in (40/50/60) is not the data's and is replaced. + st.spacegroup_hm = "P 43 21 2"; - // A frame that is neither the model's (P 1, 40/50/60) nor anything derived from it: the tetragonal - // lysozyme cell and one enantiomorph of its group, standing in for what AdoptModelFrame settled. + // The tetragonal lysozyme cell and one enantiomorph of its group, standing in for what + // AdoptModelFrame settled. const UnitCell data_cell{.a = 79, .b = 79, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}; const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 43 21 2"); REQUIRE(sg != nullptr); @@ -441,6 +444,63 @@ TEST_CASE("WriteModel_KeepsTheContentAndTakesTheGivenFrame", "[ModelValidation]" std::filesystem::remove(written); } +// A model cannot be labelled with a group its own coordinates do not obey. Where the data were merged +// in a supergroup of the model's - an over-merge across a pseudo-symmetry operation, or the hand the +// model's fit did not earn - that operation would generate atoms the model does not contain, and a +// refinement program acts on it silently. The model keeps its own group and the caller is warned. +TEST_CASE("WriteModel_DoesNotStampAGroupTheCoordinatesDoNotObey", "[ModelValidation]") { + Logger logger("WriteModel_DoesNotStampAGroupTheCoordinatesDoNotObey"); + + const auto input = WriteTemp("write_model_sg_test_input.pdb", kPdbRich); // P 1 + const gemmi::Structure st = gemmi::read_structure_gz(input, gemmi::CoorFormat::Detect); + + const UnitCell data_cell{.a = 79, .b = 79, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}; + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 43 21 2"); + REQUIRE(sg != nullptr); + + WritePlacedModel(st, data_cell, *sg, "write_model_sg_test", logger); + const std::string written = "write_model_sg_test_model.cif"; + REQUIRE(std::filesystem::exists(written)); + + const gemmi::Structure back = gemmi::read_structure_gz(written, gemmi::CoorFormat::Detect); + // The cell is still the data's - the coordinates do sit in it - but the group is the model's own. + CHECK(back.cell.a == Catch::Approx(79.0)); + CHECK(back.cell.c == Catch::Approx(38.0)); + REQUIRE(back.find_spacegroup() != nullptr); + CHECK(back.find_spacegroup()->number == 1); + + std::filesystem::remove(input); + std::filesystem::remove(written); + std::filesystem::remove("write_model_sg_test_model.pdb"); +} + +// The other direction is not a mismatch at all: a model whose coordinates obey more symmetry than the +// data were merged in obeys the data's group too - the file simply holds more than one asymmetric unit +// of it, which is what a model in P1 always does. The reflection files' group is kept, as before. +TEST_CASE("WriteModel_KeepsTheDataGroupWhenItIsASubgroup", "[ModelValidation]") { + Logger logger("WriteModel_KeepsTheDataGroupWhenItIsASubgroup"); + + const auto input = WriteTemp("write_model_sub_test_input.pdb", kPdbRich); + gemmi::Structure st = gemmi::read_structure_gz(input, gemmi::CoorFormat::Detect); + st.spacegroup_hm = "P 43 21 2"; + + const UnitCell data_cell{.a = 79, .b = 79, .c = 38, .alpha = 90, .beta = 90, .gamma = 90}; + const gemmi::SpaceGroup *p1 = gemmi::find_spacegroup_by_name("P 1"); + REQUIRE(p1 != nullptr); + + WritePlacedModel(st, data_cell, *p1, "write_model_sub_test", logger); + const std::string written = "write_model_sub_test_model.cif"; + REQUIRE(std::filesystem::exists(written)); + + const gemmi::Structure back = gemmi::read_structure_gz(written, gemmi::CoorFormat::Detect); + REQUIRE(back.find_spacegroup() != nullptr); + CHECK(back.find_spacegroup()->number == 1); + + std::filesystem::remove(input); + std::filesystem::remove(written); + std::filesystem::remove("write_model_sub_test_model.pdb"); +} + // CC(model, data) has to follow where the signal actually is, or it cannot support the one-sided // claim it exists for. The check is closed: the "observed" intensities are the model's own with // Gaussian noise added, and how much noise is chosen per shell - almost none in the first, some in