From 347e76fd8708d123ff2972c2ac977f690d151d3f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 27 Aug 2026 14:56:04 +0200 Subject: [PATCH] rugnux: read the 2Fo-Fc map at the atom centres cubically too Grid::interpolate_value defaults to trilinear, and on the d_min/3 map grid that under-reads a peak by around a tenth. The anomalous site list already asks for cubic; MEAN_ATOM_DENSITY_SIGMA was left on the default so a number already in existing reports would not move, but the same argument that makes cubic right for one reading makes it right for the other, and the two should not disagree about how a map is sampled. One named constant now carries the choice and the reason for both. Measured on two rotation datasets: 1.96 -> 2.20 and 2.69 -> 2.89 sigma. The maps, the R-factors and every other reported number are unchanged - only this reading of the map moves, so it no longer compares with what earlier versions printed. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_014zTy4Bpi4pPHw4bybf7q2R --- docs/CHANGELOG.md | 1 + rugnux/ModelValidation.cpp | 13 ++++++++----- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index dfcc7b37..b3bb229d 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.165 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. +* `MEAN_ATOM_DENSITY_SIGMA`, the mean 2Fo-Fc density at the model's atom centres that `rugnux --model` reports, is read from the map by cubic rather than linear interpolation and comes out around a tenth higher as a result. The map, the R-factors and every other number are unchanged; only this reading of the map moves, so it is no longer comparable with the figure earlier versions printed. * `rugnux --model` writes an anomalous difference map (`_anom.ccp4`) when the merge kept the Bijvoet split, and names the ten model atoms it peaks highest on - `ANOMALOUS_SITE_01`..`_10` in the results report, and one line in the log - so a run says which atoms carry the anomalous signal and how many sigma each one is. * The rugnux results report records how the run was invoked, what it cost and what it ran on: `COMMAND_LINE=` is the command line as one shell-ready line, `WALL_TIME=` the whole invocation in seconds, and `GPU_COUNT=` / `GPU=` how many GPUs were visible and which ones. The total wall time is also printed on stdout, next to the processing time it is slightly larger than. * The rugnux results report has a `10. MODEL VALIDATION` section when `--model` was given: R-work / R-free with their reflection counts, the bulk-solvent and scale parameters, the mean 2Fo-Fc density at the atom centres, the reindexing operators and where the maps were written - or `MODEL_VALIDATION= NOT_PERFORMED` and the reason, when the model could not be used. `REPORT_VERSION` is 4; `WARNINGS` moves from section 10 to 11 and no existing key changed. diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 1072650d..3387c94d 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -49,6 +49,11 @@ double write_ccp4(const gemmi::Grid &map, const std::string &path) { return ccp4.hstats.rms; } +// Cubic, not the default linear, for reading a map at a point. The maps are sampled every d_min/3, +// and a peak that sharp read by trilinear interpolation comes out up to a quarter low - unevenly +// enough to reorder the anomalous sites. +constexpr int MAP_INTERPOLATION_ORDER = 3; + // How many anomalous sites the report names. The strongest few are what says whether the anomalous // signal is there and what carries it; a full site list is what the map file is for. constexpr size_t MAX_ANOMALOUS_SITES = 10; @@ -285,7 +290,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & for (gemmi::Model &m : st.models) for (gemmi::Chain &ch : m.chains) for (gemmi::Residue &r : ch.residues) - for (gemmi::Atom &a : r.atoms) { s += grid2fofc.interpolate_value(a.pos); ++n; } + for (gemmi::Atom &a : r.atoms) { s += grid2fofc.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER); ++n; } result.mean_atom_density_sigma = (n > 0 && rms2 > 0) ? (s / n) / rms2 : 0; } @@ -347,12 +352,10 @@ ModelValidationResult ValidateAgainstModel(const std::vector & for (gemmi::Atom &a : r.atoms) { if (a.is_hydrogen()) // hydrogen scatters no anomalous signal continue; - // Cubic, not the default linear: the map is sampled every d_min/3, and a - // peak that sharp read by trilinear interpolation comes out up to a quarter - // low - unevenly enough to reorder the sites. sites.push_back({fmt::format("{} {} {}{}", a.name, r.name, ch.name, r.seqid.str()), - rms > 0 ? grid.interpolate_value(a.pos, 3) / rms : 0.0}); + rms > 0 ? grid.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER) / rms + : 0.0}); } std::sort(sites.begin(), sites.end(), [](const auto &x, const auto &y) { return x.sigma > y.sigma; });