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

model validation: the model is put into the data's description of the lattice before it is placed

Browse Source 
A deposited model arrives in the cell its depositor chose and rugnux indexes in the cell its
own reduction chose. The two are frequently different descriptions of the SAME lattice -
I-centred where the other is C-centred, unique axis c where the other took b, a cyclic
permutation of an orthorhombic cell - and the space-group NUMBER is identical in every one
of those, so the existing enantiomorph check, which compares numbers, cannot see any of them.
Re-fractionalizing straight across such a pair scrambles the model, and the rigid body cannot
undo it: six parameters about a centroid are not a change of basis. The run then reported a
placement R-free near 0.6, CC(model,data) around 0.3 and MODEL_FIT= NOT_TESTED, with no
warning, for data that were perfectly good.

Before re-fractionalizing, every integer change of basis with entries in {-1,0,1} and
determinant +1 that carries the model cell onto the data cell is enumerated, quotiented by the
model's own point-group rotations, and scored on an 8-3.5 A shell after an overall scale and an
isotropic B. The lowest wins. The identity is always among the candidates, so a probe that
finds nothing better than leaving the model alone leaves it alone; and on the ordinary
isomorphous run the only operators that map the cell are the identity and its own symmetry
equivalents, so they collapse to one candidate and the scoring is never reached.

A change of basis is a matrix AND an origin shift, and the shift is not optional. An odd
permutation of a screw-axis group lands on the same group on a moved origin, whose operator
list GEMMI cannot name because it compares those lists exactly - swapping b and c in P 21 21 21
needs (1/4, 1/4, 1/4) first. Where the bare matrix names nothing, the shift that makes it name
something is searched for on the twelfths every crystallographic origin shift lies on.

The coarse shell is what makes trying every candidate affordable: a frame that is wrong is
wrong at low resolution. The lowest resolution is excluded with it, because there an unmodelled
bulk solvent would dominate - equally for every candidate, so it would only add noise. Measured
cost 0.04-0.52 s, and only on runs whose cells actually differ.

A halved or a doubled axis is deliberately outside the search space: those are indexing errors,
not a choice of description, and a probe that silently repaired one would hide a real defect.
So is the case where the lattice is right but the symmetry was undercalled to P1, which would
need the model expanded to P1 and a centred-to-primitive basis change.

Measured over the deposited-data battery: eight runs whose model had never been placed now
place it, R-free falling by 0.14 to 0.36 and CC(model,data) rising from 0.14-0.47 to 0.51-0.96;
density at atom centres 0.35-2.1 sigma to 1.5-4.8, map mean FOM 0.13-0.28 to 0.36-0.93. Five
controls chosen to exercise the probe moved by no more than run-to-run noise and took no
operator, one of them after correctly rejecting a real alternative frame at coarse R 0.56
against 0.19 for the identity.

