RESULT-CHANGING: Refinement: waters built at difference-map peaks and taken out of negative density, between refinements

RefineModel now runs the refinement (RefineOnce, the former body), then up to
three rounds of a water build (SolventBuild.{h,cpp}): the refined model's
mFo-DFc peaks from 3 sigma (each symmetry-equivalent set once, moved to the
copy nearest the model) become waters where they have a hydrogen-bond partner
(an N/O donor or acceptor at 2.2-3.4 A, an ion within 3.2 A, a halide at
2.9-3.4 A) and nothing closer than 2.2 A to N/O or 3.0 A to anything else; the
changed model is refined again from the start and the round is kept only where
R-free is lower. Peaks that cannot be one solvent atom are reported, never
filled: 1.2-1.9 A from a model atom, beside a model atom in -3 sigma density,
where a truncated residue's missing atoms would be, a tetrahedral cluster or
bonded peaks, elongated density, denser than every S and P. Removal: any water
in mFo-DFc <= -3 sigma; built waters also when partnerless or at the B ceiling.
On by default (RefinementSettings::solvent) to 2.8 A; not with
--refine-score-only. REFINE_SOLVENT_ADDED/REMOVED/NOT_SOLVENT/BUILD_ROUNDS/
BUILD_SECONDS; per-site reasons in _solvent.json "build". The model as given
and the placement stage keep the first refinement's numbers.

Evidence for the removal rules: phenix.refine's 2mFo-DFc < 1 sigma filter raised
R-free on the deposited waters of all five sets tried, and on our own waters in
five of six; removing 1K4C's 38 partnerless deposited waters raised it too.
Candidates are not gated on 2mFo-DFc (an unmodelled atom has only part of its
density there): with that gate ~2/3 of the missed deposited waters were refused.

Battery --only vs the Mg-fix binary (refined R-free, refinement time with the
build): 36gk -0.0022, 7l6j -0.0006, 7os3 0, 8agq 0, 8sa8 -0.0016, 6f3p -0.0008,
1hr2 -0.0107, 2r8s -0.0081, 193l -0.0049, 1exr -0.0019, 1k4c 0 (its round
rejected); no set worse. Time 2-4x the single refinement (1hr2 2.4 -> 7.2 s,
8sa8 11 -> 32 s); a rejected round costs a full refinement.
AlphaFold starts placed by MR (CLI, open-arm data): 9zlo 0.2940 -> 0.2564,
9gjx 0.3051 -> 0.2896, 9s02 0.2548 -> 0.2059.
Waters stripped from the refine-arm depositions and rebuilt: R-free recovers
75-100% of what the deposited waters are worth; recall/precision within 0.5 A
of a deposited water 0.25/0.25 (1hr2, 2.25 A RNA), 0.79/0.68 (2r8s), 0.63/0.57
(193l), 0.70/0.78 (1exr), 0.56/0.76 (1k4c).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37

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,

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
 `<prefix>_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 `<prefix>_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).

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

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 `<prefix>_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 `<prefix>_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

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
         $<$<BOOL:${JFJOCH_CUDA_AVAILABLE}>:RefinementGPU.cpp>
         $<$<BOOL:${JFJOCH_CUDA_AVAILABLE}>:RefinementGPU.cu>
         RefinementGPU.h RefinementGPUEngine.h

image_analysis/structure_refinement/Refinement.cpp

@@ -5,6 +5,7 @@
 
 #include <algorithm>
 #include <chrono>
+#include <set>
 #include <future>
 #include <map>
 #include <optional>
@@ -1225,7 +1226,10 @@ std::vector<TrialOutcome> RunTrials(const std::vector<TrialSlot> &trial_slots,
 
 }  // namespace
 
-RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<MergedReflection> &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<MergedReflection> &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<M
             const SolventMaps maps{result.map.get(), result.map_rms, result.difference_map.get(),
                                    result.difference_map_rms, result.anomalous_map.get(), result.anomalous_map_rms};
             result.solvent = AnalyseSolventSites(ctx.st, result.d_min, maps);
-            result.map.reset();
         }
+        result.b_ceiling = ctx.b_ceiling;
         result.k_sol = ctx.k_sol;
         result.b_sol = ctx.b_sol;
         result.model = std::make_shared<gemmi::Structure>(ctx.st);
@@ -1569,3 +1573,87 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
     RefinementStage::Log(logger);
     return result;
 }
