diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 4b2c89dac..f24ccc0b3 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -201,8 +201,10 @@ There is no threshold on R behind it. What a model that explains nothing reaches set depends on its atom count and B-factors as much as on the data, so the same model is refitted — and re-placed as a rigid body, exactly as the real one is — from `MODEL_FIT_NULL_REPLICATES` random orientations about its own centroid, and `MODEL_FIT_SIGMA` is how far the real fit sits above that -distribution. The statistic is R-work, not R-free: nothing is refined against the working set here, so -R-work carries no optimism, and it is decided on an order of magnitude more reflections. +distribution. The statistic is R-work, not R-free - not because nothing is refined against the working +set (the placement's six parameters are), but because every null replicate is placed the same way, so +what they buy is bought on both sides and cancels; and it is decided on an order of magnitude more +reflections than R-free. That null is **only built where the model claims one of the two things it could change** — the enantiomorph, or an indexing other than the one the data were merged in. A model already in the data's diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index d1bc187a4..8548eaf3a 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -88,8 +88,9 @@ constexpr unsigned NULL_SEED = 20260902; // How far above its own null a fit has to sit before the model is allowed to decide anything. The cut // is in sigma of that null and not in R: measured, the R a model that explains nothing reaches moves // with the model's atom count and B-factors as much as with the data, so no value of R separates the -// two on its own. A model that belongs to the crystal measures around +4 sigma here and an unrelated -// one under +1. +// two on its own. Measured at nine replicates, the crystal's own model reads +17.7 sigma and an +// unrelated protein +1.7, so the cut sits in a gap an order of magnitude wider than the sd it is +// measured in. constexpr double MODEL_FIT_SIGMA = 3.0; // Mean and sample standard deviation of a small sample. @@ -366,7 +367,13 @@ ModelValidationResult ValidateAgainstModel(const std::vector & const std::vector before = ModelPositions(ms.st.models[0]); const RigidBodyRefineResult rb = RefineRigidBody(ms.st.models[0], ucell, *sg, fobs_work, d_min, logger); - if (rb.converged) { + if (!rb.converged) { + // RefineRigidBody leaves the model wherever the solver left it, usable answer or not, so + // the restore cannot be conditional on the same flag the re-fit is. + SetModelPositions(ms.st.models[0], before); + return out; + } + { compute_model_factors(ms); Fit moved = fit_model(ms, obs_in); const bool commit = moved.r_free < out.fit.r_free; @@ -451,6 +458,12 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // computed and reported exactly as in every other case. const bool decision_pending = result.model_enantiomorph_candidate || !(indexing.op == gemmi::Op::identity()); + // A null that cannot be built is a question that could not be put, and this design already has a + // state for that. Without this the run dies here - model validation runs BEFORE the reflection + // files are written, so one failed replicate would take the .mtz, .cif, .hkl and _unmerged.mtz + // with it, after the merge has already been paid for. The replicates are rotated models fed to a + // scaling path that throws on data it cannot pair up, so this is the adversarial input for it. + try { if (decision_pending) { // --- the null: this model, this data, in random orientations --- // An R-factor on its own says nothing about whether a model belongs to a crystal. What a model @@ -461,9 +474,10 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // reoriented about its own centroid and run through the identical fit and rigid-body placement, // and the real fit is asked how far above the resulting distribution it sits. // - // R-work, not R-free: nothing is refined against the working set here - the scale has four - // parameters and the placement six - so R-work carries no optimism, and it is decided on an order - // of magnitude more reflections than R-free. + // R-work, not R-free. Not because nothing is refined against it - the placement's six parameters + // are, and the scale's four - but because every null replicate is placed and scaled the same + // way, so whatever that optimism is worth is bought on both sides and cancels in + // (mean - real)/sd. R-work is then decided on an order of magnitude more reflections. // // The replicates are fitted to the data as merged even where a reindexing won above: a reindexing // is a relabelling of the same intensities, and a model in a random orientation has no more to do @@ -501,7 +515,13 @@ ModelValidationResult ValidateAgainstModel(const std::vector & result.null_r_work_sd = null_sd; auto record_fit = [&](const Placement &p) { - result.r_work_sigma = null_sd > 0 ? (null_mean - p.fit.r_work) / null_sd : 0.0; + // Guarded at a floor, not at zero: sd == 0 and sd == 1e-7 are one ulp apart and land on + // opposite verdicts, and a null with no spread - a model too small or too symmetric for a + // rotation about its own centroid to move |Fcalc| - is the case that produces it. Measured + // nulls sit near 0.60 with sd ~0.009, so this is nowhere near anything genuine. + constexpr double NULL_SD_FLOOR = 1e-3; + result.r_work_sigma = + null_sd >= NULL_SD_FLOOR ? (null_mean - p.fit.r_work) / null_sd : 0.0; result.model_fits = result.r_work_sigma >= MODEL_FIT_SIGMA; }; record_fit(real); @@ -518,8 +538,9 @@ ModelValidationResult ValidateAgainstModel(const std::vector & const auto [margin_mean, margin_sd] = mean_sd(null_margin); result.indexing_margin_null_mean = margin_mean; result.indexing_margin_null_sd = margin_sd; + constexpr double MARGIN_SD_FLOOR = 1e-4; result.indexing_margin_sigma = - margin_sd > 0 ? (indexing.margin - margin_mean) / margin_sd : 0.0; + margin_sd >= MARGIN_SD_FLOOR ? (indexing.margin - margin_mean) / margin_sd : 0.0; // The margin is judged against the margin a random placement produces, not against a value. // A random model also picks a winner, and on these data it picks one by a comparable lead // (measured), so the raw margin says nothing on its own. @@ -553,6 +574,23 @@ ModelValidationResult ValidateAgainstModel(const std::vector & "reflections and no null was run; R-work {:.4f}, R-free {:.4f}", real.fit.r_work, real.fit.r_free); } + } catch (const std::exception &e) { + // The null could not be built, so no decision is licensed and none is taken - which is the + // NOT_TESTED state, not a failure of the run. If the indexing probe had already won on a + // relabelling, the fit that survives describes those reflections and the files will not carry + // them, so it is remade on the data as merged. + logger.Warning("Model validation: the null could not be built ({}), so the model decides " + "nothing; the reflections keep the group and the indexing they were merged in", + e.what()); + result.fit_tested = false; + result.model_fits = false; + result.indexing_probed = false; + result.indexing_op = gemmi::Op::identity(); + result.r_work_sigma = 0.0; + SetModelPositions(st.models[0], as_read); + compute_model_factors(mdl); + real = place_and_fit(mdl, obs); + } Fit &best = real.fit; result.rigid_body_applied = real.rb_applied; @@ -739,7 +777,10 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // of this file works in the wrong one if the label disagrees with the .mtz beside it. mtz.spacegroup = result.adopted_model_enantiomorph ? sg - : (data_space_group != nullptr ? data_space_group : sg); + // No data group means the caller merged in P1, and P1 is what the + // reflection files beside this one carry - not the model's group. + : (data_space_group != nullptr ? data_space_group + : gemmi::find_spacegroup_by_number(1)); mtz.set_cell_for_all(ucell); mtz.add_dataset("model_validation"); mtz.datasets.back().wavelength = wavelength_A; @@ -756,7 +797,10 @@ ModelValidationResult ValidateAgainstModel(const std::vector & mtz.add_column("FREE", 'I', -1, -1, false); std::vector data; for (size_t i = 0; i < terms.size(); i++) { - const auto phi_deg = static_cast(terms[i].phi * 180.0 / PI); + // [0, 360), the convention the MTZ 'P' columns carried before sigma_A weighting: std::arg + // returns (-180, 180] and a P column is not supposed to. + const double phi_wrapped = terms[i].phi * 180.0 / PI; + const auto phi_deg = static_cast(phi_wrapped < 0 ? phi_wrapped + 360.0 : phi_wrapped); data.insert(data.end(), {static_cast(terms[i].hkl[0]), static_cast(terms[i].hkl[1]), static_cast(terms[i].hkl[2]), static_cast(terms[i].fo), static_cast(terms[i].fc), diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 6d7784ae7..16d4dc502 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -663,8 +663,9 @@ std::string RenderResultReport(const std::string &output_prefix, << " depends on its atom count and B-factors as much as on the data - so the same model was\n" << " refitted, and re-placed, from " << mv.null_replicates << " random orientations about its own centroid, and\n" << " MODEL_FIT_SIGMA is how far the real fit sits above that null. R-work carries the\n" - << " decision because nothing was refined against it and it has far more reflections than\n" - << " R-free. MODEL_DECISIONS_TAKEN names what the model was allowed to change about the\n" + << " decision because the null was placed the same way the real fit was, so what those six\n" + << " placement parameters buy is bought on both sides and cancels - and it has far more\n" + << " reflections than R-free. MODEL_DECISIONS_TAKEN names what the model was allowed to change about the\n" << " written reflections; R-factors, maps and the rigid-body placement are reported either\n" << " way, because they describe the model, not the data.\n"; else diff --git a/rugnux/RigidBodyRefine.cpp b/rugnux/RigidBodyRefine.cpp index 29f4fcc99..a2a725c55 100644 --- a/rugnux/RigidBodyRefine.cpp +++ b/rugnux/RigidBodyRefine.cpp @@ -43,8 +43,8 @@ constexpr double JACOBIAN_STEP_FRACTION = 0.01; // Parameters are carried as six lengths in angstroms - the first three are the angle-axis rotation // vector multiplied by the model's rms radius, so a unit of each of the six moves a typical atom by -// the same amount. That makes the Jacobian step and the trust region isotropic in something -// physical, rather than mixing radians with angstroms. +// the same amount. That makes the Jacobian step isotropic in something physical, rather than mixing +// radians with angstroms. struct Placement { gemmi::Position centre; // the model centroid: rotating about it decorrelates R from t double rms_radius = 1.0; // rms distance of the atoms from the centroid @@ -185,6 +185,28 @@ private: Evaluator &ev_; }; +// The directions in which this space group's origin is free. Translating the whole cell content along +// one of them multiplies every F by a phase and leaves every |F| EXACTLY unchanged, so the target +// cannot determine that component: all three directions in P1, the unique axis in a polar group. The +// R-free gate cannot stand in for this - it is a function of |F| too, so along such a direction it +// sees only grid noise and commits or not by coin flip, while the other five parameters carry the +// noise in with them. The free directions are the common fixed subspace of the group's rotation +// parts, and the projector onto it is simply their average. +gemmi::Mat33 GaugeProjector(const gemmi::SpaceGroup &sg, const gemmi::UnitCell &cell) { + const gemmi::GroupOps gops = sg.operations(); + double m[3][3] = {}; + const double n = static_cast(gops.sym_ops.size()) * gemmi::Op::DEN; + for (const gemmi::Op &op : gops.sym_ops) + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + m[i][j] += static_cast(op.rot[i][j]) / n; + const gemmi::Mat33 mean(m[0][0], m[0][1], m[0][2], + m[1][0], m[1][1], m[1][2], + m[2][0], m[2][1], m[2][2]); + // Fractional projector taken into orthogonal space, where the parameters live. + return cell.orth.mat.multiply(mean).multiply(cell.frac.mat); +} + } // namespace std::vector ModelPositions(const gemmi::Model &model) { @@ -209,11 +231,12 @@ void SetModelPositions(gemmi::Model &model, const std::vector & // bound ligand its own six parameters would refine against evidence this data does not separately // carry, and the ligand is what the difference map is meant to show rather than model away. // -// Some of the translation can be a gauge rather than a quantity - the origin is free in all three +// Some of the translation would be a gauge rather than a quantity - the origin is free in all three // directions in P1 and along the unique axis in a polar group, and |F| does not change when the whole -// content moves along it. Nothing is done about that beyond the Levenberg-Marquardt damping, which -// leaves an undetermined direction where it started, and the caller's R-free gate, which throws away -// a step that wandered along one. +// content moves along it - so that component is projected out after every zone. Neither of the two +// things that might look like they cover it actually does: the R-free gate is a function of |F| and +// therefore blind to exactly this, and the LM damping follows the gauge column of the Jacobian, which +// is not zero but noise divided by the difference step. RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, @@ -238,16 +261,18 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, if (!(placement.rms_radius > 0)) return result; + std::vector ladder; for (double zone : LADDER) if (zone >= d_min) - result.zones.push_back(zone); - if (result.zones.empty()) - result.zones.push_back(d_min); + ladder.push_back(zone); + if (ladder.empty()) + ladder.push_back(d_min); + const gemmi::Mat33 gauge = GaugeProjector(sg, cell); Evaluator ev(model, cell, sg, base, placement); double q[6] = {0, 0, 0, 0, 0, 0}; bool any_zone_solved = false; - for (double zone : result.zones) { + for (double zone : ladder) { gemmi::AsuData> zone_obs; zone_obs.unit_cell_ = fobs.unit_cell_; zone_obs.spacegroup_ = fobs.spacegroup_; @@ -256,6 +281,7 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, zone_obs.v.push_back(hv); if (zone_obs.v.size() < 50) continue; + result.zones.push_back(zone); // the ladder WALKED, which a thin zone drops out of ev.SetZone(zone_obs, zone); ceres::Problem problem; @@ -271,9 +297,12 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, const auto zone_t0 = std::chrono::steady_clock::now(); ceres::Solve(options, &problem, &summary); any_zone_solved = any_zone_solved || summary.IsSolutionUsable(); + const gemmi::Vec3 along = gauge.multiply(gemmi::Vec3(q[3], q[4], q[5])); + q[3] -= along.x; q[4] -= along.y; q[5] -= along.z; logger.Debug("Rigid body zone {:.1f} A: {} reflections, {} iterations, {} evaluations, {:.2f} s, " "rotation {:.3f} deg, translation {:.3f} A", zone, zone_obs.v.size(), - summary.iterations.size() - 1, ev.evaluations - evaluations_before, + summary.iterations.empty() ? 0 : summary.iterations.size() - 1, + ev.evaluations - evaluations_before, std::chrono::duration(std::chrono::steady_clock::now() - zone_t0).count(), std::sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]) / placement.rms_radius * 180.0 / PI, std::sqrt(q[3]*q[3] + q[4]*q[4] + q[5]*q[5])); @@ -287,9 +316,13 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, result.shift_A = std::sqrt(q[3] * q[3] + q[4] * q[4] + q[5] * q[5]); result.seconds = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); - logger.Info("Model validation: rigid body over {} resolution zone(s) down to {:.1f} A, " - "{} evaluations in {:.2f} s: rotation {:.3f} deg, translation {:.3f} A", - result.zones.size(), result.zones.back(), result.evaluations, result.seconds, - result.angle_deg, result.shift_A); + if (!result.zones.empty()) + logger.Info("Model validation: rigid body over {} resolution zone(s) down to {:.1f} A, " + "{} evaluations in {:.2f} s: rotation {:.3f} deg, translation {:.3f} A", + result.zones.size(), result.zones.back(), result.evaluations, result.seconds, + result.angle_deg, result.shift_A); + else + logger.Info("Model validation: rigid body had no resolution zone with enough reflections to " + "run in; the model is left where it arrived"); return result; } diff --git a/rugnux/WriteModel.cpp b/rugnux/WriteModel.cpp index b9e36afbf..e3719bb8b 100644 --- a/rugnux/WriteModel.cpp +++ b/rugnux/WriteModel.cpp @@ -29,12 +29,19 @@ void WritePlacedModel(const gemmi::Structure &placed, 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. - gemmi::setup_entities(st); - - std::ofstream os(path); - gemmi::cif::write_cif_to_stream(os, gemmi::make_mmcif_document(st)); - if (!os) { - logger.Error("Model validation: cannot write the placed model to {}", path); + // setup_entities and make_mmcif_document can both throw, and this runs BEFORE the reflection files + // are written: a convenience deliverable must not be able to take the run's actual output with it. + // The map writer does the same. + try { + gemmi::setup_entities(st); + std::ofstream os(path); + gemmi::cif::write_cif_to_stream(os, gemmi::make_mmcif_document(st)); + if (!os) { + logger.Error("Model validation: cannot write the placed model to {}", path); + return; + } + } catch (const std::exception &e) { + logger.Warning("Model validation: could not write the placed model to {} ({})", path, e.what()); return; } logger.Info("Model validation: the model as placed against these data written to {} "