The reflections are not touched: they keep the indexing the reduction determined, and it is the
model that is expressed in it. What this fixes is that the written model and the written
reflections now agree - Fcalc from the placed model against the merged amplitudes, both as
written, gives R 0.18-0.31 where it gave 0.54-0.65. The CRYST1 header was already correct
before, since WriteModel forces the data's cell and group, so a downstream program would have
loaded the pair without complaint and refined against scrambled atoms.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XbMaLuuwX4D58QCxZ686Ku
This commit is contained in:
leonarski_f
2026-09-08 11:06:07 +02:00
co-authored by Claude Opus 5
parent a0f000f5e9
commit 94ed864dbe
5 changed files with 312 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/ModelValidationTest.cpp
+69
View File
@@ -8,6 +8,7 @@
#include <filesystem>
#include <fstream>
#include <random>
#include <sstream>
#include <gemmi/mmread_gz.hpp>
@@ -324,6 +325,74 @@ TEST_CASE("ModelValidation_SigmaAWeightsFollowTheModelsAgreement", "[ModelValida
// 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.
// A depositor's model and rugnux's own reduction routinely describe the SAME lattice on different
// axes - I-centred against C-centred, unique axis b against c, a cyclic permutation of an
// orthorhombic cell - and the space-group NUMBER is identical in every one of those, so nothing that
// compares numbers can notice. Re-fractionalizing straight across such a pair scrambles the model,
// and the rigid body cannot undo it: six parameters about a centroid are not a change of basis. The
// check is closed - the data are the model's own structure factors - so a model that reaches the
// data's description of the lattice comes out at R near zero, and one that does not comes out near
// the 0.6 an unrelated structure gives.
TEST_CASE("ModelValidation_FindsTheDatasDescriptionOfTheLattice", "[ModelValidation]") {
Logger logger("ModelValidation_FindsTheDatasDescriptionOfTheLattice");
// The data: the cluster's own |F|^2 in its own frame, 30 x 34 x 38, P 21 21 21.
const auto data_model = WriteTemp("cob_data_test.pdb", ClusterPdb().c_str());
auto obs = ModelReferenceIntensities(data_model, {}, {}, 2.5, logger);
REQUIRE(obs.size() > 500);
for (size_t i = 0; i < obs.size(); i++) {
obs[i].F = std::sqrt(std::max(0.0f, obs[i].I));
obs[i].rfree_flag = (i % 20) == 0;
}
const UnitCell data_cell{.a = 30, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90};
const gemmi::SpaceGroup *data_sg = gemmi::find_spacegroup_by_name("P 21 21 21");
const std::string prefix = (std::filesystem::temp_directory_path() / "cob_test").string();
// The same structure as a deposition on permuted axes, a' = b, b' = c, c' = a. The cell is
// orthogonal, so the Cartesian coordinates permute with the axes.
std::string permuted = "CRYST1 34.000 38.000 30.000 90.00 90.00 90.00 P 21 21 21 4\n";
{
std::istringstream in(ClusterPdb());
std::string line;
char buf[96];
while (std::getline(in, line)) {
if (line.compare(0, 4, "ATOM") != 0)
continue;
std::snprintf(buf, sizeof buf, "%s%8.3f%8.3f%8.3f%s",
line.substr(0, 30).c_str(), std::stod(line.substr(38, 8)),
std::stod(line.substr(46, 8)), std::stod(line.substr(30, 8)),
line.substr(54).c_str());
permuted += buf;
permuted += "\n";
}
permuted += "END\n";
}
const auto other_setting = WriteTemp("cob_model_test.pdb", permuted.c_str());
const auto moved = ValidateAgainstModel(obs, data_cell, other_setting, prefix, logger, data_sg,
/*probe_indexing_ambiguity=*/false, 1, 1.0);
REQUIRE(moved.ok);
logger.Info("Change of basis test: took {} out of {}, R-work {:.4f}",
moved.change_of_basis_op.triplet(), moved.setting_as_read, moved.r_work);
CHECK_FALSE(moved.change_of_basis_op == gemmi::Op::identity());
CHECK(moved.r_work < 0.15);
// The isomorphous case, which is what nearly every run is: the model already describes the
// lattice the way the data do, the only operators that map its cell are its own symmetry, and
// nothing is scored or moved.
const auto same = ValidateAgainstModel(obs, data_cell, data_model, prefix, logger, data_sg,
false, 1, 1.0);
REQUIRE(same.ok);
CHECK(same.change_of_basis_op == gemmi::Op::identity());
CHECK(same.setting_as_read.empty());
CHECK(same.r_work < 0.15);
for (const char *suffix : {"_2fofc.ccp4", "_fofc.ccp4", "_anom.ccp4", "_maps.mtz"})
std::filesystem::remove(prefix + suffix);
std::filesystem::remove(data_model);
std::filesystem::remove(other_setting);
}
TEST_CASE("WriteModel_KeepsTheContentAndTakesTheGivenFrame", "[ModelValidation]") {
Logger logger("WriteModel_KeepsTheContentAndTakesTheGivenFrame");