+
+}  // namespace
+
+RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<MergedReflection> &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<std::string> added;   // the waters added and still there
+        std::vector<SolventEdit> 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<int>(added.size());
+        b.edits.insert(b.edits.end(), not_solvent.begin(), not_solvent.end());
+        b.seconds = std::chrono::duration<double>(std::chrono::steady_clock::now() - build_start).count();
+        result.seconds = std::chrono::duration<double>(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;
+}

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<double> 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<SolventSite> 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<MergedReflection> &merged,
                              const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double wavelength_A,

image_analysis/structure_refinement/SolventBuild.cpp

@@ -0,0 +1,464 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "SolventBuild.h"
+
+#include <algorithm>
+#include <cmath>
+#include <set>
+
+#include <spdlog/fmt/fmt.h>
+
+#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 <class F>
+void ForPointsNear(const gemmi::Grid<float> &g, const gemmi::Position &p, double radius, F &&f) {
+    const gemmi::Fractional fp = g.unit_cell.fractionalize(p);
+    const int u0 = static_cast<int>(std::round(fp.x * g.nu)), v0 = static_cast<int>(std::round(fp.y * g.nv)),
+              w0 = static_cast<int>(std::round(fp.z * g.nw));
+    const int du = static_cast<int>(std::ceil(radius / g.spacing[0])) + 1;
+    const int dv = static_cast<int>(std::ceil(radius / g.spacing[1])) + 1;
+    const int dw = static_cast<int>(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<float> &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<Peak> DifferencePeaks(const gemmi::Grid<float> &g, double rms, double cutoff) {
+    std::vector<Peak> out;
+    const std::vector<std::int8_t> 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<float> &g, const gemmi::Position &p, double cutoff) {
+    std::vector<std::pair<gemmi::Position, double>> 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<double> m{c[0] / sum, c[1] / sum, c[2] / sum, c[3] / sum, c[4] / sum, c[5] / sum};
+    std::array<double, 3> 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<const gemmi::Atom *> TruncatedEnds(const gemmi::Model &model) {
+    std::set<const gemmi::Atom *> 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<SolventSite> &sites,
+                         const std::set<std::string> &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<float> &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<const gemmi::Atom *> 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<double>(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<gemmi::Position> 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<gemmi::Position> 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<ModelNeighbour> near = AtomsNear(ns, model, peak.pos, PARTNER_MAX);
+        const ModelNeighbour *nearest = nullptr;
+        std::string clash, truncated_by;
+        std::vector<std::string> 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<SolventEdit> &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<int>(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<float>(e.b_iso);
+        r.atoms.push_back(a);
+        e.site = ResidueName(r, ch->name);
+        e.seqid = r.seqid;
+        ch->residues.push_back(r);
+    }
+}

image_analysis/structure_refinement/SolventBuild.h

@@ -0,0 +1,75 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <map>
+#include <set>
+#include <string>
+#include <vector>
+
+#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<SolventEdit> edits;
+    int added = 0, removed = 0, not_solvent = 0;
+    std::map<std::string, int> 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<SolventEdit> 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<SolventSite> &sites,
+                         const std::set<std::string> &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<SolventEdit> &edits);

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<SolventSite> &sites, const std::string &output_prefix) {
+    void AddSolventKeys(ReportSection &s, const std::vector<SolventSite> &sites, const SolventBuildResult &build,
+                        const std::string &output_prefix) {
         std::map<std::string, int> 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;

rugnux/Rugnux.h

@@ -419,6 +419,7 @@ struct RugnuxOutputs {
     std::shared_ptr<ModelMaps> refined_maps;                  // <prefix>_refined_maps.mtz
     std::optional<std::vector<DensityBlob>> blobs;            // <prefix>_blobs.json
     std::optional<std::vector<SolventSite>> solvent;          // <prefix>_solvent.json
+    std::optional<SolventBuildResult> solvent_build;         // and the water build, in it
     std::optional<MtzTable> unmerged;                         // <prefix>_unmerged.mtz
     std::optional<MtzTable> unmerged_partials;                // <prefix>_unmerged_partials.mtz
     std::string detector_jpeg;                                // <prefix>_detector.jpg; empty: none drawn

rugnux/RugnuxOutputFiles.cpp

@@ -26,8 +26,10 @@ void WriteDensityBlobs(const std::vector<DensityBlob> &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<SolventSite> &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<SolventSite> &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<SolventSite> &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");

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()));

rugnux/rugnux_cli.cpp

@@ -256,7 +256,7 @@ void print_usage() {
     std::cout << "   --scaling-iterations <num>      Cap on the per-frame scaling iterations; the loop stops earlier once the scales settle (default: 100)" << std::endl;
     std::cout << "   -z, --reference <file>          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 <label>      Reference column to use, by MTZ label; an mmCIF is read under the labels gemmi cif2mtz gives it (IMEAN, FP, I(+), FC, F-model, ...) (default: auto - F-model, else IMEAN/I, else FP/FOBS/F)" << std::endl;
-    std::cout << "   --model <file>                  After merging, validate vs this atomic model: R-free, sigma_A-weighted 2mFo-DFc/mFo-DFc maps, the model as placed (written beside them), and - where the merge kept the Bijvoet split - an anomalous difference map with the strongest anomalous scatterers named by the atom they sit on. PDB or mmCIF, gzipped or not, recognised by content rather than by extension. A prediction (AlphaFold or ModelCIF, pLDDT in the B column), or an experimental structure the rigid body does not place against its null, is placed by molecular replacement instead (see --mr). A model that is placed is then refined against the merge it was placed on - each chain as a rigid body and one B per residue, then every atom's position and B restrained by the dictionary carried in rugnux against a maximum-likelihood target, the X-ray weight chosen by R-free within geometry bounds - and the unexplained blobs of its difference map are listed: <prefix>_refined_model.cif, <prefix>_refined_maps.mtz, <prefix>_blobs.json, REFINE_ keys in the report. R_WORK, R_FREE and every other key stay those of the model as placed; the refinement decides nothing about the processing. It costs seconds on a GPU and typically 2-15 minutes on the CPU (no GPU, or --no-gpu); --no-refine skips it" << std::endl;
+    std::cout << "   --model <file>                  After merging, validate vs this atomic model: R-free, sigma_A-weighted 2mFo-DFc/mFo-DFc maps, the model as placed (written beside them), and - where the merge kept the Bijvoet split - an anomalous difference map with the strongest anomalous scatterers named by the atom they sit on. PDB or mmCIF, gzipped or not, recognised by content rather than by extension. A prediction (AlphaFold or ModelCIF, pLDDT in the B column), or an experimental structure the rigid body does not place against its null, is placed by molecular replacement instead (see --mr). A model that is placed is then refined against the merge it was placed on - each chain as a rigid body and one B per residue, then every atom's position and B restrained by the dictionary carried in rugnux against a maximum-likelihood target, the X-ray weight chosen by R-free within geometry bounds, then waters added at difference-map peaks with a hydrogen-bond partner and removed from negative difference density (to 2.8 A, up to three rounds, each kept where it lowers R-free) - and the unexplained blobs of its difference map are listed: <prefix>_refined_model.cif, <prefix>_refined_maps.mtz, <prefix>_blobs.json, REFINE_ keys in the report. R_WORK, R_FREE and every other key stay those of the model as placed; the refinement decides nothing about the processing. It costs seconds on a GPU and typically 2-15 minutes on the CPU (no GPU, or --no-gpu); --no-refine skips it" << std::endl;
     std::cout << "   --mr <file>                     As --model, but placed by molecular replacement straight away, without the rigid-body placement first: as many copies as the cell holds, in every space group the merge cannot tell apart and in the twin hypothesis' subgroup. A prediction is trimmed to the residues of pLDDT 70 and above, and where the AlphaFold DB's predicted-aligned-error JSON sits beside it under its own name it is cut into rigid domains; any other PDB or mmCIF is taken as an experimental structure, its polymer only (no waters, ligands, hydrogens or alternative conformers). Which of the two it is is read from the file. It may relabel the enantiomorph of the answer, and a different space group it chooses is written beside the answer as <prefix>_mr.mtz/.cif with the placed model <prefix>_mr_model.pdb. --mode mx only" << std::endl;
     std::cout << "   --no-refine                     Do not refine the placed model (for debugging): by default a model that --model or --mr places is refined (see --model)" << std::endl;
     std::cout << "   --export-unmerged               Write <prefix>_unmerged.mtz, an unmerged MTZ (POINTLESS column layout) of the integrated observations, for aimless / pointless / careless. On by default, whenever there is an output prefix. On a rotation run the partials of each reflection are summed into one full, written at the batch its rocking curve is centred on, with the summed rocking-curve fraction in FRACTIONCALC; an event that caught less of its rocking curve than --min-captured-fraction (or --min-partiality) is not written, exactly as in the merge, since a full is what the reading program is told these rows are. Intensities carry the deterministic per-reflection corrections and nothing else - Lorentz and polarization in the LP column, the sensor's angle-dependent efficiency in QE and the flight path's attenuation in FLIGHT, so raw counts are I / LP * QE * FLIGHT - while the partiality is not divided out and the per-image scale is not applied at all, since those programs scale the data themselves. Written in --mode mx and --mode scale, and with --no-merge" << std::endl;

tests/CMakeLists.txt

@@ -187,6 +187,7 @@ ADD_EXECUTABLE(jfjoch_test
         RefinementTest.cpp
         RefineOnlyTest.cpp
         SolventSitesTest.cpp
+        SolventBuildTest.cpp
         SearchSpaceGroupTest.cpp
         SearchSpaceGroupTwinTest.cpp
         SyntheticMergedReflections.h

tests/SolventBuildTest.cpp

@@ -0,0 +1,153 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+
+#include <cmath>
+#include <string>
+#include <vector>
+
+#include "../image_analysis/structure_refinement/SolventBuild.h"
+
+namespace {
+
+constexpr double EDGE = 24.0;   // A, a P1 box
+
+gemmi::Structure box() {
+    gemmi::Structure st;
+    st.cell = gemmi::UnitCell(EDGE, EDGE, EDGE, 90, 90, 90);
+    st.spacegroup_hm = "P 1";
+    st.setup_cell_images();
+    st.models.emplace_back(1);
+    st.models[0].chains.emplace_back("A");
+    return st;
+}
+
+void add_atom(gemmi::Structure &st, const std::string &resname, const std::string &atom, gemmi::El el,
+              const gemmi::Position &pos) {
+    gemmi::Chain &ch = st.models[0].chains[0];
+    gemmi::Residue r;
+    r.name = resname;
+    r.seqid = gemmi::SeqId(static_cast<int>(ch.residues.size()) + 1, ' ');
+    gemmi::Atom a;
+    a.name = atom;
+    a.element = gemmi::Element(el);
+    a.pos = pos;
+    a.occ = 1.0f;
+    a.b_iso = 20.0f;
+    r.atoms.push_back(a);
+    ch.residues.push_back(r);
+}
+
+// A map of the box, 0.3 A grid, in units of its rms (1), with a Gaussian peak of `height` at each of `peaks`.
+gemmi::Grid<float> map_with(const gemmi::UnitCell &cell, const std::vector<gemmi::Position> &peaks, float height,
+                            float base = 0.0f) {
+    gemmi::Grid<float> g;
+    g.unit_cell = cell;
+    g.spacegroup = gemmi::find_spacegroup_by_name("P 1");
+    g.set_size(80, 80, 80);
+    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 gemmi::Position p = g.get_position(u, v, w);
+                float value = base;
+                for (const gemmi::Position &q : peaks) {
+                    const double d = cell.find_nearest_image(p, q, gemmi::Asu::Any).dist();
+                    value += static_cast<float>(height * std::exp(-d * d / (2 * 0.35 * 0.35)));
+                }
+                g.data[g.index_q(u, v, w)] = value;
+            }
+    return g;
+}
+
+int count_kind(const SolventRound &r, SolventEdit::Kind kind) {
+    int n = 0;
+    for (const SolventEdit &e : r.edits)
+        n += e.kind == kind;
+    return n;
+}
+
+}  // namespace
+
+// A difference peak 2.8 A from a water becomes a water; one with only a carbon 3.3 A away (no hydrogen-bond partner)
+// does not, and is not reported either; the new water goes into the model as HOH.
+TEST_CASE("SolventBuild_peak_with_and_without_partner", "[Solvent]") {
+    gemmi::Structure st = box();
+    add_atom(st, "HOH", "O", gemmi::El::O, {6, 6, 6});
+    add_atom(st, "UNL", "C1", gemmi::El::C, {16, 16, 16});
+    const gemmi::Position with_partner(8.8, 6, 6), without(16, 16, 19.3);
+    const gemmi::Grid<float> diff = map_with(st.cell, {with_partner, without}, 6.0f);
+    const gemmi::Grid<float> map = map_with(st.cell, {}, 0.0f, 1.5f);
+    const SolventMaps maps{&map, 1.0, &diff, 1.0, nullptr, 0};
+
+    SolventRound r = PickSolvent(st, maps, {}, {}, 300, 1);
+    REQUIRE(r.added == 1);
+    CHECK(r.not_solvent == 0);
+    CHECK(r.refused["no hydrogen-bond partner"] == 1);
+    const SolventEdit &e = r.edits.at(0);
+    CHECK(e.kind == SolventEdit::Kind::Added);
+    CHECK(st.cell.find_nearest_image(e.pos, with_partner, gemmi::Asu::Any).dist() < 0.25);
+    CHECK(e.difference_sigma == Catch::Approx(6.0).margin(0.5));
+
+    ApplySolventEdits(st, r.edits);
+    const gemmi::Residue &added = st.models[0].chains[0].residues.back();
+    CHECK(added.name == "HOH");
+    CHECK(added.seqid.num == 3);
+    CHECK(added.atoms.at(0).b_iso == Catch::Approx(20.0));   // its partner's
+    CHECK(r.edits.at(0).site == "HOH A3");
+}
+
+// A sulfate-like cluster - a peak with four more at 1.5 A on a tetrahedron - beside a water is no water: the centre
+// is reported as a tetrahedral group, its corners as atoms of one molecule; nothing is added.
+TEST_CASE("SolventBuild_tetrahedral_cluster_refused", "[Solvent]") {
+    gemmi::Structure st = box();
+    add_atom(st, "HOH", "O", gemmi::El::O, {6, 6, 6});
+    const gemmi::Position centre(9.2, 6, 6);
+    const double t = 1.5 / std::sqrt(3.0);
+    std::vector<gemmi::Position> peaks = {centre};
+    for (const gemmi::Position &d : {gemmi::Position(t, t, t), gemmi::Position(t, -t, -t), gemmi::Position(-t, t, -t),
+                                     gemmi::Position(-t, -t, t)})
+        peaks.push_back(centre + d);
+    const gemmi::Grid<float> diff = map_with(st.cell, peaks, 6.0f);
+    const gemmi::Grid<float> map = map_with(st.cell, {}, 0.0f, 1.5f);
+    const SolventMaps maps{&map, 1.0, &diff, 1.0, nullptr, 0};
+
+    const SolventRound r = PickSolvent(st, maps, {}, {}, 300, 1);
+    CHECK(r.added == 0);
+    CHECK(count_kind(r, SolventEdit::Kind::NotSolvent) >= 1);
+    bool tetrahedral = false;
+    for (const SolventEdit &e : r.edits)
+        if (e.kind == SolventEdit::Kind::NotSolvent && st.cell.find_nearest_image(e.pos, centre, gemmi::Asu::Any).dist() < 0.3)
+            tetrahedral = e.reason.find("tetrahedral") != std::string::npos;
+    CHECK(tetrahedral);
+}
+
+// A water in negative difference density is taken out; one without a partner only where the build added it.
+TEST_CASE("SolventBuild_removal", "[Solvent]") {
+    gemmi::Structure st = box();
+    add_atom(st, "HOH", "O", gemmi::El::O, {6, 6, 6});
+    add_atom(st, "HOH", "O", gemmi::El::O, {8.8, 6, 6});
+    add_atom(st, "HOH", "O", gemmi::El::O, {16, 16, 16});
+    const gemmi::Grid<float> diff = map_with(st.cell, {}, 0.0f);
+    const gemmi::Grid<float> map = map_with(st.cell, {}, 0.0f, 1.5f);
+    const SolventMaps maps{&map, 1.0, &diff, 1.0, nullptr, 0};
+    std::vector<SolventSite> sites(3);
+    for (int k = 0; k < 3; k++) {
+        sites[k].label = "HOH";
+        sites[k].chain = "A";
+        sites[k].seqid = std::to_string(k + 1);
+        sites[k].pos = st.models[0].chains[0].residues[k].atoms[0].pos;
+        sites[k].map_sigma = 0.5;   // weak, which alone takes nothing out
+        sites[k].n_26_32 = k < 2;
+        sites[k].nearest_ion = INFINITY;
+        sites[k].b_iso = 20;
+    }
+    sites[1].difference_sigma = -4.0;
+
+    CHECK(PickSolvent(st, maps, sites, {}, 300, 1).removed == 1);
+    SolventRound r = PickSolvent(st, maps, sites, {"HOH A3"}, 300, 1);
+    CHECK(r.removed == 2);
+    ApplySolventEdits(st, r.edits);
+    REQUIRE(st.models[0].chains[0].residues.size() == 1);
+    CHECK(st.models[0].chains[0].residues[0].seqid.num == 1);
+}