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