diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 57d80ad21..e7f40b818 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -22,6 +22,7 @@ * `rugnux --model` fits the model's scale, anisotropic B and bulk-solvent parameters on the working reflections only, so the R-free it reports is measured against a model no free reflection helped scale. * The bulk-solvent parameters of `rugnux --model` are searched over their physically meaningful range instead of being fitted without bounds, so a model is never scaled with a solvent term that has silently switched itself off. * The rigid-body placement of `rugnux --model` uses the same bounded bulk solvent as the reported fit, so a model is no longer placed against a target carrying a solvent term with no physical meaning. +* `rugnux --model` puts the model into the data's own description of the lattice before placing it, so a model whose cell is written on other axes - I-centred where the run indexed C-centred, a different unique axis, a permuted orthorhombic cell - is placed rather than scored where it was read; `MODEL_CHANGE_OF_BASIS=` and `MODEL_SETTING_AS_READ=` report it when it happens. * The rugnux results report opens with a summary - `VERDICT=` (`OK`, `WARNINGS`, `UNUSABLE`, `FAILED`), `VERDICT_TEXT=`, `PATHOLOGY_FLAGS=` with one closed-vocabulary code per condition that warned, and the `WARNING:` lines, which used to close the file - and the sections after it are renumbered 1-5 with no gaps. * `rugnux --developer` writes the full results report - the pipeline-internal keys and the long explanations the default report now leaves out - and `--finalist-ledger` adds the evidence for every space group the search considered, not only the one it adopted. * The results report warns when the merged data carry no usable signal and when too little of reciprocal space was measured inside the fitted resolution, and omits `FITTED_RESOLUTION` where the CC1/2 curve it is fitted on never falls off. diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 669f83398..cc136e487 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -140,6 +140,184 @@ void reorient_about_centroid(gemmi::Model &model, const gemmi::Mat33 &rot) { SetModelPositions(model, pos); } +// Move a model into another cell by keeping its fractional coordinates. Correct only when the two +// cells describe the same axes in the same order, which is what the probe below is there to arrange. +void refractionalize_into(gemmi::Structure &st, const gemmi::UnitCell &target) { + const gemmi::UnitCell from = st.cell; + for (gemmi::Model &m : st.models) + for (gemmi::Chain &ch : m.chains) + for (gemmi::Residue &r : ch.residues) + for (gemmi::Atom &a : r.atoms) + a.pos = target.orthogonalize(from.fractionalize(a.pos)); + st.cell = target; +} + +// --- putting the model in the data's description of the lattice --------------------------------- +// +// A model arrives in the cell its depositor chose and rugnux indexes in the cell its own reduction +// chose, and the two are often 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. +// The space-group NUMBER is identical in every one of those, so no comparison of numbers can see it, +// and re-fractionalizing straight across such a pair scrambles the model. The rigid body further +// down cannot undo it either - six parameters about a centroid are not a change of basis - so the +// run would otherwise report a placement R-free near 0.6 for data that are perfectly good. + +// Do two cells describe the same box, axis for axis and in the same order? Deliberately loose: what +// this gate lets through is scored against the data afterwards, so being generous costs seconds, +// while being tight loses the case the probe exists for - a genuinely non-isomorphous pair can +// differ by a couple of degrees in beta and still be the same description. +bool cells_correspond(const gemmi::UnitCell &a, const gemmi::UnitCell &b) { + auto len = [](double x, double y) { return std::fabs(x - y) <= 0.05 * std::max(x, y); }; + auto ang = [](double x, double y) { return std::fabs(x - y) <= 3.0; }; + return len(a.a, b.a) && len(a.b, b.b) && len(a.c, b.c) + && ang(a.alpha, b.alpha) && ang(a.beta, b.beta) && ang(a.gamma, b.gamma); +} + +// Every integer change of basis with entries in {-1,0,1} and determinant +1 that carries `from` onto +// `to`. That set is small, and it is the set that matters here: an axis permutation, a sign flip, +// and the one off-diagonal that turns an I-centred monoclinic cell into a C-centred one all live in +// it. A halved or a doubled axis does not, and must not - those are indexing errors, not a choice of +// description, and a probe that silently repaired them would hide a real defect. +std::vector cell_mapping_operators(const gemmi::UnitCell &from, const gemmi::UnitCell &to) { + std::vector out; + for (int pattern = 0; pattern < 19683; pattern++) { // 3^9 matrices over {-1,0,1} + int e[9], v = pattern; + for (int &x : e) { x = v % 3 - 1; v /= 3; } + const int det = e[0] * (e[4] * e[8] - e[5] * e[7]) + - e[1] * (e[3] * e[8] - e[5] * e[6]) + + e[2] * (e[3] * e[7] - e[4] * e[6]); + if (det != 1) + continue; // +1 keeps both the volume and the hand + gemmi::Op op = gemmi::Op::identity(); + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + op.rot[i][j] = e[3 * i + j] * gemmi::Op::DEN; + gemmi::UnitCell moved = from; // changed_basis_forward is not const + if (cells_correspond(moved.changed_basis_forward(op, false), to)) + out.push_back(op); + } + return out; +} + +// Two candidates that differ by a rotation the model's own group already has describe the same +// structure, so only one of each class is worth scoring. On a holohedral cell that collapses two +// dozen candidates to one, which is what keeps the probe free on the ordinary isomorphous run. +bool same_frame(const gemmi::Op &a, const gemmi::Op &b, const gemmi::GroupOps &gops) { + const gemmi::Op::Rot d = a.inverse().combine(b).rot; + for (const gemmi::Op &s : gops.sym_ops) + if (s.rot == d) + return true; + return false; +} + +// 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 a group that is 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 is the everyday example - it needs (1/4, 1/4, 1/4) before it reads as P 21 21 21 again. +// So 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 order tries the common shifts first, +// so the search almost always ends on one of its first few candidates. +gemmi::Op with_origin_shift(const gemmi::SpaceGroup *sg, const gemmi::Op &op, + const gemmi::SpaceGroup **named) { + static const int TWELFTHS[] = {0, 6, 3, 9, 4, 8, 2, 10, 1, 5, 7, 11}; + for (int i : TWELFTHS) + for (int j : TWELFTHS) + for (int k : TWELFTHS) { + gemmi::Op shifted = op; + shifted.tran = {i * gemmi::Op::DEN / 12, j * gemmi::Op::DEN / 12, + k * gemmi::Op::DEN / 12}; + gemmi::GroupOps gops = sg->operations(); + gops.change_basis_forward(shifted); + if (const gemmi::SpaceGroup *found = gemmi::find_spacegroup_by_ops(gops)) { + *named = found; + return shifted; + } + } + *named = nullptr; + return op; +} + +// Put a model through a change of basis: coordinates, cell and space group together. Returns false - +// leaving the model untouched - when no origin shift makes the transformed group one GEMMI can name, +// which is how a basis that would leave a standard setting is refused rather than adopted. +bool change_model_basis(gemmi::Structure &st, const gemmi::SpaceGroup *&sg, gemmi::Op &op) { + const gemmi::SpaceGroup *moved = nullptr; + op = with_origin_shift(sg, op, &moved); + if (moved == nullptr) + return false; + gemmi::UnitCell old_cell = st.cell; + gemmi::Op rot_only = op; + rot_only.tran = {0, 0, 0}; // the cell follows the axes; only the atoms feel the origin shift + const gemmi::UnitCell new_cell = old_cell.changed_basis_forward(rot_only, false); + for (gemmi::Model &m : st.models) + for (gemmi::Chain &ch : m.chains) + for (gemmi::Residue &r : ch.residues) + for (gemmi::Atom &a : r.atoms) { + const gemmi::Fractional f = old_cell.fractionalize(a.pos); + const std::array t = op.apply_to_xyz({{f.x, f.y, f.z}}); + a.pos = new_cell.orthogonalize(gemmi::Fractional(t[0], t[1], t[2])); + } + st.cell = new_cell; + st.spacegroup_hm = moved->xhm(); + sg = moved; + return true; +} + +// The coarse shell the frame is decided on. A frame that is wrong is wrong at low resolution, so the +// probe never goes near the resolution the real fit uses - that is what makes trying every candidate +// affordable. The very lowest resolution is left out with it: there a bulk solvent this scorer does +// not model would dominate, equally for every candidate, and only add noise to the comparison. +constexpr double FRAME_PROBE_D_MIN = 3.5; +constexpr double FRAME_PROBE_D_MAX = 8.0; +constexpr size_t FRAME_PROBE_MIN_REFLECTIONS = 200; + +// R of the model against the observed amplitudes after an overall scale and an isotropic B, over +// that coarse shell. Only the ranking is ever used, never the value. +double frame_probe_r(const gemmi::Structure &st, const gemmi::SpaceGroup *sg, + const std::vector &obs, double d_min) { + const double probe_d_min = std::max(d_min, FRAME_PROBE_D_MIN); + gemmi::DensityCalculator dc; + dc.d_min = probe_d_min; + dc.rate = 1.5; + dc.set_grid_cell_and_spacegroup(st); + dc.set_refmac_compatible_blur(st.models[0]); + dc.put_model_density_on_grid(st.models[0]); + gemmi::AsuData> fcalc = + gemmi::transform_map_to_f_phi(dc.grid, true) + .prepare_asu_data(probe_d_min, dc.blur, false, false, false); + + gemmi::GroupOps gops = sg->operations(); + gemmi::ReciprocalAsu asu(sg); + gemmi::AsuData> fobs; + fobs.unit_cell_ = st.cell; + fobs.spacegroup_ = sg; + // The low-resolution cap is dropped, rather than the probe abandoned, when a small cell does not + // put enough reflections in the shell. + for (double d_max : {FRAME_PROBE_D_MAX, 1e9}) { + fobs.v.clear(); + for (const MergedReflection &r : obs) { + if (std::isnan(r.F) || r.d < probe_d_min || r.d > d_max) + continue; + gemmi::Miller h{{r.h, r.k, r.l}}; + if (!asu.is_in(h)) + h = asu.to_asu(h, gops).first; + fobs.v.push_back({h, {r.F, 1.0f}}); + } + if (fobs.v.size() >= FRAME_PROBE_MIN_REFLECTIONS) + break; + } + if (fobs.v.empty()) + return 1.0; + fobs.ensure_asu(); + fobs.ensure_sorted(); + gemmi::Scaling scaling(st.cell, sg); + scaling.use_solvent = false; + scaling.prepare_points(fcalc, fobs, nullptr); + scaling.fit_isotropic_b_approximately(); + return scaling.calculate_r_factor(); +} + + } // namespace ModelValidationResult ValidateAgainstModel(const std::vector &merged, @@ -227,17 +405,61 @@ ModelValidationResult ValidateAgainstModel(const std::vector & return result; } - // Re-fractionalize the model into the data cell (rigid cell adjustment; no refinement). const gemmi::UnitCell data_cell = cell; // UnitCell -> gemmi::UnitCell + + // --- same lattice, whose description of it? --- + // Scored, never asserted. The identity is always among the candidates, so a probe that finds + // nothing better than leaving the model where it is leaves it there; and the ordinary isomorphous + // run - where the only operators mapping the cell are the identity and its own symmetry + // equivalents - collapses to a single candidate and never reaches the scoring at all. if (data_cell.is_crystal()) { - gemmi::UnitCell old = st.cell; - for (gemmi::Model &m : st.models) - for (gemmi::Chain &ch : m.chains) - for (gemmi::Residue &r : ch.residues) - for (gemmi::Atom &a : r.atoms) - a.pos = data_cell.orthogonalize(old.fractionalize(a.pos)); - st.cell = data_cell; + const gemmi::GroupOps model_gops = sg->operations(); + std::vector frames{gemmi::Op::identity()}; + for (const gemmi::Op &op : cell_mapping_operators(st.cell, data_cell)) { + bool seen = false; + for (const gemmi::Op &kept : frames) + seen = seen || same_frame(kept, op, model_gops); + if (!seen) + frames.push_back(op); + } + if (frames.size() > 1) { + logger.Info("Model validation: the model's cell {:.2f} {:.2f} {:.2f} {:.1f} {:.1f} {:.1f} ({}) " + "is not how the data describe this lattice; scoring {} change(s) of basis", + st.cell.a, st.cell.b, st.cell.c, st.cell.alpha, st.cell.beta, st.cell.gamma, + sg->xhm(), frames.size() - 1); + size_t best = 0; + double best_r = 0; + for (size_t i = 0; i < frames.size(); i++) { + gemmi::Structure trial = st; + const gemmi::SpaceGroup *trial_sg = sg; + if (i > 0 && !change_model_basis(trial, trial_sg, frames[i])) + continue; // no origin shift names the transformed group; not a basis we can take + refractionalize_into(trial, data_cell); + trial.setup_cell_images(); + const double r = frame_probe_r(trial, trial_sg, obs, d_min); + logger.Info("Model validation: {:<12} -> {:<12} R {:.4f} on the coarse shell", + frames[i].triplet(), trial_sg->xhm(), r); + if (i == 0 || r < best_r) { + best_r = r; + best = i; + } + } + if (best > 0) { + result.setting_as_read = sg->xhm(); + if (change_model_basis(st, sg, frames[best])) { + result.change_of_basis_op = frames[best]; + logger.Info("Model validation: model put through {} into {} - the data's own " + "description of the same lattice", frames[best].triplet(), sg->xhm()); + } else { + result.setting_as_read.clear(); + } + } + } } + + // Re-fractionalize the model into the data cell (rigid cell adjustment; no refinement). + if (data_cell.is_crystal()) + refractionalize_into(st, data_cell); st.setup_cell_images(); const gemmi::UnitCell &ucell = st.cell; diff --git a/rugnux/ModelValidation.h b/rugnux/ModelValidation.h index 2dac0906a..21b87b706 100644 --- a/rugnux/ModelValidation.h +++ b/rugnux/ModelValidation.h @@ -110,6 +110,14 @@ struct ModelValidationResult { double r_work_sigma = 0.0; // (null mean - r_work) / null sd; positive = better than random int null_replicates = 0; + // The change of basis the model was put through to reach the data's description of the same + // lattice - an axis permutation, a sign flip, I-centred against C-centred. Identity where the two + // already agreed, which is the ordinary case. setting_as_read is the setting the model arrived + // in, kept only when it changed, so the report can say what was done rather than only that + // something was. + gemmi::Op change_of_basis_op = gemmi::Op::identity(); + std::string setting_as_read; + // The alternative-indexing operator picked by R-free, identity where none was needed or where the // choice did not beat its own null. AdoptModelFrame below applies it to the reflections that are // written out, so the file, the R-factors and the maps all describe one indexing. diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 1edcc30c4..7fc3f34d7 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -1225,6 +1225,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, {"NONE", "INDEXING", "ENANTIOMORPH", "ENANTIOMORPH+INDEXING"})); Add(s, KeyBool("MODEL_ENANTIOMORPH_ADOPTED", mv.adopted_model_enantiomorph)); Add(s, KeyText("MODEL_INDEXING_OPERATOR", mv.indexing_op.triplet())); + if (!(mv.change_of_basis_op == gemmi::Op::identity())) { + Add(s, KeyText("MODEL_CHANGE_OF_BASIS", mv.change_of_basis_op.triplet())); + Add(s, KeyText("MODEL_SETTING_AS_READ", mv.setting_as_read)); + } if (mv.indexing_probed) { Add(s, KeyReal("MODEL_INDEXING_MARGIN", mv.indexing_margin, "{:.4f}", true)); Add(s, KeyText("MODEL_INDEXING_MARGIN_NULL", diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index ed7282c2e..5d98a6318 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include @@ -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");