mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 73 Wiki Activity

model validation: write the model as it was placed

Browse Source 
`--model` re-fractionalizes the model into the data cell and then places it as
one rigid body, but the placed coordinates never reached disk. On a lysozyme
sweep against a non-isomorphous deposited model the move is 3.058 deg and
1.035 A, so a user overlaying their input model on rugnux's maps was out by
exactly that, and no file on disk corresponded to the maps at all.

`<prefix>_model.cif` is that file: the input's chains, residues, ligands,
waters, B-factors, occupancies and anisotropic Us, at the coordinates the maps
were computed from. `<prefix>.cif` is already the merged reflections, hence the
suffix.

The cell and space group come from the same two values WriteReflections is
given - the unit cell and DiffractionExperiment::GetSpaceGroupOrP1() after
AdoptModelFrame has settled the enantiomorph - so the coordinate file and the
.mtz beside it always agree. Taking them from the input model would not: with
data merged in P4(1)2(1)2 and a P4(3)2(1)2 model, the written reflections take
the model's group, which is neither the data's original label nor, when the
model is rejected, the model's own.

Written whenever the maps are, not only where the rigid-body step was
committed. The model is re-fractionalized and may be relabelled whatever the
placement decided, so an unmoved model is still not the input file; and a model
the null rejected is scored, placed and mapped like any other - the negative
result, and the case where the density is most worth looking at.

Nothing in the tree could write coordinates: gemmi_gph declared to_mmcif.hpp
but src/to_mmcif.cpp had been trimmed from the vendored subset. Both it and
to_pdb.cpp (to_mmcif.cpp calls its use_hetatm) are vendored from the same
gemmi 0.7.5 the rest of gemmi_gph comes from, unmodified, and every header they
include was already there. Same package, same MPL-2.0, same LICENSE.txt already
collected into licenses/gemmi.txt and already listed against `gemmi_gph/` in
THIRD_PARTY_NOTICES.md, so no new row and no new licence text.

Verified end to end: read back with gemmi the file differs from the
re-fractionalized input by exactly the reported 3.058 deg / 1.035 A with 0.0000 A
rms about that rigid move, and REFMAC5 at zero cycles against rugnux's own .mtz
starts at R-free 0.3667 where rugnux reports 0.3826 - against 0.5832 for the
unplaced input model, where rugnux reports 0.5911 before the placement.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
This commit is contained in:
leonarski_f
2026-09-02 14:47:38 +02:00
co-authored by Claude Opus 5
parent 8a0bbae3c4
commit 28a58e94b9
14 changed files with 2228 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/ModelValidationTest.cpp
+64
View File
@@ -15,6 +15,7 @@
#include "../rugnux/ModelValidation.h"
#include "../rugnux/RigidBodyRefine.h"
#include "../rugnux/SigmaA.h"
#include "../rugnux/WriteModel.h"
namespace {
// A synthetic P1 cell with two carbon atoms - enough for a reader to produce a Structure with
@@ -70,6 +71,16 @@ ATOM 2 C . CB GLY A 1 12.000 14.000 16.000 1.00 20.00
return pdb + "END\n";
}
// The same two atoms as kPdb, plus an anisotropic U, a partial occupancy and a water: a written
// model has to carry what the input carried, not only its coordinates.
const char *kPdbRich =
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n"
"ATOM 1 CA GLY A 1 10.000 12.000 14.000 1.00 20.00 C\n"
"ANISOU 1 CA GLY A 1 1000 1200 1400 100 200 300 C\n"
"ATOM 2 CB GLY A 1 12.000 14.000 16.000 0.60 25.00 C\n"
"HETATM 3 O HOH A 101 20.000 22.000 24.000 1.00 30.00 O\n"
"END\n";
std::string WriteTemp(const std::string &name, const char *content) {
std::ofstream f(name);
f << content;
@@ -211,3 +222,56 @@ TEST_CASE("ModelValidation_SigmaAWeightsFollowTheModelsAgreement", "[ModelValida
CHECK(unweighted.weight[0].m == 1.0);
CHECK(unweighted.weight[0].d == 1.0);
}
// The model rugnux scored has to reach disk, or a user overlaying their input model on rugnux's maps
// is wrong by the whole rigid-body shift. Two things have to hold: the file carries what the input
// carried, and it carries the cell and space group the reflection files beside it are written in -
// which, once --model has adopted the model's enantiomorph, is neither the data's original label nor
// necessarily the input model's.
TEST_CASE("WriteModel_KeepsTheContentAndTakesTheGivenFrame", "[ModelValidation]") {
Logger logger("WriteModel_KeepsTheContentAndTakesTheGivenFrame");
const auto input = WriteTemp("write_model_test_input.pdb", kPdbRich);
gemmi::Structure st = gemmi::read_structure_gz(input, gemmi::CoorFormat::Detect);
// 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.
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_test", logger);
const std::string written = "write_model_test_model.cif";
REQUIRE(std::filesystem::exists(written));
const gemmi::Structure back = gemmi::read_structure_gz(written, gemmi::CoorFormat::Detect);
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 == 96);
REQUIRE(back.models.size() == 1);
REQUIRE(back.models[0].chains.size() == 1);
const gemmi::Chain &chain = back.models[0].chains[0];
REQUIRE(chain.residues.size() == 2); // GLY and the water
REQUIRE(chain.residues[0].atoms.size() == 2);
REQUIRE(chain.residues[1].atoms.size() == 1);
// The coordinates are written as they stand - WritePlacedModel changes the label, not the atoms.
const gemmi::Atom &ca = chain.residues[0].atoms[0];
const gemmi::Atom &cb = chain.residues[0].atoms[1];
CHECK(ca.pos.x == Catch::Approx(10.0));
CHECK(ca.pos.z == Catch::Approx(14.0));
CHECK(cb.pos.y == Catch::Approx(14.0));
// B-factors, occupancies and the anisotropic U survive the round trip.
CHECK(ca.b_iso == Catch::Approx(20.0));
CHECK(cb.b_iso == Catch::Approx(25.0));
CHECK(cb.occ == Catch::Approx(0.60));
CHECK(ca.aniso.nonzero());
CHECK(ca.aniso.u11 == Catch::Approx(0.1000).margin(1e-4));
CHECK(ca.aniso.u12 == Catch::Approx(0.0100).margin(1e-4));
CHECK(chain.residues[1].name == "HOH");
std::filesystem::remove(input);
std::filesystem::remove(written);
}