diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 26414f48f..ac8c70426 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -601,6 +601,14 @@ a cutoff, found by flood fill as Coot's "unmodelled blobs" finds them, in GEMMI' (`find_blobs_by_flood_fill`): P. Emsley, B. Lohkamp, W. G. Scott and K. Cowtan, "Features and development of Coot" (2010), Acta Cryst. D66, 486-501 [doi:10.1107/S0907444910007493](https://doi.org/10.1107/S0907444910007493). +**Water build** — the waters the refinement adds and removes follow the documented ordered-solvent defaults of +phenix.refine: difference-map peaks from 3 σ, hydrogen-bond partners to 3.2 Å and picking only to 2.8 Å, with +the waters added, filtered and refined between macro-cycles: +Afonine et al. (2012, above); +the distances to the model that Coot's Find Waters keeps (2.4-3.2 Å), from the Coot manual and Emsley et al. +(2010, above); and the removal of a water in negative difference density at -3 σ, the cut Echols et al. (2014, +below) use to exclude a site. + **Solvent audit** — the waters and ions of the refined model are scored by their surroundings with the stereochemical criteria of the surveys of ions at nucleic acids by P. Auffinger's group: the coordination bands 1.9-2.3 / 2.3-2.6 / 2.6-3.2 Å, the Mg²⁺ exclusion zone and the clash shells of N. Naleem, A. Henning-Knechtel, diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index 7a6f9336b..0c9d74b4f 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -332,6 +332,40 @@ is read only from the residual difference density, which a B factor or an occupa values are not yet compared with the neighbouring atoms'. Per site, everything above is in `_solvent.json`; the counts are the `REFINE_SOLVENT_` keys of the report. +**The water build.** After the refinement, and to 2.8 Å (phenix.refine's ordered-solvent limit), waters are +built and removed between refinements, as phenix.refine does in its macro-cycles: a round reads the refined +model's maps, changes its waters, and refines the changed model again from the start; the round is kept where +R-free is lower than before it, and the rounds stop when one changes nothing or after three. The atom +list of a refinement is fixed, so each round is a refinement of its own. Waters only: no ion is placed or +relabelled, and `--refine-score-only` builds nothing. + +- *Candidates* are the local maxima of mFo-DFc from 3 σ (phenix.refine's peak cutoff), each symmetry-equivalent + set once, moved to the copy nearest the model. A peak closer than 1.2 Å to a model atom is that atom's own + residual and is passed over. +- *Not solvent* - reported in `_solvent.json` with the reason, never filled: a peak 1.2-1.9 Å from a model + atom (a missing or displaced atom); one beside a model atom (within 2.6 Å) that sits in mFo-DFc of -3 σ or + lower (the atom is displaced); one within 3 Å of where a truncated residue's missing atoms would be (an atom the + dictionary bonds to an atom the model lacks); the centre of three or more other peaks at 1.3-1.7 Å (a + tetrahedral group, sulfate or phosphate) or any peak with another 1.1-1.9 Å away (atoms of one molecule); + elongated density (the rms length of the density above 3 σ within 1.6 Å along its longest axis over its shortest + above 3); and 2mFo-DFc above every S and P of the model (a heavier atom). +- *A water* is placed where the peak has a hydrogen-bond partner - an N or O that can donate or accept at + 2.2-3.4 Å (phenix.refine's 3.2 Å, widened for the coordinate error), a metal within 2.6 Å (its ligand), or a + halide at 2.9-3.4 Å - and nothing closer than 2.2 Å to an N or O, nor than 3.0 Å to any other atom (Coot keeps + waters 2.4-3.2 Å from the model; a little wider here for the peak's position), strongest first, none within + 2.2 Å of one already placed. It goes into the chain of its nearest atom, numbered on from that chain, with + occupancy 1 and the mean B of its partners. +- *Removed* is any water whose mFo-DFc at its centre is -3 σ or lower (the cut Echols et al. use), and a water the + build itself added that has no N, O or ion within 3.4 Å any more or whose B reached the refinement's ceiling. + phenix.refine's other filter, 2mFo-DFc below 1 σ at a refined water, is not used: on the refine arm it raised + R-free wherever it was tried, on the deposited waters (five sets of five) and on the waters built here, and so + did taking out deposited waters without a partner. A weak water still explains density, and a deposited one + carries its depositor's judgement. + +`REFINE_SOLVENT_ADDED` and `REFINE_SOLVENT_REMOVED` count the waters added and still there and those of the model +as given taken out; `REFINE_SOLVENT_NOT_SOLVENT` the not-solvent peaks of the final model's difference map; +`REFINE_SOLVENT_BUILD_ROUNDS` the rounds kept. The model as given (`REFINE_R_FREE_START`, the `REFINE_SCORE_` terms) and the placement stage are the first refinement's; the refined model, its R-factors, geometry, maps and blobs are the last kept round's. + ### 15.11 GPU and CPU On a CUDA build the target, the gradient and the restraints are evaluated on the card - the density gather, @@ -355,8 +389,8 @@ What the refinement does not do yet - each a known gap, not a design choice: refined. - **NCS restraints**: copies are refined independently. - **Twin refinement**: the target is not a twinned one, also where the twin hypothesis' merge was written. -- **Water building and ion identification**: waters are neither picked nor removed and no ion is relabelled; - the blob finder points at candidates and the solvent audit (§15.10) scores the waters and ions there are. +- **Ion identification**: waters are built and removed (§15.10), but no ion is placed or relabelled; the solvent + audit scores the waters and ions there are. - **Alternative conformations and occupancies**: existing conformers are refined in place; no occupancy is refined and no conformer is added. - **Rotamer and peptide-flip search**, and **Ramachandran restraints**. @@ -385,6 +419,7 @@ What the refinement does not do yet - each a known gap, not a design choice: - A. T. Brünger, "Free R value: a novel statistical quantity for assessing the accuracy of crystal structures", *Nature* **355** (1992), 472-475 (R-free). - P. Emsley, B. Lohkamp, W. G. Scott & K. Cowtan, "Features and development of Coot", *Acta Cryst.* **D66** (2010), 486-501 (unmodelled blobs, §15.9). - N. Naleem, A. Henning-Knechtel, S. Kirmizialtin & P. Auffinger, "Cat_Wiz: a stereochemistry-guided toolkit for locating, diagnosing, and annotating Mg2+ ions in RNA structures", *Nucleic Acids Res.* **54** (2026), gkag120; F. Leonarski, L. D'Ascenzo & P. Auffinger, *Nucleic Acids Res.* **45** (2017), 987-1004 and *RNA* **25** (2019), 173-192; F. Leonarski, A. Henning-Knechtel, S. Kirmizialtin, E. Ennifar & P. Auffinger, *Nucleic Acids Res.* **53** (2025), gkae1148 (the coordination bands, distances and exclusion zone of the solvent audit, §15.10). +- P. V. Afonine et al., "Towards automated crystallographic structure refinement with phenix.refine", *Acta Cryst.* **D68** (2012), 352-367, and P. Emsley, B. Lohkamp, W. G. Scott & K. Cowtan, "Features and development of Coot", *Acta Cryst.* **D66** (2010), 486-501 (the peak cutoff, the partner and clash distances and the resolution limit of the water build, §15.10). - N. Echols et al., "Automated identification of elemental ions in macromolecular crystal structures", *Acta Cryst.* **D70** (2014), 1104-1114, and N. Morshed, N. Echols & P. D. Adams, "Using support vector machines to improve elemental ion identification in macromolecular crystal structures", *Acta Cryst.* **D71** (2015), 1147-1158 (the ligand sphere and the halide rules, §15.10). - H. Zheng et al., "CheckMyMetal: a macromolecular metal-binding validation tool", *Acta Cryst.* **D73** (2017), 223-233 (the vector sum, §15.10); O. Carugo, "Buried chloride stereochemistry in the Protein Data Bank", *BMC Struct. Biol.* **14** (2014), 19 (the chloride's shell, §15.10). - I. D. Brown & D. Altermatt, *Acta Cryst.* **B41** (1985), 244-247, and N. E. Brese & M. O'Keeffe, "Bond-valence parameters for solids", *Acta Cryst.* **B47** (1991), 192-197 (the bond-valence sums, §15.10). diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 28c5b0eba..e2011cc93 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -137,7 +137,7 @@ myrun_model.cif / .pdb the model as placed, in the cell and space group of my myrun_refined_model.cif the refined model (and .pdb) myrun_refined_maps.mtz map coefficients on the refined model myrun_blobs.json the unexplained difference-density blobs -myrun_solvent.json the refined model's waters and ions, each scored by its surroundings +myrun_solvent.json the refined model's waters and ions, each scored by its surroundings, and the waters the build added and removed ``` A stills run writes the same merged files (its `.hkl` holds the merged reflections) and, on a diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index dbe130610..8865f17d4 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -1206,6 +1206,11 @@ Where a model was placed and refined ([Refinement](RUGNUX_MODEL_AND_REFERENCE.md | `REFINE_SOLVENT_LABELS`, `REFINE_SOLVENT_BEST` | How many of them carry each label (`HOH 265 MG 42`), and how many each species is the best-supported one for | | `REFINE_SOLVENT_DISAGREE`, `REFINE_SOLVENT_AMBIGUOUS` | The sites whose best-supported species is not their label (water, Mg, Na, K, Ca or Cl), not ambiguously; the sites where the runner-up is within 2 of the best | | `REFINE_SOLVENT_IONS`, `REFINE_SOLVENT_IONS_DISAGREE` | The same for the sites not labelled water | +| `REFINE_SOLVENT_ADDED` | Waters the water build added and still in the refined model (CPU_DATA_ANALYSIS_REFINEMENT §15.10). Present where the build ran: a restrained refinement to 2.8 A, not `--refine-score-only` | +| `REFINE_SOLVENT_REMOVED` | Waters of the model as given that the build took out (mFo-DFc at -3 sigma or lower at the water) | +| `REFINE_SOLVENT_NOT_SOLVENT` | Difference-map peaks of the final model that cannot be a single solvent atom (bonded to the model, beside a displaced atom, a truncated side chain, a cluster or elongated density, denser than S and P): reported, not filled; each with its reason in `_solvent.json` | +| `REFINE_SOLVENT_BUILD_ROUNDS` | The build's rounds kept, each a refinement of its own kept where it lowered R-free | +| `REFINE_SOLVENT_BUILD_SECONDS` | The build's time, its refinements included (REFINE_SECONDS includes it too) | A prose paragraph after the blob keys lists the five strongest blobs with their reading, peak, volume, position and nearest atom; all of them are in `_blobs.json`, strongest first, one object per blob: `kind`, @@ -1220,7 +1225,11 @@ in model order, one object per site: `label`, `chain`, `seqid`, `altloc`, `xyz`, `donor_acceptors` (closer than 3.4 Å), `min_ligand_distance_A`, `nearest_positive_A`, `nearest_negative_A`, `positive_within_3.5`, `negative_within_3.2`, `nearest_ion_A`, `ions_within_4.2`, `ions_within_6`, `bvs` and `nvecsum` (per cation: `MG`, `NA`, `K`, `CA`), `map_sigma`, `difference_sigma`, `anomalous_sigma` (null where -there is no such map), `score` (per species), `best`, `ambiguous`, `agrees` and `neighbours`. +there is no such map), `score` (per species), `best`, `ambiguous`, `agrees` and `neighbours`. Where the water +build ran, a `build` object follows: `rounds`, `added`, `removed`, `not_solvent`, `stop` (why it stopped), +`seconds`, and `edits`, one per water added or removed in the rounds kept and per not-solvent peak of the final +model: `kind` (`ADDED`, `REMOVED`, `NOT_SOLVENT`), `round`, `site` (the water's name, empty for a peak), +`xyz`, `difference_sigma`, `map_sigma`, `b_iso`, `elongation` and `reason`. **The score of the model as given** follows, every term under `REFINE_SCORE_`: `REFINE_SCORE_R_WORK`, `REFINE_SCORE_R_FREE`, `REFINE_SCORE_N_WORK` and `REFINE_SCORE_N_FREE` by default, and with `--developer` the diff --git a/image_analysis/structure_refinement/CMakeLists.txt b/image_analysis/structure_refinement/CMakeLists.txt index a3d9334e7..d7a5b82ef 100644 --- a/image_analysis/structure_refinement/CMakeLists.txt +++ b/image_analysis/structure_refinement/CMakeLists.txt @@ -24,6 +24,7 @@ ADD_LIBRARY(JFJochStructureRefinement STATIC Refinement.cpp Refinement.h DensityBlobs.cpp DensityBlobs.h SolventSites.cpp SolventSites.h + SolventBuild.cpp SolventBuild.h $<$:RefinementGPU.cpp> $<$:RefinementGPU.cu> RefinementGPU.h RefinementGPUEngine.h diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index 89c02e5f1..ed119cd3f 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -1225,7 +1226,10 @@ std::vector RunTrials(const std::vector &trial_slots, } // namespace -RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector &merged, +namespace { + +// One refinement of a model whose atoms stay as they are (RefineModel adds and removes waters between calls). +RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector &merged, const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double wavelength_A, Logger &logger, const RefinementSettings &settings) { const auto start = std::chrono::steady_clock::now(); @@ -1309,8 +1313,8 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector(ctx.st); @@ -1569,3 +1573,87 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector &merged, + const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double wavelength_A, + Logger &logger, const RefinementSettings &settings) { + const auto start = std::chrono::steady_clock::now(); + RefinementResult result = RefineOnce(placed, merged, cell, sg, wavelength_A, logger, settings); + // The water build: a round of waters added and removed on the refined model's maps (SolventBuild.h), then a + // refinement of the model so changed, kept where its R-free is lower; until a round changes nothing. + const bool build = settings.solvent && settings.restrained && !settings.score_only && result.ok && + result.restrained && !result.cancelled && result.d_min <= SOLVENT_BUILD_D_MIN; + if (build) { + const auto build_start = std::chrono::steady_clock::now(); + SolventBuildResult b; + b.ran = true; + std::set added; // the waters added and still there + std::vector not_solvent; // the not-solvent peaks of the model last picked on + for (int round = 1;; round++) { + const SolventMaps maps{result.map.get(), result.map_rms, result.difference_map.get(), + result.difference_map_rms, nullptr, 0}; + SolventRound picked = PickSolvent(*result.model, maps, result.solvent, added, result.b_ceiling, round); + not_solvent.clear(); + for (const SolventEdit &e : picked.edits) + if (e.kind == SolventEdit::Kind::NotSolvent) + not_solvent.push_back(e); + std::string refused; + for (const auto &[why, n] : picked.refused) + refused += fmt::format("{}{} {}", refused.empty() ? "" : ", ", n, why); + logger.Info("Solvent build, round {}: {} waters to add, {} to remove, {} peaks not solvent{}", round, + picked.added, picked.removed, picked.not_solvent, + refused.empty() ? std::string() : "; peaks refused: " + refused); + b.not_solvent = picked.not_solvent; + if (picked.added + picked.removed == 0) { + b.stop = "nothing to change"; + break; + } + if (round > SOLVENT_BUILD_ROUNDS) { + b.stop = fmt::format("{} rounds", SOLVENT_BUILD_ROUNDS); + break; + } + gemmi::Structure changed = *result.model; + ApplySolventEdits(changed, picked.edits); + RefinementResult next = RefineOnce(changed, merged, cell, sg, wavelength_A, logger, settings); + if (!next.ok || !next.restrained || next.cancelled || !(next.r_free < result.r_free)) { + logger.Info("Solvent build, round {}: not kept (R-free {:.4f}, before {:.4f})", round, next.r_free, + result.r_free); + b.stop = fmt::format("round {} raised R-free", round); + break; + } + for (const SolventEdit &e : picked.edits) { + if (e.kind == SolventEdit::Kind::Added) { + added.insert(e.site); + b.edits.push_back(e); + } else if (e.kind == SolventEdit::Kind::Removed) { + if (added.erase(e.site) == 0) + b.removed++; + b.edits.push_back(e); + } + } + b.rounds = round; + // the model as given and its placement stay those of the first refinement + next.r_work_start = result.r_work_start; + next.r_free_start = result.r_free_start; + next.score_start = std::move(result.score_start); + next.geometry_start = result.geometry_start; + next.rigid_bodies = result.rigid_bodies; + next.group_b = result.group_b; + next.r_work_placement = result.r_work_placement; + next.r_free_placement = result.r_free_placement; + result = std::move(next); + } + b.added = static_cast(added.size()); + b.edits.insert(b.edits.end(), not_solvent.begin(), not_solvent.end()); + b.seconds = std::chrono::duration(std::chrono::steady_clock::now() - build_start).count(); + result.seconds = std::chrono::duration(std::chrono::steady_clock::now() - start).count(); + logger.Info("Solvent build: {} waters added, {} removed in {} round{} kept ({}), {} peaks not solvent; " + "R-free {:.4f}; {:.1f} s", b.added, b.removed, b.rounds, b.rounds == 1 ? "" : "s", b.stop, + b.not_solvent, result.r_free, b.seconds); + result.solvent_build = std::move(b); + } + result.map.reset(); + return result; +} diff --git a/image_analysis/structure_refinement/Refinement.h b/image_analysis/structure_refinement/Refinement.h index c33d2e1e4..cefde1799 100644 --- a/image_analysis/structure_refinement/Refinement.h +++ b/image_analysis/structure_refinement/Refinement.h @@ -20,6 +20,7 @@ #include "DensityBlobs.h" #include "ModelValidation.h" // ModelMaps #include "RefinementRestraints.h" // GeometryStatistics +#include "SolventBuild.h" #include "SolventSites.h" class Logger; @@ -58,6 +59,9 @@ struct RefinementSettings { std::vector final_ladder = {0.5, 1.0, 2.0}; double max_weight = 1.0 / 8; // the positions' X-ray weight at most, times the balanced one bool riding_hydrogens = true; // the dictionary's hydrogens added, riding on their parents + // Waters added at the difference map's peaks and removed where the maps no longer hold them, between refinements + // (SolventBuild.h), to SOLVENT_BUILD_D_MIN. + bool solvent = true; bool gpu = true; size_t nthreads = 1; size_t gpu_resident_bytes = 0; // held on the card by what runs beside @@ -159,9 +163,13 @@ struct RefinementResult { // The refined model's waters and ions and the species their surroundings support (AnalyseSolventSites): an // audit, nothing in the model is changed by it. std::vector solvent; + double b_ceiling = 0; // the highest B the refinement allowed + SolventBuildResult solvent_build; // the waters added and removed (settings.solvent) }; -// `placed`: the model in the data's cell and space group. `merged`: the merge, with French-Wilson amplitudes +// `placed`: the model in the data's cell and space group. With settings.solvent (and restrained, not score_only), +// the refinement is followed by up to SOLVENT_BUILD_ROUNDS rounds of the water build, each a new refinement of the +// model with its waters changed, kept where it lowers R-free. `merged`: the merge, with French-Wilson amplitudes // and free flags, in the same group. RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector &merged, const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double wavelength_A, diff --git a/image_analysis/structure_refinement/SolventBuild.cpp b/image_analysis/structure_refinement/SolventBuild.cpp new file mode 100644 index 000000000..755eb675e --- /dev/null +++ b/image_analysis/structure_refinement/SolventBuild.cpp @@ -0,0 +1,464 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "SolventBuild.h" + +#include +#include +#include + +#include + +#include "gemmi/asumask.hpp" +#include "gemmi/math.hpp" +#include "gemmi/neighbor.hpp" + +#include "MonomerLibrary.h" +#include "RefinementRestraints.h" // AtomsNear + +namespace { + +// The ordered-solvent defaults of phenix.refine: peaks of mFo-DFc from 3 sigma, hydrogen-bond partners to 3.2 A +// (taken here to 3.4 A, the outer edge of the coordinate error). +// Following Afonine et al. (2012) Acta Cryst. D68, 352-367 +constexpr double PEAK_SIGMA = 3.0; +constexpr double PARTNER_MAX = 3.4; // A +// A water no nearer than 2.2 A to an N or O (one that close is a cation's ligand or an error), nor than 3.0 A +// to a C or any other atom that makes no hydrogen bond; Coot's Find Waters keeps waters 2.4-3.2 A from the +// model, Emsley et al. (2010) Acta Cryst. D66, 486-501, here a little wider for the peak's own position error. +constexpr double POLAR_CLASH = 2.2; // A +constexpr double APOLAR_CLASH = 3.0; // A +// A: a water this near a metal is its ligand or in its shell (K+ holds its waters at up to 3.2 A, Cat_Wiz's band) +constexpr double METAL_LIGAND_MAX = 3.2; +constexpr double HALIDE_PARTNER_MIN = 2.9; // A: a water donates to a halide from here +// Nearer than ON_ATOM to a model atom, the peak is that atom's own residual (its B, its anisotropy) and no candidate; +// out to BONDED it is bonded to it (a missing or displaced atom), not a water. +constexpr double ON_ATOM = 1.2, BONDED = 1.9; // A +// A difference peak beside a model atom that sits in negative difference density: the atom is displaced. +// A water in negative density at its centre has more electrons modelled than the data show. +// Following Echols et al. (2014) Acta Cryst. D70, 1104-1114 (their -3 sigma) +constexpr double NEGATIVE_SIGMA = -3.0; +constexpr double NEGATIVE_REACH = 2.6; // A: the model atom the peak is beside +// Peaks of one molecule: local maxima this far apart are bonded atoms, nearer they are one peak. +constexpr double CLUSTER_MIN = 1.1, CLUSTER_MAX = 1.9; // A +// A tetrahedral cluster: three or more peaks 1.3-1.7 A from one (sulfate, phosphate). +constexpr double TETRAHEDRAL_MIN = 1.3, TETRAHEDRAL_MAX = 1.7; // A +// The density a peak is measured in for its shape: the points above the cutoff within this distance. +constexpr double SHAPE_RADIUS = 1.6; // A +// rms along the longest axis over the shortest; the peaks at deposited waters of five refine-arm sets (stripped and +// rebuilt) reach 3.0 at their 90th percentile, the measure being coarse on a few grid points. +constexpr double ELONGATION_MAX = 3.0; +// A peak this near the model's truncated residues' last atoms is where their missing atoms would be. +constexpr double TRUNCATED_REACH = 3.0; // A + +bool Polar(gemmi::El el) { + return el == gemmi::El::N || el == gemmi::El::O; +} + +bool IsHalogen(gemmi::El el) { + return el == gemmi::El::F || el == gemmi::El::Cl || el == gemmi::El::Br || el == gemmi::El::I; +} + +std::string ResidueName(const gemmi::Residue &r, const std::string &chain) { + return r.name + " " + chain + r.seqid.str(); +} + +std::string AtomName(const gemmi::Residue &r, const gemmi::Atom &a, const std::string &chain) { + return a.name + " " + ResidueName(r, chain); +} + +// Each grid point within `radius` of `p`, its position and value, the cell's periodicity taken care of. +template +void ForPointsNear(const gemmi::Grid &g, const gemmi::Position &p, double radius, F &&f) { + const gemmi::Fractional fp = g.unit_cell.fractionalize(p); + const int u0 = static_cast(std::round(fp.x * g.nu)), v0 = static_cast(std::round(fp.y * g.nv)), + w0 = static_cast(std::round(fp.z * g.nw)); + const int du = static_cast(std::ceil(radius / g.spacing[0])) + 1; + const int dv = static_cast(std::ceil(radius / g.spacing[1])) + 1; + const int dw = static_cast(std::ceil(radius / g.spacing[2])) + 1; + for (int w = w0 - dw; w <= w0 + dw; w++) + for (int v = v0 - dv; v <= v0 + dv; v++) + for (int u = u0 - du; u <= u0 + du; u++) { + const gemmi::Position q = g.get_position(u, v, w); + if (q.dist(p) <= radius) + f(q, g.get_value(u, v, w), u, v, w); + } +} + +bool IsLocalMaximum(const gemmi::Grid &g, int u, int v, int w) { + const float value = g.get_value(u, v, w); + for (int dw = -1; dw <= 1; dw++) + for (int dv = -1; dv <= 1; dv++) + for (int du = -1; du <= 1; du++) + if ((du || dv || dw) && g.get_value(u + du, v + dv, w + dw) > value) + return false; + return true; +} + +struct Peak { + gemmi::Position pos; + double sigma; +}; + +// The local maxima of `g` (in rms `rms`) from `cutoff`, each symmetry-equivalent set once, strongest first. +std::vector DifferencePeaks(const gemmi::Grid &g, double rms, double cutoff) { + std::vector out; + const std::vector asu = gemmi::get_asu_mask(g); + for (int w = 0; w < g.nw; w++) + for (int v = 0; v < g.nv; v++) + for (int u = 0; u < g.nu; u++) { + const size_t idx = g.index_q(u, v, w); + if (asu[idx] != 0 || g.data[idx] < cutoff * rms || !IsLocalMaximum(g, u, v, w)) + continue; + out.push_back({g.get_position(u, v, w), g.data[idx] / rms}); + } + std::sort(out.begin(), out.end(), [](const Peak &a, const Peak &b) { return a.sigma > b.sigma; }); + return out; +} + +// The copy of `p` nearest to the model's own atoms (not their symmetry mates); false where no atom is within +// `reach`. +bool NearestCopy(gemmi::NeighborSearch &ns, const gemmi::UnitCell &cell, double reach, gemmi::Position &p) { + gemmi::NeighborSearch::Mark *m = ns.find_nearest_atom(p, reach); + if (m == nullptr) + return false; + const gemmi::FTransform t = ns.get_image_transformation(m->image_idx); + const gemmi::Fractional fc = cell.fractionalize(p); + const gemmi::Fractional fi = t.apply(cell.fractionalize(m->pos)); + const gemmi::Fractional shifted(fc.x - std::round(fc.x - fi.x), fc.y - std::round(fc.y - fi.y), + fc.z - std::round(fc.z - fi.z)); + p = cell.orthogonalize(gemmi::FTransform(t.inverse()).apply(shifted)); + return true; +} + +// The longest over the shortest rms axis of the density above `cutoff` (in map units) within SHAPE_RADIUS of +// `p`, each point weighted by how far it is above; 1 where there are too few points to tell. +double Elongation(const gemmi::Grid &g, const gemmi::Position &p, double cutoff) { + std::vector> points; + double sum = 0; + gemmi::Vec3 mean; + ForPointsNear(g, p, SHAPE_RADIUS, [&](const gemmi::Position &q, float value, int, int, int) { + if (value > cutoff) { + points.emplace_back(q, value - cutoff); + sum += value - cutoff; + mean += q * (value - cutoff); + } + }); + if (points.size() < 6) + return 1.0; + mean /= sum; + double c[6] = {0, 0, 0, 0, 0, 0}; + for (const auto &[q, wt] : points) { + const gemmi::Vec3 d = gemmi::Vec3(q) - mean; + c[0] += wt * d.x * d.x; + c[1] += wt * d.y * d.y; + c[2] += wt * d.z * d.z; + c[3] += wt * d.x * d.y; + c[4] += wt * d.x * d.z; + c[5] += wt * d.y * d.z; + } + const gemmi::SMat33 m{c[0] / sum, c[1] / sum, c[2] / sum, c[3] / sum, c[4] / sum, c[5] / sum}; + std::array ev = m.calculate_eigenvalues(); + std::sort(ev.begin(), ev.end()); + return ev[0] > 0 ? std::sqrt(ev[2] / ev[0]) : INFINITY; +} + +// The atoms of residues the dictionary describes with heavy atoms the model lacks, that the dictionary bonds to +// one of the missing atoms: the ends a truncated side chain stops at. +std::set TruncatedEnds(const gemmi::Model &model) { + std::set out; + const gemmi::MonLib &monlib = EmbeddedMonomerLibrary(); + for (const gemmi::Chain &ch : model.chains) + for (const gemmi::Residue &r : ch.residues) { + auto cc = monlib.monomers.find(r.name); + if (cc == monlib.monomers.end() || r.is_water()) + continue; + auto missing = [&](const std::string &name) { + auto a = cc->second.find_atom(name); + return a != cc->second.atoms.end() && !a->is_hydrogen() && name != "OXT" && name != "OP3" && + r.find_atom(name, '*') == nullptr; + }; + for (const gemmi::Restraints::Bond &b : cc->second.rt.bonds) + for (int k = 0; k < 2; k++) { + const std::string &here = k == 0 ? b.id1.atom : b.id2.atom, &there = k == 0 ? b.id2.atom : b.id1.atom; + if (missing(there)) + if (const gemmi::Atom *a = r.find_atom(here, '*')) + out.insert(a); + } + } + return out; +} + +} // namespace + +const char *SolventEditKindName(SolventEdit::Kind kind) { + switch (kind) { + case SolventEdit::Kind::Added: return "ADDED"; + case SolventEdit::Kind::Removed: return "REMOVED"; + case SolventEdit::Kind::NotSolvent: return "NOT_SOLVENT"; + } + return "?"; +} + +SolventRound PickSolvent(const gemmi::Structure &st, const SolventMaps &maps, const std::vector &sites, + const std::set &built, double b_ceiling, int round) { + SolventRound out; + if (st.models.empty() || maps.difference == nullptr || maps.map == nullptr || !(maps.difference_rms > 0) || + !(maps.map_rms > 0)) + return out; + const gemmi::Grid &diff = *maps.difference, &map = *maps.map; + gemmi::Model model = st.models[0]; // NeighborSearch takes it non-const + gemmi::UnitCell cell = st.cell; + cell.set_cell_images_from_spacegroup(st.find_spacegroup()); + gemmi::NeighborSearch ns(model, cell, 5.0); + ns.populate(false); + const gemmi::MonLib &monlib = EmbeddedMonomerLibrary(); + const std::set truncated = TruncatedEnds(model); + + // The densest S and P of the model: a peak denser than every one of them is no light atom. + double heaviest = -INFINITY; + for (const gemmi::Chain &ch : model.chains) + for (const gemmi::Residue &r : ch.residues) + for (const gemmi::Atom &a : r.atoms) + if (a.element == gemmi::El::S || a.element == gemmi::El::P) + heaviest = std::max(heaviest, map.interpolate_value(a.pos) / maps.map_rms); + if (!std::isfinite(heaviest)) + heaviest = INFINITY; + + // The waters the maps no longer support: any in negative difference density, and those built here also where they + // have no partner left or their B reached the ceiling. phenix.refine's other filter, 2mFo-DFc below 1 sigma, is + // left out: on the refine arm it raised R-free wherever it was tried - on the depositions' waters (five sets of + // five) and on the waters built here (the next round rejected in five of six) - and so did taking out the + // depositions' waters without a partner (1K4C's 38). A weak water still explains density. + for (const SolventSite &s : sites) { + if (!SolventLabelIsWater(s)) + continue; + std::string reason; + const std::string name = s.label + " " + s.chain + s.seqid; + const bool own = built.count(name) > 0; + if (s.difference_sigma <= NEGATIVE_SIGMA) + reason = fmt::format("mFo-DFc {:.1f} sigma", s.difference_sigma); + else if (own && s.n_below_19 + s.n_19_23 + s.n_23_26 + s.n_26_32 + s.n_32_34 == 0 && !(s.nearest_ion < PARTNER_MAX)) + reason = fmt::format("no N, O or ion within {} A", PARTNER_MAX); + else if (own && s.b_iso >= b_ceiling - 0.5) + reason = fmt::format("B {:.0f} at the refinement's ceiling", s.b_iso); + if (reason.empty()) + continue; + SolventEdit e; + e.kind = SolventEdit::Kind::Removed; + e.round = round; + e.site = name; + e.chain = s.chain; + e.seqid = gemmi::SeqId(s.seqid); + e.altloc = s.altloc; + e.pos = s.pos; + e.difference_sigma = s.difference_sigma; + e.map_sigma = s.map_sigma; + e.b_iso = s.b_iso; + e.reason = reason; + out.edits.push_back(e); + out.removed++; + } + // A removed water is no partner of a new one, nor in its way. + std::vector removed; + for (const SolventEdit &e : out.edits) + removed.push_back(e.pos); + auto is_removed = [&](const gemmi::Position &p) { + return std::any_of(removed.begin(), removed.end(), [&](const gemmi::Position &q) { return q.dist(p) < 0.01; }); + }; + + // The difference map's peaks, strongest first. + std::vector added; + for (Peak peak : DifferencePeaks(diff, maps.difference_rms, PEAK_SIGMA)) { + if (!NearestCopy(ns, cell, PARTNER_MAX + 0.1, peak.pos)) + continue; // nothing of the model near enough to hold it + SolventEdit e; + e.round = round; + e.pos = peak.pos; + e.difference_sigma = peak.sigma; + e.map_sigma = map.interpolate_value(peak.pos) / maps.map_rms; + e.elongation = Elongation(diff, peak.pos, PEAK_SIGMA * maps.difference_rms); + auto refuse = [&](const std::string &why) { out.refused[why]++; }; + auto not_solvent = [&](const std::string &why) { + e.kind = SolventEdit::Kind::NotSolvent; + e.reason = why; + out.edits.push_back(e); + out.not_solvent++; + }; + + const std::vector near = AtomsNear(ns, model, peak.pos, PARTNER_MAX); + const ModelNeighbour *nearest = nullptr; + std::string clash, truncated_by; + std::vector partners; + double partner_b = 0; + for (const ModelNeighbour &nb : near) { + const gemmi::Chain &ch = model.chains[nb.chain]; + const gemmi::Residue &r = ch.residues[nb.residue]; + const gemmi::Atom &a = r.atoms[nb.atom]; + if (a.is_hydrogen() || (r.is_water() && is_removed(a.pos))) + continue; + if (nearest == nullptr) + nearest = &nb; + const double d = nb.distance; + const gemmi::El el = a.element.elem; + if (d < TRUNCATED_REACH && truncated_by.empty() && truncated.count(&a)) + truncated_by = ResidueName(r, ch.name); + bool partner = false; + if (a.element.is_metal()) + partner = d < METAL_LIGAND_MAX; + else if (Polar(el)) { + if (d < POLAR_CLASH && clash.empty()) + clash = AtomName(r, a, ch.name); + const AtomChemistry chem = ChemistryOf(monlib, r, a, false); + partner = d >= POLAR_CLASH && (chem.donor || chem.acceptor); + } else if (IsHalogen(el)) + partner = d >= HALIDE_PARTNER_MIN; + else if (d < APOLAR_CLASH && clash.empty()) + clash = AtomName(r, a, ch.name); + if (partner) { + partners.push_back(fmt::format("{} {:.2f}", AtomName(r, a, ch.name), d)); + partner_b += a.b_iso; + } + } + if (nearest == nullptr) + continue; + const gemmi::Chain &nch = model.chains[nearest->chain]; + const gemmi::Residue &nres = nch.residues[nearest->residue]; + const gemmi::Atom &natom = nres.atoms[nearest->atom]; + + if (nearest->distance < ON_ATOM) { + refuse("on a model atom"); + continue; + } + // Not a solvent atom: reported, never filled. + if (nearest->distance < BONDED) { + not_solvent(fmt::format("{:.2f} A from {}: a missing or displaced atom of the model", nearest->distance, + AtomName(nres, natom, nch.name))); + continue; + } + if (nearest->distance < NEGATIVE_REACH && + diff.interpolate_value(nearest->pos) / maps.difference_rms <= NEGATIVE_SIGMA) { + not_solvent(fmt::format("beside {} in negative density: the atom is displaced", + AtomName(nres, natom, nch.name))); + continue; + } + if (!truncated_by.empty()) { + not_solvent(fmt::format("where the missing atoms of {} would be", truncated_by)); + continue; + } + int bonded = 0, tetrahedral = 0; + ForPointsNear(diff, peak.pos, CLUSTER_MAX, [&](const gemmi::Position &q, float value, int u, int v, int w) { + const double d = q.dist(peak.pos); + if (d >= CLUSTER_MIN && value >= PEAK_SIGMA * maps.difference_rms && IsLocalMaximum(diff, u, v, w)) { + bonded++; + tetrahedral += d >= TETRAHEDRAL_MIN && d <= TETRAHEDRAL_MAX; + } + }); + if (tetrahedral >= 3) { + not_solvent(fmt::format("the centre of {} peaks at {}-{} A: a tetrahedral group (sulfate, phosphate)", + tetrahedral, TETRAHEDRAL_MIN, TETRAHEDRAL_MAX)); + continue; + } + if (bonded > 0) { + not_solvent(fmt::format("{} other peak{} at {}-{} A: atoms of one molecule", bonded, bonded > 1 ? "s" : "", + CLUSTER_MIN, CLUSTER_MAX)); + continue; + } + if (e.elongation > ELONGATION_MAX) { + not_solvent(fmt::format("elongated ({:.1f}:1)", e.elongation)); + continue; + } + if (e.map_sigma > heaviest) { + not_solvent(fmt::format("2mFo-DFc {:.1f} sigma, above every S and P of the model ({:.1f})", e.map_sigma, + heaviest)); + continue; + } + + // A water, where it has a partner and nothing in its way. + if (!clash.empty()) { + refuse("too close to the model"); + continue; + } + if (partners.empty()) { + refuse("no hydrogen-bond partner"); + continue; + } + if (std::any_of(added.begin(), added.end(), [&](const gemmi::Position &q) { + return cell.find_nearest_image(peak.pos, q, gemmi::Asu::Any).dist() < POLAR_CLASH; + })) { + refuse("too close to a water added beside it"); + continue; + } + e.kind = SolventEdit::Kind::Added; + e.chain = nch.name; + e.b_iso = partner_b / partners.size(); + e.reason = "partners"; + for (size_t k = 0; k < partners.size(); k++) + e.reason += (k == 0 ? " " : ", ") + partners[k]; + out.edits.push_back(e); + out.added++; + added.push_back(peak.pos); + } + return out; +} + +void ApplySolventEdits(gemmi::Structure &st, std::vector &edits) { + gemmi::Model &model = st.models.at(0); + // A chain name can stand for several chains of the model (a polymer and its waters read as two); each is meant. + for (const SolventEdit &e : edits) { + if (e.kind != SolventEdit::Kind::Removed) + continue; + for (gemmi::Chain &ch : model.chains) { + if (ch.name != e.chain) + continue; + for (gemmi::Residue &r : ch.residues) + if (r.is_water() && r.seqid == e.seqid) + std::erase_if(r.atoms, [&](const gemmi::Atom &a) { return a.altloc == e.altloc && !a.is_hydrogen(); }); + std::erase_if(ch.residues, [](const gemmi::Residue &r) { + return r.is_water() && std::none_of(r.atoms.begin(), r.atoms.end(), [](const gemmi::Atom &a) { + return !a.is_hydrogen(); + }); + }); + } + } + for (SolventEdit &e : edits) { + if (e.kind != SolventEdit::Kind::Added) + continue; + // to the last chain of the name, numbered on from the last residue of any of them + gemmi::Chain *ch = nullptr; + int last = 0; + const gemmi::Residue *water = nullptr; + for (gemmi::Chain &c : model.chains) { + if (c.name != e.chain) + continue; + ch = &c; + for (const gemmi::Residue &r : c.residues) { + if (r.seqid.num != gemmi::SeqId::OptionalNum::None) + last = std::max(last, static_cast(r.seqid.num)); + if (r.is_water()) + water = &r; + } + } + if (ch == nullptr) + continue; + gemmi::Residue r; + r.name = "HOH"; + r.seqid = gemmi::SeqId(last + 1, ' '); + r.het_flag = 'H'; + r.entity_type = gemmi::EntityType::Water; + if (water != nullptr) { + r.subchain = water->subchain; + r.entity_id = water->entity_id; + } + gemmi::Atom a; + a.name = "O"; + a.element = gemmi::El::O; + a.pos = e.pos; + a.occ = 1; + a.b_iso = static_cast(e.b_iso); + r.atoms.push_back(a); + e.site = ResidueName(r, ch->name); + e.seqid = r.seqid; + ch->residues.push_back(r); + } +} diff --git a/image_analysis/structure_refinement/SolventBuild.h b/image_analysis/structure_refinement/SolventBuild.h new file mode 100644 index 000000000..1ab90dc49 --- /dev/null +++ b/image_analysis/structure_refinement/SolventBuild.h @@ -0,0 +1,75 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include + +#include "gemmi/model.hpp" + +#include "SolventSites.h" // SolventSite, SolventMaps + +// One round of the water build that RefineModel runs between refinements (Refinement.h): waters placed at the +// peaks of the difference map a refined model leaves, and the model's waters the maps no longer support taken +// out. Waters only - no ion is placed or relabelled here. +// +// A peak of mFo-DFc becomes a water where it has a hydrogen-bond partner and nothing too close; a peak that +// cannot be a single solvent atom - one bonded to the model, beside a negative peak on a model atom, where a +// truncated side chain's missing atoms would be, one of a cluster of peaks at bonded distances, elongated, or +// denser than the model's S and P - is reported as not solvent and left alone. +struct SolventEdit { + enum class Kind { Added, Removed, NotSolvent }; + Kind kind = Kind::NotSolvent; + int round = 0; + std::string site; // the water as modelled ("HOH A401"); empty for a peak not taken + std::string chain; // where an added water goes: the chain of the atom nearest to it + gemmi::SeqId seqid; // a removed water's residue + char altloc = '\0'; + gemmi::Position pos; + double difference_sigma = 0; // mFo-DFc at the site, in its rms + double map_sigma = 0; // 2mFo-DFc + double b_iso = 0; // an added water's starting B, a removed one's + double elongation = 0; // a peak's density: the longest axis over the shortest, rms + std::string reason; // why: the partners an added water has, what a removed or refused one lacks +}; + +const char *SolventEditKindName(SolventEdit::Kind kind); + +struct SolventRound { + std::vector edits; + int added = 0, removed = 0, not_solvent = 0; + std::map refused; // the peaks neither taken nor reported, by what they lack +}; + +// The resolution waters are built to at most, following the ordered-solvent defaults of phenix.refine +// (Afonine et al. (2012) Acta Cryst. D68, 352-367): no picking where the data stop short of 2.8 A. +constexpr double SOLVENT_BUILD_D_MIN = 2.8; + +// At most this many rounds are kept, each a refinement of its own. +constexpr int SOLVENT_BUILD_ROUNDS = 3; + +// What the build did, over the rounds kept. +struct SolventBuildResult { + bool ran = false; + int rounds = 0; // the rounds kept + int added = 0; // waters added and still in the model + int removed = 0; // waters of the model as given taken out + int not_solvent = 0; // peaks of the final model's difference map that are no solvent + std::vector edits; // the kept rounds' additions and removals, then the final not-solvent peaks + std::string stop; // why it stopped + double seconds = 0; +}; + +// A round on `model` (the refined model the maps were made from, in their cell and space group): `maps` its +// 2mFo-DFc and mFo-DFc (the anomalous one unused), `sites` its solvent audit (AnalyseSolventSites on the same +// maps), `built` the waters added by earlier rounds (by name, "HOH A401"), `b_ceiling` the highest B its refinement +// allowed. +SolventRound PickSolvent(const gemmi::Structure &model, const SolventMaps &maps, const std::vector &sites, + const std::set &built, double b_ceiling, int round); + +// The round's waters taken out and put in: each added one to the chain of its nearest atom, numbered on from that +// chain's last residue (its `site` is set to the name it gets). +void ApplySolventEdits(gemmi::Structure &st, std::vector &edits); diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 00c10b6cd..ec448e6b6 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -459,7 +459,8 @@ namespace { // The solvent audit of the refined model (SolventSites.h): REFINE_SOLVENT_ keys. Counts by label as modelled // and by the species the surroundings support best; a site disagrees where the best is another species, not // ambiguously. - void AddSolventKeys(ReportSection &s, const std::vector &sites, const std::string &output_prefix) { + void AddSolventKeys(ReportSection &s, const std::vector &sites, const SolventBuildResult &build, + const std::string &output_prefix) { std::map labels, best; constexpr int LISTED_IONS = 40; // the rest are in the JSON int disagree = 0, ambiguous = 0, ions = 0, ions_disagree = 0; @@ -500,6 +501,21 @@ namespace { " model is not changed). Per site in {}_solvent.json; the ions, label -> best-supported\n" " species, N/O in the bands 1.9-2.3/2.3-2.6/2.6-3.2 A, and their neighbours:{}", output_prefix, ion_rows.empty() ? std::string(" none") : ion_rows))); + if (build.ran) { + Add(s, KeyInt("REFINE_SOLVENT_ADDED", build.added)); + Add(s, KeyInt("REFINE_SOLVENT_REMOVED", build.removed)); + Add(s, KeyInt("REFINE_SOLVENT_NOT_SOLVENT", build.not_solvent)); + Add(s, KeyInt("REFINE_SOLVENT_BUILD_ROUNDS", build.rounds)); + Add(s, KeyReal("REFINE_SOLVENT_BUILD_SECONDS", build.seconds, "{:.1f}", true)); + Add(s, Prose(fmt::format( + " Waters were added at difference-map peaks with a hydrogen-bond partner and removed where the\n" + " maps no longer hold them, each round a refinement kept only where it lowered R-free: {}\n" + " added, {} of the model's taken out, in {} round{} ({}). {} peaks are not solvent (bonded to\n" + " the model, beside a displaced atom, a truncated side chain, a molecule's or elongated density,\n" + " denser than S and P): reported, not filled. Each, with its reason, in {}_solvent.json.", + build.added, build.removed, build.rounds, build.rounds == 1 ? "" : "s", build.stop, + build.not_solvent, output_prefix))); + } } // The refinement's keys (Refinement.h): REFINE_ of the refinement of a placed model and of --mode refine, with every term of the @@ -565,7 +581,7 @@ namespace { " R_FREE above remain the model as given. The refinement decided nothing about the\n" " processing. The strongest unexplained blobs of its difference map ({}_blobs.json):{}", output_prefix, output_prefix, top.empty() ? std::string(" none") : top))); - AddSolventKeys(s, rf.solvent, output_prefix); + AddSolventKeys(s, rf.solvent, rf.solvent_build, output_prefix); } if (rf.ok) { const RefinementScore &sc = rf.score_start; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 9c8277b65..660e1a954 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -419,6 +419,7 @@ struct RugnuxOutputs { std::shared_ptr refined_maps; // _refined_maps.mtz std::optional> blobs; // _blobs.json std::optional> solvent; // _solvent.json + std::optional solvent_build; // and the water build, in it std::optional unmerged; // _unmerged.mtz std::optional unmerged_partials; // _unmerged_partials.mtz std::string detector_jpeg; // _detector.jpg; empty: none drawn diff --git a/rugnux/RugnuxOutputFiles.cpp b/rugnux/RugnuxOutputFiles.cpp index 596fe2385..992b2c15b 100644 --- a/rugnux/RugnuxOutputFiles.cpp +++ b/rugnux/RugnuxOutputFiles.cpp @@ -26,8 +26,10 @@ void WriteDensityBlobs(const std::vector &blobs, const std::string logger.Info("Difference-map blobs written to {}", path); } -// The solvent audit of the refined model (SolventSites.h), one entry per site in model order. -void WriteSolventSites(const std::vector &sites, const std::string &path, Logger &logger) { +// The solvent audit of the refined model (SolventSites.h), one entry per site in model order, and what the water +// build (SolventBuild.h) added, removed and reported as not solvent, where it ran. +void WriteSolventSites(const std::vector &sites, const SolventBuildResult *build, const std::string &path, + Logger &logger) { std::ofstream f(path, std::ofstream::out | std::ofstream::trunc); if (!f) { logger.Warning("Cannot open {} for writing", path); @@ -62,7 +64,19 @@ void WriteSolventSites(const std::vector &sites, const std::string {"best", SolventSpeciesName(s.best)}, {"ambiguous", s.ambiguous}, {"agrees", SolventLabelAgrees(s)}, {"neighbours", s.neighbours}}); } - f << nlohmann::json{{"sites", list}}.dump(1) << "\n"; + nlohmann::json out{{"sites", list}}; + if (build != nullptr && build->ran) { + nlohmann::json edits = nlohmann::json::array(); + for (const SolventEdit &e : build->edits) + edits.push_back({{"kind", SolventEditKindName(e.kind)}, {"round", e.round}, {"site", e.site}, + {"xyz", {e.pos.x, e.pos.y, e.pos.z}}, {"difference_sigma", num(e.difference_sigma)}, + {"map_sigma", num(e.map_sigma)}, {"b_iso", num(e.b_iso)}, + {"elongation", num(e.elongation)}, {"reason", e.reason}}); + out["build"] = {{"rounds", build->rounds}, {"added", build->added}, {"removed", build->removed}, + {"not_solvent", build->not_solvent}, {"stop", build->stop}, {"seconds", build->seconds}, + {"edits", edits}}; + } + f << out.dump(1) << "\n"; logger.Info("Solvent sites written to {}", path); } @@ -92,7 +106,7 @@ void WriteOutputs(const RugnuxOutputs &o, const std::string &prefix, const Outpu if (o.blobs) WriteDensityBlobs(*o.blobs, prefix + "_blobs.json", logger); if (o.solvent) - WriteSolventSites(*o.solvent, prefix + "_solvent.json", logger); + WriteSolventSites(*o.solvent, o.solvent_build ? &*o.solvent_build : nullptr, prefix + "_solvent.json", logger); if (o.unmerged) { WriteMtz(*o.unmerged, prefix + "_unmerged.mtz"); logger.Info("Unmerged observations written to {}_unmerged.mtz", prefix); @@ -239,7 +253,7 @@ ReportDocument WriteRefineOutputs(const RefineOnlyResult &result, const std::str WritePlacedModel(*rf.model, cell, *sg, prefix + "_refined", logger); WriteModelMaps(*rf.maps, prefix + "_refined", logger); WriteDensityBlobs(rf.blobs, prefix + "_blobs.json", logger); - WriteSolventSites(rf.solvent, prefix + "_solvent.json", logger); + WriteSolventSites(rf.solvent, &rf.solvent_build, prefix + "_solvent.json", logger); } try { std::ofstream file(prefix + "_report.txt"); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 542a7462e..962a70049 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -5115,6 +5115,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { result.outputs->refined_maps = refined.maps; result.outputs->blobs = refined.blobs; result.outputs->solvent = refined.solvent; + result.outputs->solvent_build = refined.solvent_build; Note(fmt::format("refined{}: R-work {:.3f}, R-free {:.3f}, rmsd bonds {:.3f} A, {} blobs", refined.cancelled ? " (cut short)" : "", refined.r_work, refined.r_free, refined.geometry.rmsd_bonds, refined.blobs.size())); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 938e87c92..ebaa49208 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -256,7 +256,7 @@ void print_usage() { std::cout << " --scaling-iterations Cap on the per-frame scaling iterations; the loop stops earlier once the scales settle (default: 100)" << std::endl; std::cout << " -z, --reference Reference data of the same crystal form: an MTZ, or a PDB structure-factor mmCIF (e.g. 1abc-sf.cif, first merged block with a usable column; R-free set from _refln.status f/o, else pdbx_r_free_flag). Gzipped or not, recognised by content. --reference-mtz is the same option" << std::endl; std::cout << " --reference-column