RESULT-CHANGING: Refinement: the solvent build within the weight search, two rounds each judged by refining the solvent

Each round of the build was a refinement of its own - the whole weight search again, and a rejected round cost a
full refinement - so the build took 2-4 times the refinement. Now the atom list carries placeholder waters from the
start (occupancy 0: no density, no mask, no contact, not refined; a fifteenth of the heavy atoms, at least 64), and
the build runs where the weight trajectory ends: a round reads the maps and the audit of that model, puts waters into
placeholders, takes waters out (occupancy 0) and relabels sites in place (RetypeAtom: element, flags, metal bonds,
exclusions; the GPU engine refreshes an atom's density on an occupancy or element change), in every trial slot; then
the solvent alone is refined at the trajectory's weights (B refined, 10 iterations a cycle) and the round is kept
where that lowers R-free, else undone. At most two rounds; the final ladder then refines what they kept with the
rest. The output model is the view with the placeholders moved into their chains (SolventView) and the relabelled
metals' inner-sphere records. One topology for the whole refinement, so the metal records are decided once by
construction (the outer_records plumbing of the previous commit is gone again). Contacts() leaves out atoms at
occupancy 0; FindSpecialPositions() too (placeholders at the origin sat on the 2-folds of I222/P43212/C2 and were
projected there - R-free 0.44 on 36gk until fixed); waters put on a special position get their site constraint.
The audit runs on the threads (its density tables made up front), and MakeMaps makes its three maps at once.

Timing (--mode refine on the depositions, GPU, each set's three runs under one lock, load 15-24 from other agents),
refinement seconds without the build / with the old build / with this one, R-free old / new:
1hr2 3.5/11.3/4.6 0.2462/0.2490, 2r8s 2.8/11.2/4.6 0.2238/0.2244, 193l 1.0/3.1/1.8 0.1713/0.1713,
1exr 1.6/5.2/1.9 0.2076/0.2084, 1k4c 2.5/4.9/3.2 0.2106/0.2106, 36gk 3.7/12.4/6.7 0.2559/0.2543,
7l6j 3.6/16.3/6.5 0.1530/0.1546, 8agq 2.1/6.7/2.7 0.1753/0.1761, 6hwj 1.9/8.1/3.6 0.2222/0.2232,
8sa8 8.5/36.2/14.1 0.1714/0.1710, 6f3p 4.8/22.3/8.3 0.1629/0.1630, 7os3 1.4/4.9/2.6 0.2122/0.2115.
Build overhead 0.2-0.9x the refinement (median 0.7), was 1.0-3.7x (median 2.3); R-free within 0.001 on 8/12.
Tried and dropped (P3.md): the trajectory again after the rounds (never gained), a weight ladder after the build
(never displaced), building after the ladder instead (same), accepting rounds within the R-free tie (worse on 4),
a full step as the round's judge (noisier), one macro-cycle for the solvent (6f3p's rounds refused).

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

docs/CPU_DATA_ANALYSIS_REFINEMENT.md

@@ -367,12 +367,22 @@ a site denser than any of them (a cacodylate's As, a heavier metal) is called th
 Per site, everything above is in `<prefix>_solvent.json` (with `electrons`, `b_env`, `peak_shift_A`,
 `anomalous_expected`, `n7_ligands`, `n7_off_plane`); the counts are the `REFINE_SOLVENT_` keys of the report.
 
-**The solvent build.** After the refinement, and to 2.8 Å (phenix.refine's ordered-solvent limit), waters are
-built and removed and waters and ions relabelled between refinements, as phenix.refine builds waters in its
-macro-cycles: a round reads the refined model's maps and audit, changes its solvent, 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.
-`--refine-score-only` builds nothing, and `--no-refine-solvent` turns the build off (the audit is still made).
+**The solvent build.** Within the refinement, and to 2.8 Å (phenix.refine's ordered-solvent limit), waters are
+built and removed and waters and ions relabelled, as phenix.refine builds waters within its macro-cycles. Where the
+weight trajectory (§15.7) ends, a round reads the maps and the audit of that model and changes its solvent; then the
+solvent alone - its waters and single-atom ions, the rest held - is refined at the weights the trajectory ended at
+(its B too, at its balanced weight where the trajectory held B, and 10 iterations a macro-cycle), and the round is
+kept where that lowers R-free and undone otherwise. The rounds stop at one that changes nothing, at one undone, or
+after two (a third gained at most 0.0004 in R-free where it was kept, on twelve sets, and was refused on the others);
+the final ladder of §15.7 then refines what they kept with the rest. The atom list of a refinement is fixed - its
+engines are made for it - so the waters a round may add are in the model from the start as placeholders: one water
+each (as many as a fifteenth of the model's heavy atoms, at least 64) at occupancy 0, which puts no density, no mask
+and no contact anywhere and is not refined until a round puts one at a peak; a water taken out goes to occupancy 0
+the same way, and a relabelled site keeps its place. A round thus costs a map, an audit and a refinement of the
+solvent: on twelve refine-arm and open-arm depositions the build adds 0.2-0.9 times the refinement's time (median
+0.7), where a refinement of its own per round added 1-3.7 times (median 2.3), at an R-free within 0.001 of that on
+eight of them (two 0.0016-0.0028 higher, one 0.0016 lower). `--refine-score-only` builds nothing, and `--no-refine-solvent`
+turns the build off (the audit is still made).
 
 - *Relabelled* is a water or a single-atom ion (no alternative conformer) whose best-supported species is another
   one, not ambiguously - a water becomes an ion, an ion water or another ion - once per site. It takes the

docs/RUGNUX_REPORT.md

@@ -1210,8 +1210,8 @@ Where a model was placed and refined ([Refinement](RUGNUX_MODEL_AND_REFERENCE.md
 | `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_RELABELLED` | Waters and single-atom ions the build relabelled as another species (water, Mg, Na, K, Ca or Cl) where their surroundings and density support it, not ambiguously; each, with the scores, in the report's prose and in `<prefix>_solvent.json`. Not with `--no-refine-solvent` |
 | `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) |
+| `REFINE_SOLVENT_BUILD_ROUNDS` | The build's rounds kept (at most two, within the weight search), each kept where refining the solvent it changed lowered R-free |
+| `REFINE_SOLVENT_BUILD_SECONDS` | The build's time: its maps, audits and refinements of the solvent (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`,

image_analysis/structure_refinement/Refinement.cpp

@@ -54,6 +54,8 @@ constexpr double REPULSION_MARGIN = 0.5;   // A: pairs this much beyond their co
 // The geometry the weight is chosen within (Afonine et al. 2012 select among weights by R-free with the
 // geometry held to a target in the same way).
 constexpr double MAX_RMSD_BONDS = 0.02, MAX_RMSD_ANGLES = 2.0;
+// The minimiser's iterations per macro-cycle where the solvent build refines the solvent alone.
+constexpr int SOLVENT_ITERATIONS = 10;
 // R-free differences smaller than this do not decide between two weights: the tighter geometry is taken.
 constexpr double R_FREE_TIE = 0.001;
 // The minimiser (RefinementMinimiser): a solve stops where an iteration lowers the target by less than this,
@@ -925,11 +927,13 @@ void MakeMaps(Context &ctx, const std::vector<MergedReflection> &merged, double
         fwt_data.v.push_back({ctx.refl.hkl[m], std::polar(static_cast<float>(fwt), static_cast<float>(phi))});
     }
     result.maps = maps;
-    // The difference map beside the anomalous one below: neither reads the other.
-    std::future<void> difference_map = std::async(std::launch::async, [&diff, &fwt_data, &result] {
-        const std::array<int, 3> size = gemmi::get_size_for_hkl(diff, {{0, 0, 0}}, 3.0);
+    // The difference map and the 2mFo-DFc map beside the anomalous one below: none reads another.
+    const std::array<int, 3> size = gemmi::get_size_for_hkl(diff, {{0, 0, 0}}, 3.0);
+    std::future<void> difference_map = std::async(std::launch::async, [&diff, &size, &result] {
         result.difference_map = std::make_shared<gemmi::Grid<float>>(MapFromFPhi(gemmi::get_f_phi_on_grid<float>(diff, size, true)));
         result.difference_map_rms = gemmi::calculate_data_statistics(result.difference_map->data).rms;
+    });
+    std::future<void> map = std::async(std::launch::async, [&fwt_data, &size, &result] {
         result.map = std::make_shared<gemmi::Grid<float>>(MapFromFPhi(gemmi::get_f_phi_on_grid<float>(fwt_data, size, true)));
         result.map_rms = gemmi::calculate_data_statistics(result.map->data).rms;
     });
@@ -965,6 +969,7 @@ void MakeMaps(Context &ctx, const std::vector<MergedReflection> &merged, double
         }
     }
     difference_map.get();
+    map.get();
 }
 
 // The engine for the model as it is now: on the GPU where it is wanted and takes the case within its budget,
@@ -1052,8 +1057,7 @@ struct OwnSlot {
 // A copy of the model as it was before the topology was made (the same topology comes out of it, hydrogens
 // included), with an engine of the main one's kind - a GPU's on `device`. Null where the card cannot make one.
 std::unique_ptr<OwnSlot> MakeOwnSlot(const gemmi::Structure &before_topology, const Context &main, int device,
-                                     size_t nthreads, const RefinementSettings &settings,
-                                     const std::set<std::string> *outer_records) {
+                                     size_t nthreads, const RefinementSettings &settings) {
     auto slot = std::make_unique<OwnSlot>();
     Context &ctx = slot->ctx;
     ctx.st = before_topology;
@@ -1066,7 +1070,7 @@ std::unique_ptr<OwnSlot> MakeOwnSlot(const gemmi::Structure &before_topology, co
     ctx.cancel = main.cancel;
     ctx.b_ceiling = main.b_ceiling;
     ctx.mask = std::make_unique<BulkSolventMask>(ctx.cell, *ctx.sg, ctx.d_min, nthreads);
-    slot->topo = BuildRefinementTopology(ctx.st, settings.riding_hydrogens, outer_records);
+    slot->topo = BuildRefinementTopology(ctx.st, settings.riding_hydrogens);
     ctx.atoms = AtomList(ctx.Model());
     ctx.SetRiderParents(slot->topo.rider_parent);
 #ifdef JFJOCH_USE_CUDA
@@ -1262,12 +1266,110 @@ std::vector<TrialOutcome> RunTrials(const std::vector<TrialSlot> &trial_slots,
 
 namespace {
 
-// One refinement of a model whose atoms stay as they are (RefineModel adds and removes waters between calls).
-// `outer_records`: the metal records to leave as contacts, decided by an earlier round (null: decide here).
+// The solvent build inside the refinement (SolventBuild.h): the atom list of a refinement is fixed - its engines are
+// made for it - so the waters a round may add are there from the start, as placeholders: one water each in a chain
+// of their own at the end of the model, at occupancy 0, which puts no density, no mask and no contact anywhere and is
+// not refined, until a round puts one at a peak. A water taken out goes to occupancy 0 the same way, and a site
+// relabelled keeps its place in the list.
+const char *const PLACEHOLDER_CHAIN = "~";
+
+// As many placeholders as a fifteenth of the model's heavy atoms, at least 64: the depositions' builds add up to a
+// tenth (2R8S, 1HR2), most of it in the first round.
+int PlaceholderCount(const gemmi::Model &model) {
+    int heavy = 0;
+    for (const gemmi::Chain &ch : model.chains)
+        for (const gemmi::Residue &r : ch.residues)
+            for (const gemmi::Atom &a : r.atoms)
+                heavy += !a.is_hydrogen();
+    return std::max(64, heavy / 15);
+}
+
+void AddPlaceholders(gemmi::Structure &st, int n) {
+    gemmi::Chain ch(PLACEHOLDER_CHAIN);
+    for (int k = 0; k < n; k++) {
+        gemmi::Residue r;
+        r.name = "HOH";
+        r.seqid = gemmi::SeqId(k + 1, ' ');
+        r.het_flag = 'H';
+        r.entity_type = gemmi::EntityType::Water;
+        gemmi::Atom a;
+        a.name = "O";
+        a.element = gemmi::El::O;
+        a.occ = 0;
+        a.b_iso = 20;
+        r.atoms.push_back(a);
+        ch.residues.push_back(r);
+    }
+    st.models[0].chains.push_back(ch);
+}
+
+// The placeholders of a model (the last `n` atoms) and where those put in go: the chain of their nearest atom and the
+// number they take there.
+struct SolventPool {
+    int first = 0, n = 0;
+    std::vector<std::string> chain;   // empty: not put in
+    std::vector<gemmi::SeqId> seqid;
+    std::map<std::string, int> last;  // per chain name, the last residue number taken
+
+    int Free() const {
+        for (int k = 0; k < n; k++)
+            if (chain[k].empty())
+                return k;
+        return -1;
+    }
+};
+
+// The model as the solvent build and the output see it: the placeholders put in moved into their chains (the last of
+// the name) under their numbers, those not put in and the waters taken out (occupancy 0) left out. Each atom's serial
+// is its index in `st`.
+gemmi::Structure SolventView(const gemmi::Structure &st, const SolventPool &pool) {
+    gemmi::Structure view = st;
+    gemmi::Model &model = view.models[0];
+    int index = 0;
+    for (gemmi::Chain &ch : model.chains)
+        for (gemmi::Residue &r : ch.residues)
+            for (gemmi::Atom &a : r.atoms)
+                a.serial = index++;
+    std::vector<gemmi::Residue> moved;
+    std::vector<std::string> into;
+    for (gemmi::Chain &ch : model.chains)
+        if (ch.name == PLACEHOLDER_CHAIN)
+            for (size_t k = 0; k < ch.residues.size(); k++)
+                if (!pool.chain[k].empty() && ch.residues[k].atoms[0].occ > 0) {
+                    moved.push_back(ch.residues[k]);
+                    moved.back().seqid = pool.seqid[k];
+                    into.push_back(pool.chain[k]);
+                }
+    std::erase_if(model.chains, [](const gemmi::Chain &ch) { return ch.name == PLACEHOLDER_CHAIN; });
+    for (size_t k = 0; k < moved.size(); k++) {
+        gemmi::Chain *target = nullptr;
+        for (gemmi::Chain &ch : model.chains)
+            if (ch.name == into[k])
+                target = &ch;
+        if (target == nullptr) {
+            model.chains.emplace_back(into[k]);
+            target = &model.chains.back();
+        }
+        for (const gemmi::Residue &w : target->residues)
+            if (w.is_water()) {
+                moved[k].subchain = w.subchain;
+                moved[k].entity_id = w.entity_id;
+            }
+        target->residues.push_back(moved[k]);
+    }
+    for (gemmi::Chain &ch : model.chains) {
+        for (gemmi::Residue &r : ch.residues)
+            if (r.is_water())
+                std::erase_if(r.atoms, [](const gemmi::Atom &a) { return !(a.occ > 0); });
+        std::erase_if(ch.residues, [](const gemmi::Residue &r) { return r.atoms.empty(); });
+    }
+    return view;
+}
+
+// One refinement of a model, with the solvent build in it where settings.solvent asks for it.
 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 std::set<std::string> *outer_records = nullptr) {
+                             Logger &logger, const RefinementSettings &settings) {
     const auto start = std::chrono::steady_clock::now();
     RefinementResult result;
     Context ctx;
@@ -1306,6 +1408,21 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
 
     ctx.mask = std::make_unique<BulkSolventMask>(cell, sg, ctx.d_min, settings.nthreads);
 
+    // The solvent build's placeholders, before the slots copy the model.
+    const bool build = settings.solvent && settings.restrained && !settings.score_only &&
+                       result.d_min <= SOLVENT_BUILD_D_MIN;
+    SolventPool pool;
+    if (build) {
+        pool.n = PlaceholderCount(ctx.Model());
+        for (const gemmi::Chain &ch : ctx.Model().chains)
+            for (const gemmi::Residue &r : ch.residues)
+                if (r.seqid.num != gemmi::SeqId::OptionalNum::None)
+                    pool.last[ch.name] = std::max(pool.last[ch.name], static_cast<int>(r.seqid.num));
+        AddPlaceholders(ctx.st, pool.n);
+        pool.chain.assign(pool.n, "");
+        pool.seqid.assign(pool.n, gemmi::SeqId());
+    }
+
     // The dictionary's view of the model, and its hydrogens riding on their parents from the start, so that the
     // model is scored as the restrained stage refines it.
     const bool restrained = settings.restrained || settings.score_only;
@@ -1314,13 +1431,13 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
     RefinementTopology topo;
     if (restrained) {
         std::optional<RefinementStage> topo_stage("topology");
-        topo = BuildRefinementTopology(ctx.st, settings.riding_hydrogens, outer_records);
+        topo = BuildRefinementTopology(ctx.st, settings.riding_hydrogens);
         topo_stage.reset();
-        result.outer_records = topo.outer_records;
         ctx.atoms = AtomList(ctx.Model());
         ctx.SetRiderParents(topo.rider_parent);
         PlaceRiders(ctx, topo);
     }
+    pool.first = static_cast<int>(ctx.atoms.size()) - pool.n;   // the placeholders are the last atoms (no hydrogens)
 
     MakeEngine(ctx, settings.gpu, settings, logger);
     result.on_gpu = ctx.engine->OnGPU();
@@ -1337,12 +1454,32 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
                 result.on_gpu ? "GPU" : "CPU", r.r_work, r.r_free, ctx.k_sol, ctx.b_sol, result.score_start.b_iso,
                 result.score_start.mask_spacing);
 
+    // The waters and ions the build relabelled and kept (SolventEdit::site their old name), for the metal records the
+    // refined model is written with.
+    std::vector<SolventEdit> relabels;
     // The maps, the blobs and the model of where the refinement ended.
     auto finish = [&]() {
         MakeMaps(ctx, merged, wavelength_A, result);
+        // the model as the build left it: the waters put in in their chains, those taken out out, and each metal it
+        // relabelled with the records of its inner sphere in place of its old ones
+        gemmi::Structure model = build ? SolventView(ctx.st, pool) : ctx.st;
+        int record = 0;
+        for (const SolventEdit &e : relabels) {
+            const std::string from = e.site.substr(0, e.site.find(' '));
+            std::erase_if(model.connections, [&](const gemmi::Connection &c) {
+                auto is_site = [&](const gemmi::AtomAddress &a) {
+                    return a.chain_name == e.chain && a.res_id.seqid == e.seqid && a.res_id.name == from;
+                };
+                return is_site(c.partner1) || is_site(c.partner2);
+            });
+            for (const gemmi::Chain &ch : model.models[0].chains)
+                for (const gemmi::Residue &r : ch.residues)
+                    if (ch.name == e.chain && r.seqid == e.seqid && r.name == e.species)
+                        AddMetalRecords(model, ch, r, record);
+        }
         {
             const RefinementStage stage("blobs");
-            result.blobs = FindDensityBlobs(*result.difference_map, result.difference_map_rms, ctx.st, result.d_min,
+            result.blobs = FindDensityBlobs(*result.difference_map, result.difference_map_rms, model, result.d_min,
                                             result.anomalous_map.get(), result.anomalous_map_rms);
         }
         {
@@ -1350,12 +1487,17 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
             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,
                                    wavelength_A};
-            result.solvent = AnalyseSolventSites(ctx.st, result.d_min, maps);
+            result.solvent = AnalyseSolventSites(model, result.d_min, maps, ctx.nthreads);
         }
         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);
+        if (build)
+            for (gemmi::Chain &ch : model.models[0].chains)
+                for (gemmi::Residue &r : ch.residues)
+                    for (gemmi::Atom &a : r.atoms)
+                        a.serial = 0;   // SolventView's indices
+        result.model = std::make_shared<gemmi::Structure>(std::move(model));
         result.ok = true;
         result.seconds = std::chrono::duration<double>(std::chrono::steady_clock::now() - start).count();
     };
@@ -1410,7 +1552,7 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
 #endif
             // a slot that cannot be made is left out, as where the card cannot take its engine
             try {
-                own[k] = MakeOwnSlot(before_topology, ctx, devices[k], slot_threads, settings, outer_records);
+                own[k] = MakeOwnSlot(before_topology, ctx, devices[k], slot_threads, settings);
             } catch (const std::exception &e) {
                 logger.Warning("Refinement: a trial slot could not be made ({}), the trials share fewer", e.what());
             }
@@ -1427,7 +1569,7 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
 
     if (settings.restrained) {
         ctx.Refit();
-        result.restrained_residues = topo.restrained_residues;
+        result.restrained_residues = topo.restrained_residues - pool.n;   // not the placeholders
         result.fixed_residues = topo.fixed_residues;
         result.fixed_names = topo.fixed_names;
         std::vector<int> refinable;
@@ -1438,7 +1580,7 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
         logger.Info("Refinement: {} residues restrained, {} held fixed{}; {} bonds, {} metal bonds ({} links made, {} records "
                     "beyond the inner shell left as contacts), {} "
                     "angles, {} torsions, {} chiral volumes, {} planes",
-                    topo.restrained_residues, topo.fixed_residues,
+                    result.restrained_residues, topo.fixed_residues,
                     topo.fixed_names.empty() ? std::string() : " (" + Joined(topo.fixed_names) + ")",
                     topo.bonds, topo.metal_bonds.records.size(), topo.metal_links, topo.metal_links_refused, topo.angles, topo.torsions, topo.chirals, topo.planes);
 
@@ -1516,41 +1658,268 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
         //    A step is taken while it lowers R-free or R-work by more than the tie and its R-free stays within the tie
         //    of the best so far; the trajectory ends at its best step (on a deposited model, whose free set it was
         //    fitted to, R-work falls while R-free rises, and it ends where it began).
-        TrialOutcome best = state;
         // The trials of the step that went from the best so far: the final ladder goes from there too, and a trial
         // of the same weights from the same model is the same computation.
         std::vector<TrialOutcome> from_best;
-        bool state_is_best = true;
-        for (int step = 0; step < settings.max_steps && state.ran && !ctx.Cancelled(); step++) {
-            const double w = state.trial.weight, adp = state.trial.adp_weight;
-            std::vector<std::pair<double, double>> trials = {{w, adp}};
-            if (step % 2 == 1 && adp > 0) {
-                trials.emplace_back(w, adp / 4);
-                trials.emplace_back(w, 4 * adp);
-            } else {
-                trials.emplace_back(w / 2, adp);
-                if (2 * w <= w_max)
-                    trials.emplace_back(2 * w, adp);
+        // `fresh`: the first step fits the scale, the solvent and the likelihood anew (the model changed since).
+        auto trajectory = [&](TrialOutcome state, bool fresh, std::vector<TrialOutcome> &from_best) {
+            TrialOutcome best = state;
+            from_best.clear();
+            bool state_is_best = true;
+            for (int step = 0; step < settings.max_steps && state.ran && !ctx.Cancelled(); step++) {
+                const double w = state.trial.weight, adp = state.trial.adp_weight;
+                std::vector<std::pair<double, double>> trials = {{w, adp}};
+                if (step % 2 == 1 && adp > 0) {
+                    trials.emplace_back(w, adp / 4);
+                    trials.emplace_back(w, 4 * adp);
+                } else {
+                    trials.emplace_back(w / 2, adp);
+                    if (2 * w <= w_max)
+                        trials.emplace_back(2 * w, adp);
+                }
+                std::vector<TrialOutcome> outcomes = RunTrials(trial_slots, trials, refinable, state.pos, state.b,
+                                                               fresh && step == 0 ? nullptr : &state.fit,
+                                                               settings.macro_cycles, settings);
+                if (state_is_best)
+                    from_best = outcomes;
+                const int k = choose(outcomes, fmt::format("step {}", step + 1).c_str());
+                if (k < 0)
+                    break;
+                const RefinementWeightTrial &t = outcomes[k].trial;
+                // after a change of the model the first step is taken whatever it gains: it is where the change is
+                // judged from
+                const bool gains = (fresh && step == 0) || displaces(t, state.trial) ||
+                                   (t.within_bounds && t.r_work < state.trial.r_work - R_FREE_TIE);
+                if (!gains || !((fresh && step == 0) || (t.within_bounds && t.r_free < best.trial.r_free + R_FREE_TIE)))
+                    break;
+                state = std::move(outcomes[k]);
+                state_is_best = (fresh && step == 0) || displaces(state.trial, best.trial);
+                if (state_is_best) {
+                    best = state;
+                    from_best.clear();
+                }
             }
-            outcomes = RunTrials(trial_slots, trials, refinable, state.pos, state.b, &state.fit, settings.macro_cycles,
-                                 settings);
-            if (state_is_best)
-                from_best = outcomes;
-            k = choose(outcomes, fmt::format("step {}", step + 1).c_str());
-            if (k < 0)
-                break;
-            const RefinementWeightTrial &t = outcomes[k].trial;
-            const bool gains = displaces(t, state.trial) || (t.within_bounds && t.r_work < state.trial.r_work - R_FREE_TIE);
-            if (!gains || !(t.within_bounds && t.r_free < best.trial.r_free + R_FREE_TIE))
-                break;
-            state = std::move(outcomes[k]);
-            state_is_best = displaces(state.trial, best.trial);
-            if (state_is_best) {
-                best = state;
+            return best;
+        };
+        state = trajectory(std::move(state), false, from_best);
+
+        // 2a. The solvent build (SolventBuild.h), from where the trajectory ended: waters put in at the difference
+        //     map's peaks, taken out where the maps no longer hold them, and waters and ions relabelled, on the maps
+        //     and the solvent audit of that model; then the solvent alone refined, the rest held, at the weights the
+        //     trajectory ended at. A round is kept where that lowers R-free, and undone otherwise (the build ends
+        //     there); at most SOLVENT_BUILD_ROUNDS of them. What they kept is then refined with the rest by the final
+        //     ladder (the trajectory run again from there gained nothing on twelve sets). Each round thus costs a map,
+        //     an audit and a refinement of the solvent, where a refinement of its own cost the whole weight search
+        //     (the build was 2-4 times the refinement).
+        if (build && state.ran && !ctx.Cancelled()) {
+            const auto build_start = std::chrono::steady_clock::now();
+            const RefinementStage build_stage("solvent_build");
+            SolventBuildResult &b = result.solvent_build;
+            b.ran = true;
+            std::set<std::string> added;        // the waters put in and still there, by name
+            std::set<std::string> relabelled;   // the sites relabelled, by their new names
+            std::vector<SolventEdit> not_solvent;
+            // per slot, the residue of each atom (the atoms stay where they are in memory)
+            std::vector<std::vector<gemmi::Residue *>> residue_of(trial_slots.size());
+            for (size_t sl = 0; sl < trial_slots.size(); sl++)
+                for (gemmi::Chain &ch : trial_slots[sl].ctx->Model().chains)
+                    for (gemmi::Residue &r : ch.residues)
+                        for (size_t a = 0; a < r.atoms.size(); a++)
+                            residue_of[sl].push_back(&r);
+            for (int round = 1; !ctx.Cancelled(); round++) {
+                if (round > SOLVENT_BUILD_ROUNDS) {
+                    b.stop = fmt::format("{} rounds", SOLVENT_BUILD_ROUNDS);
+                    break;
+                }
+                // the maps and the audit of where the trajectory ended
+                ctx.SetPositions(state.pos);
+                ctx.SetB(state.b);
+                PlaceRiders(ctx, topo);
+                RefinementResult scratch;
+                MakeMaps(ctx, merged, wavelength_A, scratch);
+                std::optional<RefinementStage> pick_stage("build.view");
+                const gemmi::Structure view = SolventView(ctx.st, pool);
+                const SolventMaps maps{scratch.map.get(), scratch.map_rms, scratch.difference_map.get(),
+                                       scratch.difference_map_rms, scratch.anomalous_map.get(),
+                                       scratch.anomalous_map_rms, wavelength_A};
+                pick_stage.emplace("build.audit");
+                const std::vector<SolventSite> sites = AnalyseSolventSites(view, result.d_min, maps, settings.nthreads);
+                pick_stage.emplace("build.pick");
+                SolventRound picked = PickSolvent(view, maps, sites, added, relabelled, ctx.b_ceiling, round);
+                pick_stage.reset();
+                not_solvent.clear();
+                for (const SolventEdit &e : picked.edits)
+                    if (e.kind == SolventEdit::Kind::NotSolvent)
+                        not_solvent.push_back(e);
+                b.not_solvent = picked.not_solvent;
+                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, {} sites to relabel, {} peaks "
+                            "not solvent{}", round, picked.added, picked.removed, picked.relabelled,
+                            picked.not_solvent, refused.empty() ? std::string() : "; peaks refused: " + refused);
+                if (picked.added + picked.removed + picked.relabelled == 0) {
+                    b.stop = "nothing to change";
+                    break;
+                }
+                // the atoms the round edits (an index in the model of every slot), each with what it was
+                std::map<std::string, int> by_name;   // "HOH A401" (and "|altloc" where it has one) -> atom index
+                for (const gemmi::Chain &ch : view.models[0].chains)
+                    for (const gemmi::Residue &r : ch.residues)
+                        for (const gemmi::Atom &a : r.atoms)
+                            if (!a.is_hydrogen())
+                                by_name[r.name + " " + ch.name + r.seqid.str() +
+                                        (a.altloc ? std::string("|") + a.altloc : std::string())] = a.serial;
+                struct Undo {
+                    int index;
+                    gemmi::Atom atom;
+                    std::string resname;
+                    gemmi::EntityType entity;
+                };
+                std::vector<std::vector<Undo>> undo(trial_slots.size());
+                std::vector<RefinementTopology> topo_before;
+                for (const TrialSlot &slot : trial_slots)
+                    topo_before.push_back(*slot.topo);
+                const SolventPool pool_before = pool;
+                auto edit_atom = [&](int index, const std::function<void(gemmi::Atom &, gemmi::Residue &, size_t)> &f) {
+                    for (size_t sl = 0; sl < trial_slots.size(); sl++) {
+                        gemmi::Residue &r = *residue_of[sl][index];
+                        gemmi::Atom &a = *trial_slots[sl].ctx->atoms[index];
+                        undo[sl].push_back({index, a, r.name, r.entity_type});
+                        f(a, r, sl);
+                    }
+                };
+                std::vector<SolventEdit> kept;
+                for (SolventEdit &e : picked.edits) {
+                    if (e.kind == SolventEdit::Kind::Added) {
+                        const int k = pool.Free();
+                        if (k < 0)
+                            continue;   // no placeholder left: the weaker peaks wait
+                        pool.chain[k] = e.chain;
+                        pool.seqid[k] = gemmi::SeqId(++pool.last[e.chain], ' ');
+                        e.site = "HOH " + e.chain + pool.seqid[k].str();
+                        // on a special position (PickSolvent put it exactly there), held on it
+                        SpecialPosition site = SiteSymmetry(view.cell, e.pos);
+                        site.atom = pool.first + k;
+                        edit_atom(pool.first + k, [&](gemmi::Atom &a, gemmi::Residue &, size_t sl) {
+                            a.pos = e.pos;
+                            a.occ = static_cast<float>(e.occupancy);
+                            a.b_iso = static_cast<float>(e.b_iso);
+                            RefinementTopology &t = *trial_slots[sl].topo;
+                            t.refinable[pool.first + k] = 1;
+                            if (!site.ops.empty()) {
+                                t.special_index[pool.first + k] = static_cast<int>(t.special.size());
+                                t.special.push_back(site);
+                            }
+                        });
+                    } else if (e.kind == SolventEdit::Kind::Removed) {
+                        const auto it = by_name.find(e.site + (e.altloc ? std::string("|") + e.altloc : std::string()));
+                        if (it == by_name.end())
+                            continue;
+                        edit_atom(it->second, [&](gemmi::Atom &a, gemmi::Residue &, size_t sl) {
+                            a.occ = 0;
+                            trial_slots[sl].topo->refinable[it->second] = 0;
+                        });
+                    } else if (e.kind == SolventEdit::Kind::Relabelled) {
+                        const auto it = by_name.find(e.site);
+                        if (it == by_name.end())
+                            continue;
+                        const int i = it->second;
+                        // its inner sphere, for a metal's bonds, read in the view and numbered as the slots number it
+                        std::vector<int> donors;
+                        if (gemmi::Element(SolventSpeciesElement(e.species)).is_metal())
+                            for (const gemmi::Chain &ch : view.models[0].chains)
+                                for (const gemmi::Residue &r : ch.residues)
+                                    if (ch.name == e.chain && r.seqid == e.seqid && r.atoms.size() == 1) {
+                                        gemmi::Residue as_metal = r;
+                                        as_metal.atoms[0].element = gemmi::Element(SolventSpeciesElement(e.species));
+                                        for (const InnerSphereDonor &d : InnerSphereDonors(view.models[0], r, as_metal.atoms[0]))
+                                            donors.push_back(view.models[0].chains[d.chain].residues[d.residue].atoms[d.atom].serial);
+                                    }
+                        edit_atom(i, [&](gemmi::Atom &a, gemmi::Residue &r, size_t sl) {
+                            const bool water = e.species == "HOH";
+                            r.name = e.species;
+                            r.entity_type = water ? gemmi::EntityType::Water : gemmi::EntityType::NonPolymer;
+                            a.element = gemmi::Element(SolventSpeciesElement(e.species));
+                            a.name = water ? "O" : e.species;
+                            a.charge = ScatteringCharge(r, a);
+                            RetypeAtom(*trial_slots[sl].topo, trial_slots[sl].ctx->Model(), i, donors);
+                        });
+                    } else
+                        continue;
+                    kept.push_back(e);
+                }
+                refinable.clear();
+                for (size_t i = 0; i < ctx.atoms.size(); i++)
+                    if (topo.refinable[i])
+                        refinable.push_back(static_cast<int>(i));
+                // the solvent alone refined, the rest held, at the weights the trajectory ended at: what the round's
+                // changes are judged by
+                std::vector<int> solvent;
+                for (int i : refinable) {
+                    const gemmi::Residue &r = *residue_of[0][i];
+                    if (r.is_water() || (r.atoms.size() == 1 && (r.atoms[0].element.is_metal() || r.atoms[0].element == gemmi::El::Cl)))
+                        solvent.push_back(i);
+                }
+                const std::vector<gemmi::Position> changed_pos = ctx.Positions();
+                const std::vector<double> changed_b = ctx.B();
+                // (B refined, at its balanced weight where the trajectory held B: a new water's starting B is a
+                // guess; fewer iterations, the waters' few parameters converging fast)
+                RefinementSettings solvent_settings = settings;
+                solvent_settings.xyz_iterations = SOLVENT_ITERATIONS;
+                solvent_settings.b_iterations = SOLVENT_ITERATIONS;
+                std::vector<TrialOutcome> local =
+                        RunTrials(trial_slots, {{state.trial.weight, state.trial.adp_weight > 0 ? state.trial.adp_weight : adp_w0}},
+                                  solvent, changed_pos, changed_b, nullptr, settings.macro_cycles, solvent_settings);
+                logger.Info("Solvent build, round {}: the solvent refined alone: {}", round, local[0].line);
+                TrialOutcome after = std::move(local[0]);
+                if (!after.ran || !(after.trial.r_free < state.trial.r_free)) {
+                    logger.Info("Solvent build, round {}: not kept (R-free {:.4f}, before {:.4f})", round,
+                                after.trial.r_free, state.trial.r_free);
+                    b.stop = fmt::format("round {} raised R-free", round);
+                    // undone, in every slot
+                    for (size_t sl = 0; sl < trial_slots.size(); sl++) {
+                        for (auto u = undo[sl].rbegin(); u != undo[sl].rend(); ++u) {
+                            *trial_slots[sl].ctx->atoms[u->index] = u->atom;
+                            residue_of[sl][u->index]->name = u->resname;
+                            residue_of[sl][u->index]->entity_type = u->entity;
+                        }
+                        *trial_slots[sl].topo = topo_before[sl];
+                    }
+                    pool = pool_before;
+                    refinable.clear();
+                    for (size_t i = 0; i < ctx.atoms.size(); i++)
+                        if (topo.refinable[i])
+                            refinable.push_back(static_cast<int>(i));
+                    break;
+                }
+                for (const SolventEdit &e : kept) {
+                    if (e.kind == SolventEdit::Kind::Added)
+                        added.insert(e.site);
+                    else if (e.kind == SolventEdit::Kind::Removed) {
+                        if (added.erase(e.site) == 0)
+                            b.removed++;
+                    } else {
+                        const std::string now = e.species + e.site.substr(e.site.find(' '));
+                        if (added.erase(e.site) > 0)
+                            added.insert(now);
+                        relabelled.insert(now);
+                        relabels.push_back(e);
+                        b.relabelled++;
+                    }
+                    b.edits.push_back(e);
+                }
+                b.rounds = round;
+                state = std::move(after);
                 from_best.clear();
             }
+            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();
+            logger.Info("Solvent build: {} waters added, {} removed, {} sites relabelled in {} round{} kept ({}), {} "
+                        "peaks not solvent; R-free {:.4f}; {:.1f} s", b.added, b.removed, b.relabelled, b.rounds,
+                        b.rounds == 1 ? "" : "s", b.stop, b.not_solvent, state.trial.r_free, b.seconds);
         }
-        state = std::move(best);
 
         // 3. The final weights: a short ladder from where the trajectory ended, tighter geometry first, its choice
         //    kept where it displaces that. A rung the step from there tried already is taken from that step.
@@ -1617,92 +1986,7 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
 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;
-        // the metal records left as contacts on the model as given stay so in every round
-        const std::set<std::string> outer_records = result.outer_records;
-        std::set<std::string> added;   // the waters added and still there
-        std::set<std::string> relabelled;   // the sites relabelled, by their new names
-        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, relabelled, 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, {} sites to relabel, {} peaks not "
-                        "solvent{}", round, picked.added, picked.removed, picked.relabelled, picked.not_solvent,
-                        refused.empty() ? std::string() : "; peaks refused: " + refused);
-            b.not_solvent = picked.not_solvent;
-            if (picked.added + picked.removed + picked.relabelled == 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, &outer_records);
-            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);
-                } else if (e.kind == SolventEdit::Kind::Relabelled) {
-                    const std::string now = e.species + e.site.substr(e.site.find(' '));
-                    if (added.erase(e.site) > 0)
-                        added.insert(now);
-                    relabelled.insert(now);
-                    b.relabelled++;
-                    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, {} sites relabelled in {} round{} kept ({}), {} peaks "
-                    "not solvent; R-free {:.4f}; {:.1f} s", b.added, b.removed, b.relabelled, 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

@@ -60,9 +60,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; off (rugnux --no-refine-solvent) the solvent stays as given and is
-    // only audited.
+    // Waters added at the difference map's peaks and removed where the maps no longer hold them, and waters and ions
+    // relabelled, within the weight trajectory (SolventBuild.h), to SOLVENT_BUILD_D_MIN; off (rugnux
+    // --no-refine-solvent) the solvent stays as given and is only audited.
     bool solvent = true;
     bool gpu = true;
     size_t nthreads = 1;
@@ -166,14 +166,13 @@ struct RefinementResult {
     // audit, nothing in the model is changed by it.
     std::vector<SolventSite> solvent;
     double b_ceiling = 0;                      // the highest B the refinement allowed
-    std::set<std::string> outer_records;       // the metal records left as contacts (RefinementTopology)
-    SolventBuildResult solvent_build;          // the waters added and removed (settings.solvent)
+    SolventBuildResult solvent_build;          // the solvent build's rounds (settings.solvent)
 };
 
 // `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.
+// the weight trajectory is followed by up to SOLVENT_BUILD_ROUNDS rounds of the solvent build, each the trajectory
+// again from the model with its solvent changed, kept where it ends at a lower R-free (Refinement.cpp). `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,
                              Logger &logger, const RefinementSettings &settings);

image_analysis/structure_refinement/RefinementGPU.cpp

@@ -303,6 +303,7 @@ private:
         const int n = static_cast<int>(atoms.size());
         if (last_b_.size() != atoms.size()) {
             last_b_.assign(atoms.size(), NAN);
+            last_key_.assign(atoms.size(), AtomKey());
             atoms_.resize(atoms.size());
             pos_.resize(atoms.size());
         }
@@ -310,9 +311,12 @@ private:
         const int n3[3] = {setup_.grid.nu, setup_.grid.nv, setup_.grid.nw};
         ParallelChunks(n, ThreadsForWork(atoms.size(), nthreads_, 4096), [&](int lo, int hi) {
             for (int i = lo; i < hi; i++) {
-                if (!(atoms[i]->b_iso == last_b_[i])) {
+                // made again where its B moved, or where the solvent build put it in, took it out or relabelled it
+                const AtomKey key{atoms[i]->occ, static_cast<int>(atoms[i]->element.elem), atoms[i]->charge};
+                if (!(atoms[i]->b_iso == last_b_[i]) || !(key == last_key_[i])) {
                     atoms_[i] = DensityAtom(*atoms[i], atoms[i]->b_iso, dc_);
                     last_b_[i] = atoms[i]->b_iso;
+                    last_key_[i] = key;
                     changed_[i] = 1;
                 }
                 const gemmi::Fractional f = cell_.fractionalize(atoms[i]->pos);
@@ -338,6 +342,13 @@ private:
     std::unique_ptr<RefinementGPUEngine> engine_;
     std::vector<const gemmi::Atom *> atom_ptrs_;   // the model's atoms at the last SetAtoms()
     std::vector<float> last_b_;
+    struct AtomKey {
+        float occ = -1;
+        int element = -1;
+        signed char charge = 0;
+        bool operator==(const AtomKey &o) const { return occ == o.occ && element == o.element && charge == o.charge; }
+    };
+    std::vector<AtomKey> last_key_;
     PinnedVector<ModelDensityAtom> atoms_;   // pinned: a B rung sends them every evaluation (37 MB on a large model)
     PinnedVector<std::array<float, 4>> pos_;
     std::vector<char> changed_;

image_analysis/structure_refinement/RefinementRestraints.cpp

@@ -141,13 +141,10 @@ std::string MetalDonorModification(const std::string &resname, const std::string
 // 2.9 A) would otherwise pull the ion onto that O; it stays a contact (BuildRepulsions, METAL_CONTACT), and the
 // coordination counts only the records taken.
 // Records to a symmetry mate are left out, as every bond to a symmetry mate is (BuildRefinementTopology).
-// With `outer_given`, the records left as contacts are those it names instead, whatever the distances now (a
-// refinement's later rounds keep the first round's decision). Returns the number of links made; `outer` the
-// records left as contacts.
+// Returns the number of links made; `refused` the records left as contacts for their distance.
 // Distances as surveyed by Harding (2001) Acta Cryst. D57, 401-411 and (2006) Acta Cryst. D62, 678-682,
 // and for Mg, Na and K at nucleic acids by Leonarski, D'Ascenzo & Auffinger (2017) Nucleic Acids Res. 45, 987-1004
-int AddMetalLinks(gemmi::Structure &st, gemmi::MonLib &monlib, const std::set<std::string> *outer_given,
-                  std::set<std::string> &outer) {
+int AddMetalLinks(gemmi::Structure &st, gemmi::MonLib &monlib, int &refused) {
     gemmi::Model &model = st.models.at(0);
     struct Site {
         gemmi::Connection *conn;
@@ -175,7 +172,7 @@ int AddMetalLinks(gemmi::Structure &st, gemmi::MonLib &monlib, const std::set<st
     // as they are): the site's coordination is then every record taken, a symmetry mate's or a described one too.
     std::vector<Site> sites;
     std::unordered_map<const gemmi::Atom *, int> coordination;
-    outer.clear();
+    refused = 0;
     for (const Site &rec : records) {
         const gemmi::Connection &conn = *rec.conn;
         const gemmi::CRA &a = rec.metal_first ? rec.metal : rec.donor, &b = rec.metal_first ? rec.donor : rec.metal;
@@ -185,17 +182,13 @@ int AddMetalLinks(gemmi::Structure &st, gemmi::MonLib &monlib, const std::set<st
                                                              b.atom->name, conn.partner2.altloc)) != nullptr;
         double distance = 0, sigma = 0;
         const int claimed_n = DictionaryMetalBonds(monlib, rec.metal) + claimed[rec.metal.atom];
-        const std::string key = MetalRecordKey(conn);
-        const bool beyond = outer_given != nullptr
-                            ? outer_given->count(key) > 0
-                            : !described && EmbeddedMetalDistance(rec.metal.atom->element.elem,
-                                                                  rec.donor.atom->element.elem, claimed_n, distance,
-                                                                  sigma) && rec.d > distance + INNER_SHELL_MARGIN;
-        if (beyond) {
+        if (!described &&
+            EmbeddedMetalDistance(rec.metal.atom->element.elem, rec.donor.atom->element.elem, claimed_n, distance, sigma) &&
+            rec.d > distance + INNER_SHELL_MARGIN) {
             // named so that prepare_topology() matches it to no link, not even the one made for another site
             // of the same residues and atoms (BuildRefinementTopology restores the record's own name after it)
             rec.conn->link_id = OUTER_SPHERE_RECORD;
-            outer.insert(key);
+            refused++;
             continue;
         }
         if (coordination.count(rec.metal.atom) == 0)
@@ -230,8 +223,66 @@ int AddMetalLinks(gemmi::Structure &st, gemmi::MonLib &monlib, const std::set<st
 
 }  // namespace
 
-std::string MetalRecordKey(const gemmi::Connection &conn) {
-    return conn.partner1.str() + " - " + conn.partner2.str();
+void RetypeAtom(RefinementTopology &topo, const gemmi::Model &model, int i, const std::vector<int> &donors) {
+    const gemmi::Residue *res = nullptr;
+    const gemmi::Atom *atom = nullptr;
+    int k = 0;
+    for (const gemmi::Chain &ch : model.chains)
+        for (const gemmi::Residue &r : ch.residues)
+            for (const gemmi::Atom &a : r.atoms)
+                if (k++ == i) {
+                    res = &r;
+                    atom = &a;
+                }
+    if (res == nullptr)
+        return;
+    const gemmi::Atom &a = *atom;
+    const gemmi::MonLib &monlib = *topo.monlib;
+    topo.element[i] = a.element.elem;
+    topo.metal[i] = IsMetal(a.element);
+    topo.vdw_radius[i] = VdwRadius(monlib, res->name, a);
+    topo.hbond[i] = res->is_water() || CanHydrogenBond(monlib, res->name, a);
+    // its bonds go, with what they excluded from the repulsions
+    for (int j : topo.excluded[i])
+        if (j != i)
+            std::erase(topo.excluded[j], i);
+    topo.excluded[i] = {i};
+    auto without_i = [&](const RestraintSet &set) {
+        RestraintSet out;
+        for (RestraintRecord r : set.records)
+            if (r.atom[0] != i && r.atom[1] != i) {
+                r.slot = 0;
+                out.Add(r, 2);
+            }
+        return out;
+    };
+    topo.metal_bonds = without_i(topo.metal_bonds);
+    RestraintSet adp;
+    for (const RestraintRecord &r : topo.adp.records)
+        if (r.atom[0] != i && r.atom[1] != i)
+            adp.Add(r, 2);
+    topo.adp = adp;
+    // a metal: a bond to each of its donors, at the coordination they make
+    if (topo.metal[i]) {
+        for (int j : donors) {
+            double distance = 0, sigma = 0;
+            if (!EmbeddedMetalDistance(a.element.elem, topo.element[j], static_cast<int>(donors.size()), distance, sigma))
+                continue;
+            RestraintRecord r;
+            r.kind = static_cast<int>(RestraintKind::Bond);
+            r.atom[0] = j;
+            r.atom[1] = i;
+            r.ideal = static_cast<float>(distance);
+            r.sigma = static_cast<float>(sigma);
+            topo.metal_bonds.Add(r, 2);
+            topo.excluded[i].push_back(j);
+            topo.excluded[j].push_back(i);
+            std::sort(topo.excluded[j].begin(), topo.excluded[j].end());
+        }
+        std::sort(topo.excluded[i].begin(), topo.excluded[i].end());
+    }
+    topo.metal_bonds.Finalize(topo.element.size());
+    topo.adp.Finalize(topo.element.size());
 }
 
 void RestraintSet::Add(RestraintRecord r, int n_atoms_of_r) {
@@ -435,14 +486,13 @@ std::vector<SpecialPosition> FindSpecialPositions(const gemmi::Structure &st) {
             for (const gemmi::Atom &a : r.atoms) {
                 SpecialPosition sp = SiteSymmetry(cell, a.pos);
                 sp.atom = i++;
-                if (!sp.ops.empty())
+                if (!sp.ops.empty() && a.occ > 0)   // an atom that is not there (occupancy 0) is on no site
                     out.push_back(sp);
             }
     return out;
 }
 
-RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens,
-                                           const std::set<std::string> *outer_records) {
+RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens) {
     RefinementTopology out;
     gemmi::Model &model = st.models.at(0);
     out.monlib = std::make_shared<gemmi::MonLib>(EmbeddedMonomerLibrary());
@@ -453,8 +503,7 @@ RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens,
     std::vector<std::string> link_ids;
     for (const gemmi::Connection &conn : st.connections)
         link_ids.push_back(conn.link_id);
-    out.metal_links = AddMetalLinks(st, monlib, outer_records, out.outer_records);
-    out.metal_links_refused = static_cast<int>(out.outer_records.size());
+    out.metal_links = AddMetalLinks(st, monlib, out.metal_links_refused);
     gemmi::Logger quiet;   // missing atoms and unknown residues are what the fixed residues report
     quiet.threshold = 0;
     out.topo = gemmi::prepare_topology(st, monlib, 0,
@@ -718,11 +767,14 @@ RestraintSet Contacts(const RefinementTopology &topo, gemmi::Model &model, const
 
     std::vector<std::vector<int>> first(model.chains.size());
     std::vector<gemmi::Position> pos;
+    std::vector<char> absent;   // occupancy 0: an atom that is not there makes no contact
     for (size_t c = 0; c < model.chains.size(); c++)
         for (const gemmi::Residue &r : model.chains[c].residues) {
             first[c].push_back(static_cast<int>(pos.size()));
-            for (const gemmi::Atom &a : r.atoms)
+            for (const gemmi::Atom &a : r.atoms) {
                 pos.push_back(a.pos);
+                absent.push_back(!(a.occ > 0));
+            }
         }
 
     // Per atom on a special position, the symmetry images of it the search lists: of the copies that coincide (the
@@ -758,12 +810,12 @@ RestraintSet Contacts(const RefinementTopology &topo, gemmi::Model &model, const
     ParallelFor(BLOCKS, nthreads, [&](int block) {
         const size_t lo = pos.size() * block / BLOCKS, hi = pos.size() * (block + 1) / BLOCKS;
         for (size_t i = lo; i < hi; i++) {
-            if (topo.rider_parent[i] >= 0 || topo.element[i] == gemmi::El::H)
+            if (topo.rider_parent[i] >= 0 || topo.element[i] == gemmi::El::H || absent[i])
                 continue;
             const gemmi::Fractional fi = cell.fractionalize(pos[i]);
             ns.for_each(pos[i], '\0', reach, [&](gemmi::NeighborSearch::Mark &m, double) {
                 const int j = first[m.chain_idx][m.residue_idx] + m.atom_idx;
-                if (topo.rider_parent[j] >= 0 || topo.element[j] == gemmi::El::H)
+                if (topo.rider_parent[j] >= 0 || topo.element[j] == gemmi::El::H || absent[j])
                     return;
                 if (!topo.refinable[i] && !topo.refinable[j])
                     return;

image_analysis/structure_refinement/RefinementRestraints.h

@@ -102,7 +102,6 @@ struct RefinementTopology {
     int bonds = 0, angles = 0, torsions = 0, chirals = 0, planes = 0;
     int metal_links = 0;                   // metal-donor bonds made from the model's LINK records (AddMetalLinks)
     int metal_links_refused = 0;           // LINK records at an outer-sphere distance, left as contacts
-    std::set<std::string> outer_records;   // and which (MetalRecordKey)
     // The atoms on special positions (FindSpecialPositions), placed exactly on them; per atom its entry, -1 none.
     std::vector<SpecialPosition> special;
     std::vector<int> special_index;
@@ -119,17 +118,19 @@ struct RefinementTopology {
 // With `hydrogens`, the dictionary's hydrogens are added to the residues it describes (waters excepted) and
 // placed; they ride on their parent atoms and take part in no restraint.
 // The atoms on special positions are moved onto them (SpecialPosition::OnSite).
-// `outer_records`: where given, the metal records left as contacts are these (MetalRecordKey) whatever the model's
-// distances now - a decision taken once, on the model as given, and not again on distances its refinement moved.
-RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens = false,
-                                           const std::set<std::string> *outer_records = nullptr);
+RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens = false);
 
 // A metal's inner sphere ends this far beyond the library's distance: Mg-O 2.07 A + 0.33 = 2.4 A, the inner-sphere
 // bound of Leonarski, D'Ascenzo & Auffinger (2017) Nucleic Acids Res. 45, 987-1004, taken for every metal.
 constexpr double INNER_SHELL_MARGIN = 0.33;           // A
 
-// The name a metal LINK / struct_conn record is known by: its two partners.
-std::string MetalRecordKey(const gemmi::Connection &conn);
+// Atom `i` of `model` (in model order) retyped where its residue was relabelled (a water become an ion, an ion a
+// water): its element, van der Waals radius, hydrogen-bond and metal flags taken again, its bonds and B similarities
+// removed with the repulsions they excluded, and - a metal - a bond made to each of `donors` (atom indices), at the
+// library's distance for the coordination they make. The atoms stay as they are numbered, so an engine made for them
+// goes on.
+void RetypeAtom(RefinementTopology &topo, const gemmi::Model &model, int i, const std::vector<int> &donors);
+
 
 // The riding hydrogens placed again on their parents as they are now (gemmi's place_hydrogens_on_all_atoms).
 void PlaceRidingHydrogens(RefinementTopology &topo);

image_analysis/structure_refinement/SolventBuild.cpp

@@ -451,12 +451,68 @@ SolventRound PickSolvent(const gemmi::Structure &st, const SolventMaps &maps, co
     return out;
 }
 
-namespace {
+std::vector<InnerSphereDonor> InnerSphereDonors(const gemmi::Model &model, const gemmi::Residue &own, const gemmi::Atom &metal) {
+    std::vector<InnerSphereDonor> donors;
+    double widest = 0, sigma = 0;
+    for (gemmi::El el : {gemmi::El::O, gemmi::El::N, gemmi::El::S}) {
+        double d0 = 0;
+        if (EmbeddedMetalDistance(metal.element.elem, el, 6, d0, sigma))
+            widest = std::max(widest, d0 + INNER_SHELL_MARGIN + 0.2);
+    }
+    for (size_t c = 0; c < model.chains.size(); c++)
+        for (size_t r = 0; r < model.chains[c].residues.size(); r++) {
+            const gemmi::Residue &res = model.chains[c].residues[r];
+            for (size_t k = 0; k < res.atoms.size(); k++) {
+                const gemmi::Atom &d = res.atoms[k];
+                if (&res != &own && d.occ > 0 &&
+                    (d.element == gemmi::El::O || d.element == gemmi::El::N || d.element == gemmi::El::S))
+                    if (const double dist = d.pos.dist(metal.pos); dist < widest)
+                        donors.push_back({static_cast<int>(c), static_cast<int>(r), static_cast<int>(k), dist});
+            }
+        }
+    auto ideal = [&](const InnerSphereDonor &d, int cn, double &value) {
+        const gemmi::Atom &a = model.chains[d.chain].residues[d.residue].atoms[d.atom];
+        return EmbeddedMetalDistance(metal.element.elem, a.element.elem, cn, value, sigma);
+    };
+    // judged at the coordination of those within the inner sphere of an octahedral site, as AddMetalLinks judges
+    // a model's records at the coordination they claim
+    int cn = 0;
+    for (const InnerSphereDonor &d : donors) {
+        double value = 0;
+        cn += ideal(d, 6, value) && d.distance <= value + INNER_SHELL_MARGIN;
+    }
+    std::erase_if(donors, [&](const InnerSphereDonor &d) {
+        double value = 0;
+        return !ideal(d, std::max(cn, 1), value) || d.distance > value + INNER_SHELL_MARGIN;
+    });
+    return donors;
+}
 
-gemmi::El SpeciesElement(const std::string &species) {
+void AddMetalRecords(gemmi::Structure &st, const gemmi::Chain &chain, const gemmi::Residue &res, int &record) {
+    const gemmi::Model &model = st.models.at(0);
+    const gemmi::Atom &a = res.atoms.at(0);
+    if (!a.element.is_metal())
+        return;
+    for (const InnerSphereDonor &d : InnerSphereDonors(model, res, a)) {
+        const gemmi::Residue &dr = model.chains[d.chain].residues[d.residue];
+        const gemmi::Atom &da = dr.atoms[d.atom];
+        gemmi::Connection c;
+        c.name = "rugnux_metalc" + std::to_string(++record);
+        c.type = gemmi::Connection::MetalC;
+        c.asu = gemmi::Asu::Same;
+        c.partner1 = gemmi::AtomAddress(model.chains[d.chain].name, dr.seqid, dr.name, da.name, da.altloc);
+        c.partner2 = gemmi::AtomAddress(chain.name, res.seqid, res.name, a.name);
+        c.reported_distance = d.distance;
+        st.connections.push_back(c);
+    }
+}
+
+gemmi::El SolventSpeciesElement(const std::string &species) {
     return species == "HOH" ? gemmi::El::O : gemmi::Element(species).elem;
 }
 
+namespace {
+
 // The site of `e` relabelled as e.species, its metal records replaced by its inner sphere's.
 void Relabel(gemmi::Structure &st, const SolventEdit &e, int &record) {
     gemmi::Model &model = st.models.at(0);
@@ -488,53 +544,9 @@ void Relabel(gemmi::Structure &st, const SolventEdit &e, int &record) {
     res->subchain = like != nullptr ? like->subchain : "";
     res->entity_id = like != nullptr ? like->entity_id : "";
     gemmi::Atom &a = res->atoms[0];
-    a.element = gemmi::Element(SpeciesElement(e.species));
+    a.element = gemmi::Element(SolventSpeciesElement(e.species));
     a.name = water ? "O" : e.species;
-    a.occ = std::max(a.occ, 0.0f);
-    if (!a.element.is_metal())
-        return;
-    // the N, O and S of its inner sphere, within the model's own copy (a record to a symmetry mate is left out)
-    struct Donor {
-        const gemmi::Chain *ch;
-        const gemmi::Residue *r;
-        const gemmi::Atom *atom;
-        double d;
-    };
-    std::vector<Donor> donors;
-    double widest = 0, sigma = 0;
-    for (gemmi::El el : {gemmi::El::O, gemmi::El::N, gemmi::El::S}) {
-        double d0 = 0;
-        if (EmbeddedMetalDistance(a.element.elem, el, 6, d0, sigma))
-            widest = std::max(widest, d0 + INNER_SHELL_MARGIN + 0.2);
-    }
-    for (const gemmi::Chain &ch : model.chains)
-        for (const gemmi::Residue &r : ch.residues)
-            for (const gemmi::Atom &d : r.atoms)
-                if (&r != res && (d.element == gemmi::El::O || d.element == gemmi::El::N || d.element == gemmi::El::S))
-                    if (const double dist = d.pos.dist(a.pos); dist < widest)
-                        donors.push_back({&ch, &r, &d, dist});
-    // judged at the coordination of those within the inner sphere of an octahedral site, as AddMetalLinks judges
-    // a model's records at the coordination they claim
-    int cn = 0;
-    for (const Donor &d : donors) {
-        double ideal = 0;
-        cn += EmbeddedMetalDistance(a.element.elem, d.atom->element.elem, 6, ideal, sigma) &&
-              d.d <= ideal + INNER_SHELL_MARGIN;
-    }
-    for (const Donor &d : donors) {
-        double ideal = 0;
-        if (!EmbeddedMetalDistance(a.element.elem, d.atom->element.elem, std::max(cn, 1), ideal, sigma) ||
-            d.d > ideal + INNER_SHELL_MARGIN)
-            continue;
-        gemmi::Connection c;
-        c.name = "rugnux_metalc" + std::to_string(++record);
-        c.type = gemmi::Connection::MetalC;
-        c.asu = gemmi::Asu::Same;
-        c.partner1 = gemmi::AtomAddress(d.ch->name, d.r->seqid, d.r->name, d.atom->name, d.atom->altloc);
-        c.partner2 = gemmi::AtomAddress(chain->name, res->seqid, res->name, a.name);
-        c.reported_distance = d.d;
-        st.connections.push_back(c);
-    }
+    AddMetalRecords(st, *chain, *res, record);
 }
 
 }  // namespace

image_analysis/structure_refinement/SolventBuild.h

@@ -12,7 +12,7 @@
 
 #include "SolventSites.h"   // SolventSite, SolventMaps
 
-// One round of the solvent build that RefineModel runs between refinements (Refinement.h): waters placed at the
+// One round of the solvent build that RefineModel runs within its weight trajectory (Refinement.h): waters placed at the
 // peaks of the difference map a refined model leaves, the model's waters the maps no longer support taken out, and
 // waters and ions relabelled.
 //
@@ -54,8 +54,9 @@ struct SolventRound {
 // (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;
+// At most this many rounds, each judged by refining the solvent it changed (a third one gained at most 0.0004 in
+// R-free where it was kept, on twelve sets, and was refused on the others).
+constexpr int SOLVENT_BUILD_ROUNDS = 2;
 
 // What the build did, over the rounds kept.
 struct SolventBuildResult {
@@ -78,6 +79,23 @@ SolventRound PickSolvent(const gemmi::Structure &model, const SolventMaps &maps,
                          const std::set<std::string> &built, const std::set<std::string> &relabelled, double b_ceiling,
                          int round);
 
+// The N, O and S of a metal's inner sphere in the model's own copy (no symmetry mate): within the library's distance
+// plus INNER_SHELL_MARGIN, at the coordination of those within an octahedral site's - the atoms a metal's records
+// name. `own` is the metal's residue.
+struct InnerSphereDonor {
+    int chain = 0, residue = 0, atom = 0;   // indices into the model
+    double distance = 0;
+};
+std::vector<InnerSphereDonor> InnerSphereDonors(const gemmi::Model &model, const gemmi::Residue &own,
+                                                const gemmi::Atom &metal);
+
+// The element of a species by its residue name (HOH: O, MG: Mg, ...).
+gemmi::El SolventSpeciesElement(const std::string &species);
+
+// A record (struct_conn metalc, named rugnux_metalcN with N on from `record`) for each of InnerSphereDonors() of the
+// single-atom ion `res` of `chain`, where it is a metal.
+void AddMetalRecords(gemmi::Structure &st, const gemmi::Chain &chain, const gemmi::Residue &res, int &record);
+
 // 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). A relabelled site gets its new residue name,
 // element and atom name; its metal records go, and a metal gets one for each N, O or S of its inner sphere (within

image_analysis/structure_refinement/SolventSites.cpp

@@ -6,6 +6,7 @@
 #include <algorithm>
 #include <cmath>
 #include <map>
+#include <mutex>
 #include <set>
 #include <tuple>
 
@@ -16,6 +17,7 @@
 #include "gemmi/neighbor.hpp"
 #include "gemmi/resinfo.hpp"
 
+#include "../../common/ParallelFor.h"
 #include "ModelGrid.h"             // IonFormFactor
 #include "MonomerLibrary.h"
 #include "RefinementRestraints.h"   // AtomsNear, SiteSymmetry
@@ -145,16 +147,25 @@ class RadialDensity {
 public:
     explicit RadialDensity(double d_min) : d_min_(d_min) {}
     double operator()(const Scatterer &sc, double b, double r) { return At(Table(sc, b), r); }
+    // The tables of these scatterers at these B made up front, on `nthreads` threads: Table() then only reads them.
+    void Prepare(const std::vector<std::pair<Scatterer, double>> &wanted, size_t nthreads) {
+        std::map<Key, std::vector<double>> made;
+        for (const auto &[sc, b] : wanted)
+            made[KeyOf(sc, b)];
+        std::vector<std::pair<const Key, std::vector<double>> *> todo;
+        for (auto &entry : made)
+            if (!tables_.count(entry.first))
+                todo.push_back(&entry);
+        ParallelFor(static_cast<int>(todo.size()), nthreads, [&](int k) { todo[k]->second = Make(todo[k]->first); });
+        for (auto *entry : todo)
+            tables_.emplace(entry->first, std::move(entry->second));
+    }
     const std::vector<double> &Table(const Scatterer &sc, double b) {
-        const int bi = static_cast<int>(std::lround(std::clamp(b, 0.0, 500.0)));
-        const auto key = std::make_tuple(static_cast<int>(sc.el), static_cast<int>(sc.charge), sc.hydrogens, bi);
+        const std::lock_guard lock(m_);   // the sites are analysed on several threads; a table, once made, stays
+        const Key key = KeyOf(sc, b);
         auto it = tables_.find(key);
-        if (it == tables_.end()) {
-            std::vector<double> t(static_cast<size_t>(TAIL_RADIUS / STEP) + 2);
-            for (size_t k = 0; k < t.size(); k++)
-                t[k] = ScattererHeight(sc, bi, d_min_, k * STEP);
-            it = tables_.emplace(key, std::move(t)).first;
-        }
+        if (it == tables_.end())
+            it = tables_.emplace(key, Make(key)).first;
         return it->second;
     }
     static double At(const std::vector<double> &t, double r) {
@@ -166,9 +177,23 @@ public:
     }
 
 private:
+    using Key = std::tuple<int, int, int, int>;   // element, charge, hydrogens, B to the nearest A^2
+    static Key KeyOf(const Scatterer &sc, double b) {
+        return {static_cast<int>(sc.el), static_cast<int>(sc.charge), sc.hydrogens,
+                static_cast<int>(std::lround(std::clamp(b, 0.0, 500.0)))};
+    }
+    std::vector<double> Make(const Key &key) const {
+        const Scatterer sc{static_cast<gemmi::El>(std::get<0>(key)), static_cast<signed char>(std::get<1>(key)),
+                           std::get<2>(key)};
+        std::vector<double> t(static_cast<size_t>(TAIL_RADIUS / STEP) + 2);
+        for (size_t k = 0; k < t.size(); k++)
+            t[k] = ScattererHeight(sc, std::get<3>(key), d_min_, k * STEP);
+        return t;
+    }
     static constexpr double STEP = 0.05;
     double d_min_;
-    std::map<std::tuple<int, int, int, int>, std::vector<double>> tables_;
+    std::mutex m_;
+    std::map<Key, std::vector<double>> tables_;
 };
 
 // An atom near a site, as it scatters (for the density it puts there).
@@ -494,7 +519,8 @@ bool SolventLabelAgrees(const SolventSite &s) {
     return LabelSpecies(s.label) == static_cast<int>(s.best);
 }
 
-std::vector<SolventSite> AnalyseSolventSites(const gemmi::Structure &st, double d_min, const SolventMaps &maps) {
+std::vector<SolventSite> AnalyseSolventSites(const gemmi::Structure &st, double d_min, const SolventMaps &maps,
+                                             size_t nthreads) {
     std::vector<SolventSite> out;
     if (st.models.empty())
         return out;
@@ -523,6 +549,14 @@ std::vector<SolventSite> AnalyseSolventSites(const gemmi::Structure &st, double
     };
     const bool density = maps.map != nullptr && d_min > 0;
     RadialDensity radial(d_min);
+    if (density) {
+        std::vector<std::pair<Scatterer, double>> wanted;
+        for (const gemmi::Chain &ch : model.chains)
+            for (const gemmi::Residue &r : ch.residues)
+                for (const gemmi::Atom &a : r.atoms)
+                    wanted.emplace_back(ScattererOf(r, a), a.b_iso);
+        radial.Prepare(wanted, nthreads);
+    }
     // The model's atoms near the one at (ci, ri, ai) at `own` - in `near`, AtomsNear's list - other than it and its
     // coinciding copies, as they scatter.
     auto others_of = [&](size_t ci, size_t ri, size_t ai, const gemmi::Position &own,
@@ -552,7 +586,7 @@ std::vector<SolventSite> AnalyseSolventSites(const gemmi::Structure &st, double
     std::map<const gemmi::Atom *, double> atom_calibration;
     std::pair<double, double> global_calibration{NAN, NAN};
     if (density) {
-        std::vector<double> c;
+        std::vector<std::array<size_t, 3>> atoms;
         for (size_t ci = 0; ci < model.chains.size(); ci++)
             for (size_t ri = 0; ri < model.chains[ci].residues.size(); ri++) {
                 const gemmi::Residue &r = model.chains[ci].residues[ri];
@@ -560,18 +594,27 @@ std::vector<SolventSite> AnalyseSolventSites(const gemmi::Structure &st, double
                     continue;
                 for (size_t ai = 0; ai < r.atoms.size(); ai++) {
                     const gemmi::Atom &a = r.atoms[ai];
-                    if (!(a.element == gemmi::El::C || Polar(a.element.elem)) || a.occ != 1 || a.altloc)
-                        continue;
-                    const double model_height =
-                            radial(ScattererOf(r, a), a.b_iso, 0.0) +
-                            tails(others_of(ci, ri, ai, a.pos, AtomsNear(ns, model, a.pos, TAIL_RADIUS)), a.pos);
-                    const double v = maps.map->interpolate_value(a.pos) / model_height;
-                    if (v > 0) {
-                        atom_calibration[&a] = std::log(v);
-                        c.push_back(std::log(v));
-                    }
+                    if ((a.element == gemmi::El::C || Polar(a.element.elem)) && a.occ == 1 && !a.altloc)
+                        atoms.push_back({ci, ri, ai});
                 }
             }
+        std::vector<double> v(atoms.size());
+        ParallelFor(static_cast<int>(atoms.size()), nthreads, [&](int k) {
+            const auto [ci, ri, ai] = atoms[k];
+            const gemmi::Residue &r = model.chains[ci].residues[ri];
+            const gemmi::Atom &a = r.atoms[ai];
+            const double model_height =
+                    radial(ScattererOf(r, a), a.b_iso, 0.0) +
+                    tails(others_of(ci, ri, ai, a.pos, AtomsNear(ns, model, a.pos, TAIL_RADIUS)), a.pos);
+            v[k] = maps.map->interpolate_value(a.pos) / model_height;
+        });
+        std::vector<double> c;
+        for (size_t k = 0; k < atoms.size(); k++)
+            if (v[k] > 0) {
+                const auto [ci, ri, ai] = atoms[k];
+                atom_calibration[&model.chains[ci].residues[ri].atoms[ai]] = std::log(v[k]);
+                c.push_back(std::log(v[k]));
+            }
         global_calibration = MedianSpread(c);
     }
     // f'' by Cromer & Liberman (1981) Acta Cryst. A37, 267-268, corrected as Kissel & Pratt (1990) Acta Cryst. A46,
@@ -603,181 +646,189 @@ std::vector<SolventSite> AnalyseSolventSites(const gemmi::Structure &st, double
             anomalous_scale = num / den;
     }
 
+    // the sites, in model order, each measured on a thread of its own
+    std::vector<std::array<size_t, 3>> where;
     for (size_t ci = 0; ci < model.chains.size(); ci++)
         for (size_t ri = 0; ri < model.chains[ci].residues.size(); ri++) {
             const gemmi::Residue &res = model.chains[ci].residues[ri];
-            if (!IsWater(res) && IonAtom(res) == nullptr)
+            if (IsWater(res) || IonAtom(res) != nullptr)
+                for (size_t ai = 0; ai < res.atoms.size(); ai++)
+                    if (!res.atoms[ai].is_hydrogen())
+                        where.push_back({ci, ri, ai});
+        }
+    out.resize(where.size());
+    site_ligands.resize(where.size());
+    site_expected.resize(where.size());
+    ParallelFor(static_cast<int>(where.size()), nthreads, [&](int index) {
+        const auto [ci, ri, ai] = where[index];
+        const gemmi::Residue &res = model.chains[ci].residues[ri];
+        const gemmi::Atom &a = res.atoms[ai];
+        SolventSite &site = out[index];
+        site.label = res.name;
+        site.chain = model.chains[ci].name;
+        site.seqid = res.seqid.str();
+        site.altloc = a.altloc;
+        site.pos = a.pos;
+        site.occupancy = a.occ;
+        site.b_iso = a.b_iso;
+        site.sigma_xyz = CoordinateError(a.b_iso, d_min);
+        site.min_ligand_distance = NAN;
+        site.nearest_positive = site.nearest_negative = site.nearest_ion = INFINITY;
+        site.map_sigma = map_sigma(maps.map, maps.map_rms, a.pos);
+        site.difference_sigma = map_sigma(maps.difference, maps.difference_rms, a.pos);
+        site.anomalous_sigma = map_sigma(maps.anomalous, maps.anomalous_rms, a.pos);
+        const std::vector<ModelNeighbour> near = AtomsNear(ns, model, a.pos, SEARCH_RADIUS + SOLVENT_PEAK_REACH);
+        // The environment's B, and the map's scale from the macromolecule's atoms nearby where there are
+        // enough of them.
+        double b_sum = 0;
+        int b_n = 0;
+        std::vector<double> calibration;
+        for (const ModelNeighbour &nb : near) {
+            const gemmi::Atom &pa = model.chains[nb.chain].residues[nb.residue].atoms[nb.atom];
+            if (nb.distance > SEARCH_RADIUS || (nb.chain == static_cast<int>(ci) && nb.residue == static_cast<int>(ri)))
                 continue;
-            for (size_t ai = 0; ai < res.atoms.size(); ai++) {
-                const gemmi::Atom &a = res.atoms[ai];
-                if (a.is_hydrogen())
-                    continue;
-                SolventSite site;
-                site.label = res.name;
-                site.chain = model.chains[ci].name;
-                site.seqid = res.seqid.str();
-                site.altloc = a.altloc;
-                site.pos = a.pos;
-                site.occupancy = a.occ;
-                site.b_iso = a.b_iso;
-                site.sigma_xyz = CoordinateError(a.b_iso, d_min);
-                site.min_ligand_distance = NAN;
-                site.nearest_positive = site.nearest_negative = site.nearest_ion = INFINITY;
-                site.map_sigma = map_sigma(maps.map, maps.map_rms, a.pos);
-                site.difference_sigma = map_sigma(maps.difference, maps.difference_rms, a.pos);
-                site.anomalous_sigma = map_sigma(maps.anomalous, maps.anomalous_rms, a.pos);
-                const std::vector<ModelNeighbour> near = AtomsNear(ns, model, a.pos, SEARCH_RADIUS + SOLVENT_PEAK_REACH);
-                // The environment's B, and the map's scale from the macromolecule's atoms nearby where there are
-                // enough of them.
-                double b_sum = 0;
-                int b_n = 0;
-                std::vector<double> calibration;
-                for (const ModelNeighbour &nb : near) {
-                    const gemmi::Atom &pa = model.chains[nb.chain].residues[nb.residue].atoms[nb.atom];
-                    if (nb.distance > SEARCH_RADIUS || (nb.chain == static_cast<int>(ci) && nb.residue == static_cast<int>(ri)))
-                        continue;
-                    if (const auto it = atom_calibration.find(&pa); it != atom_calibration.end())
-                        calibration.push_back(it->second);
-                    if (Polar(pa.element.elem) && nb.distance < BVS_RADIUS && (!a.altloc || !pa.altloc || a.altloc == pa.altloc)) {
-                        b_sum += pa.b_iso;
-                        b_n++;
-                    }
-                }
-                site.b_env = b_n > 0 ? b_sum / b_n : NAN;
-                const double b_env = b_n > 0 ? site.b_env : a.b_iso;
-                const double copies = a.occ * SiteSymmetry(cell, a.pos).multiplicity;
-                const auto [log_k, spread] = static_cast<int>(calibration.size()) >= CALIBRATION_NEIGHBOURS
-                                             ? MedianSpread(calibration) : global_calibration;
-                // The site's own density fitted with its neighbours' held, at the environment's B; an ion, whose
-                // label's restraints may hold it off its density, where that fit puts it.
-                const bool ion = !IsWater(res);
-                gemmi::Position q = a.pos;
-                site.electrons = site.electron_spread = NAN;
-                if (density && std::isfinite(log_k) && copies > 0) {
-                    const LocalFit fit = FitSite(*maps.map, std::exp(log_k), a.pos, b_env, others_of(ci, ri, ai, a.pos, near),
-                                                 radial, ion);
-                    q = fit.p;
-                    site.electrons = 10.0 * fit.a / copies;
-                    site.electron_spread = std::max(ELECTRON_SPREAD_MIN, spread);
-                }
-                site.measured_at = q;
-                site.peak_shift = q.dist(a.pos);
+            if (const auto it = atom_calibration.find(&pa); it != atom_calibration.end())
+                calibration.push_back(it->second);
+            if (Polar(pa.element.elem) && nb.distance < BVS_RADIUS && (!a.altloc || !pa.altloc || a.altloc == pa.altloc)) {
+                b_sum += pa.b_iso;
+                b_n++;
+            }
+        }
+        site.b_env = b_n > 0 ? b_sum / b_n : NAN;
+        const double b_env = b_n > 0 ? site.b_env : a.b_iso;
+        const double copies = a.occ * SiteSymmetry(cell, a.pos).multiplicity;
+        const auto [log_k, spread] = static_cast<int>(calibration.size()) >= CALIBRATION_NEIGHBOURS
+                                     ? MedianSpread(calibration) : global_calibration;
+        // The site's own density fitted with its neighbours' held, at the environment's B; an ion, whose
+        // label's restraints may hold it off its density, where that fit puts it.
+        const bool ion = !IsWater(res);
+        gemmi::Position q = a.pos;
+        site.electrons = site.electron_spread = NAN;
+        if (density && std::isfinite(log_k) && copies > 0) {
+            const LocalFit fit = FitSite(*maps.map, std::exp(log_k), a.pos, b_env, others_of(ci, ri, ai, a.pos, near),
+                                         radial, ion);
+            q = fit.p;
+            site.electrons = 10.0 * fit.a / copies;
+            site.electron_spread = std::max(ELECTRON_SPREAD_MIN, spread);
+        }
+        site.measured_at = q;
+        site.peak_shift = q.dist(a.pos);
 
-                std::vector<Ligand> ligands;
-                std::vector<gemmi::Position> bound_cations;   // other cations closer than 2.6 A
-                std::vector<gemmi::Position> inner;   // the N and O closer than 2.6 A
-                gemmi::Vec3 vec[BVS_CATIONS];
-                for (ModelNeighbour nb : near) {
-                    // the site itself, its own hydrogens and its other conformers; its symmetry mates stay
-                    if (nb.chain == static_cast<int>(ci) && nb.residue == static_cast<int>(ri) &&
-                        (nb.atom != static_cast<int>(ai) || nb.pos.dist(a.pos) < 0.5))
-                        continue;
-                    nb.distance = q.dist(nb.pos);
-                    if (!(nb.distance < SEARCH_RADIUS))
-                        continue;
-                    const gemmi::Chain &pch = model.chains[nb.chain];
-                    const gemmi::Residue &pr = pch.residues[nb.residue];
-                    const gemmi::Atom &pa = pr.atoms[nb.atom];
-                    if (pa.is_hydrogen() || (a.altloc && pa.altloc && a.altloc != pa.altloc))
-                        continue;
-                    const double d = nb.distance;
-                    const int charge = IonCharge(pr);
-                    if (charge != 0) {
-                        site.ions_within_6++;
-                        site.ions_within_42 += d < 4.2;
-                        site.nearest_ion = std::min(site.nearest_ion, d);
-                        if (charge > 0) {
-                            site.nearest_positive = std::min(site.nearest_positive, d);
-                            site.positive_contacts += d < LIKE_CHARGE_DISTANCE;
-                            if (d < 2.6)
-                                bound_cations.push_back(nb.pos);
-                        } else {
-                            site.nearest_negative = std::min(site.nearest_negative, d);
-                            site.negative_contacts += d < HALIDE_NEGATIVE_DISTANCE;
-                        }
-                        continue;
-                    }
-                    const AtomChemistry chem = ChemistryOf(monlib, pr, pa, is_n_terminus(pr));
-                    if (chem.positive) {
-                        site.nearest_positive = std::min(site.nearest_positive, d);
-                        site.positive_contacts += d < LIKE_CHARGE_DISTANCE;
-                    }
-                    if (chem.negative) {
-                        site.nearest_negative = std::min(site.nearest_negative, d);
-                        site.negative_contacts += d < HALIDE_NEGATIVE_DISTANCE;
-                    }
-                    if (!Polar(pa.element.elem))
-                        continue;
-
-                    if (d < 1.9) site.n_below_19++;
-                    else if (d < 2.3) site.n_19_23++;
-                    else if (d < 2.6) site.n_23_26++;
-                    else if (d < 3.2) site.n_26_32++;
-                    else if (d < 3.4) site.n_32_34++;
-                    else if (d < 3.6) site.n_34_36++;
-                    if (d < SHELL_RADIUS) {
-                        if (chem.donor && chem.acceptor)
-                            site.donor_acceptors++;
-                        else if (chem.donor)
-                            site.donors++;
-                        else if (chem.acceptor)
-                            site.acceptors++;
-                        site.neighbours += fmt::format("{}{} {:.2f}", site.neighbours.empty() ? "" : ", ",
-                                                       AtomName(pr, pa, pch.name), d);
-                    }
+        std::vector<Ligand> ligands;
+        std::vector<gemmi::Position> bound_cations;   // other cations closer than 2.6 A
+        std::vector<gemmi::Position> inner;   // the N and O closer than 2.6 A
+        gemmi::Vec3 vec[BVS_CATIONS];
+        for (ModelNeighbour nb : near) {
+            // the site itself, its own hydrogens and its other conformers; its symmetry mates stay
+            if (nb.chain == static_cast<int>(ci) && nb.residue == static_cast<int>(ri) &&
+                (nb.atom != static_cast<int>(ai) || nb.pos.dist(a.pos) < 0.5))
+                continue;
+            nb.distance = q.dist(nb.pos);
+            if (!(nb.distance < SEARCH_RADIUS))
+                continue;
+            const gemmi::Chain &pch = model.chains[nb.chain];
+            const gemmi::Residue &pr = pch.residues[nb.residue];
+            const gemmi::Atom &pa = pr.atoms[nb.atom];
+            if (pa.is_hydrogen() || (a.altloc && pa.altloc && a.altloc != pa.altloc))
+                continue;
+            const double d = nb.distance;
+            const int charge = IonCharge(pr);
+            if (charge != 0) {
+                site.ions_within_6++;
+                site.ions_within_42 += d < 4.2;
+                site.nearest_ion = std::min(site.nearest_ion, d);
+                if (charge > 0) {
+                    site.nearest_positive = std::min(site.nearest_positive, d);
+                    site.positive_contacts += d < LIKE_CHARGE_DISTANCE;
                     if (d < 2.6)
-                        inner.push_back(nb.pos);
-
-                    Ligand l{};
-                    l.pos = nb.pos;
-                    l.distance = d;
-                    const double e = CoordinateError(pa.b_iso, d_min);
-                    l.sigma = std::sqrt(site.sigma_xyz * site.sigma_xyz + e * e);
-                    l.chem = chem;
-                    l.nitrogen = pa.element == gemmi::El::N;
-                    l.n7_off_plane = false;
-                    if (pa.name == "N7" && d < 2.6) {
-                        // the purine's ring, as the image the neighbour is, and its O6 or N6
-                        std::vector<gemmi::Atom> ring;
-                        const gemmi::Atom *o6n6 = nullptr;
-                        for (const gemmi::Atom &ra : pr.atoms) {
-                            static const std::set<std::string> RING = {"N9", "C8", "N7", "C5", "C6", "N1", "C2", "N3", "C4"};
-                            if (RING.count(ra.name) && (!a.altloc || !ra.altloc || ra.altloc == a.altloc)) {
-                                ring.push_back(ra);
-                                ring.back().pos = gemmi::Position(nb.image.apply(ra.pos));
-                            }
-                            if (ra.name == "O6" || ra.name == "N6")
-                                o6n6 = &ra;
-                        }
-                        if (ring.size() == 9 && o6n6 != nullptr) {
-                            site.n7_ligands++;
-                            std::vector<gemmi::Atom *> ptrs;
-                            for (gemmi::Atom &ra : ring)
-                                ptrs.push_back(&ra);
-                            const double plane = std::fabs(gemmi::get_distance_from_plane(q, gemmi::find_best_plane(ptrs)));
-                            const double side = q.dist(gemmi::Position(nb.image.apply(o6n6->pos)));
-                            l.n7_off_plane = plane > N7_PLANE_MAX + l.sigma || side < N7_O6N6_MIN - l.sigma ||
-                                             side > N7_O6N6_MAX + l.sigma;
-                            site.n7_off_plane += l.n7_off_plane;
-                        }
-                    }
-                    if (d < BVS_RADIUS)
-                        for (int c = 0; c < BVS_CATIONS; c++) {
-                            // Following Brown & Altermatt (1985) Acta Cryst. B41, 244-247, occupancy weighted as
-                            // Muller et al. (2003) Acta Cryst. D59, 32-37
-                            l.valence[c] = pa.occ * std::exp(((l.nitrogen ? R0_N[c] : R0_O[c]) - d) / BV_B);
-                            site.bvs[c] += l.valence[c];
-                            vec[c] += (nb.pos - q) * (l.valence[c] / d);
-                        }
-                    ligands.push_back(l);
+                        bound_cations.push_back(nb.pos);
+                } else {
+                    site.nearest_negative = std::min(site.nearest_negative, d);
+                    site.negative_contacts += d < HALIDE_NEGATIVE_DISTANCE;
                 }
-                // A site bound to a cation shares that cation's other ligands as neighbours: the cis waters of a
-                // Mg(H2O)6 are 2.9 A from each other, and are no hydrogen bonds of theirs. The scores leave them out.
-                std::erase_if(ligands, [&](const Ligand &l) {
-                    for (const gemmi::Position &c : bound_cations)
-                        if (l.pos.dist(c) < 2.6)
-                            return true;
-                    return false;
+                continue;
+            }
+            const AtomChemistry chem = ChemistryOf(monlib, pr, pa, is_n_terminus(pr));
+            if (chem.positive) {
+                site.nearest_positive = std::min(site.nearest_positive, d);
+                site.positive_contacts += d < LIKE_CHARGE_DISTANCE;
+            }
+            if (chem.negative) {
+                site.nearest_negative = std::min(site.nearest_negative, d);
+                site.negative_contacts += d < HALIDE_NEGATIVE_DISTANCE;
+            }
+            if (!Polar(pa.element.elem))
+                continue;
+
+            if (d < 1.9) site.n_below_19++;
+            else if (d < 2.3) site.n_19_23++;
+            else if (d < 2.6) site.n_23_26++;
+            else if (d < 3.2) site.n_26_32++;
+            else if (d < 3.4) site.n_32_34++;
+            else if (d < 3.6) site.n_34_36++;
+            if (d < SHELL_RADIUS) {
+                if (chem.donor && chem.acceptor)
+                    site.donor_acceptors++;
+                else if (chem.donor)
+                    site.donors++;
+                else if (chem.acceptor)
+                    site.acceptors++;
+                site.neighbours += fmt::format("{}{} {:.2f}", site.neighbours.empty() ? "" : ", ",
+                                               AtomName(pr, pa, pch.name), d);
+            }
+            if (d < 2.6)
+                inner.push_back(nb.pos);
+
+            Ligand l{};
+            l.pos = nb.pos;
+            l.distance = d;
+            const double e = CoordinateError(pa.b_iso, d_min);
+            l.sigma = std::sqrt(site.sigma_xyz * site.sigma_xyz + e * e);
+            l.chem = chem;
+            l.nitrogen = pa.element == gemmi::El::N;
+            l.n7_off_plane = false;
+            if (pa.name == "N7" && d < 2.6) {
+                // the purine's ring, as the image the neighbour is, and its O6 or N6
+                std::vector<gemmi::Atom> ring;
+                const gemmi::Atom *o6n6 = nullptr;
+                for (const gemmi::Atom &ra : pr.atoms) {
+                    static const std::set<std::string> RING = {"N9", "C8", "N7", "C5", "C6", "N1", "C2", "N3", "C4"};
+                    if (RING.count(ra.name) && (!a.altloc || !ra.altloc || ra.altloc == a.altloc)) {
+                        ring.push_back(ra);
+                        ring.back().pos = gemmi::Position(nb.image.apply(ra.pos));
+                    }
+                    if (ra.name == "O6" || ra.name == "N6")
+                        o6n6 = &ra;
+                }
+                if (ring.size() == 9 && o6n6 != nullptr) {
+                    site.n7_ligands++;
+                    std::vector<gemmi::Atom *> ptrs;
+                    for (gemmi::Atom &ra : ring)
+                        ptrs.push_back(&ra);
+                    const double plane = std::fabs(gemmi::get_distance_from_plane(q, gemmi::find_best_plane(ptrs)));
+                    const double side = q.dist(gemmi::Position(nb.image.apply(o6n6->pos)));
+                    l.n7_off_plane = plane > N7_PLANE_MAX + l.sigma || side < N7_O6N6_MIN - l.sigma ||
+                                     side > N7_O6N6_MAX + l.sigma;
+                    site.n7_off_plane += l.n7_off_plane;
+                }
+            }
+            if (d < BVS_RADIUS)
+                for (int c = 0; c < BVS_CATIONS; c++) {
+                    // Following Brown & Altermatt (1985) Acta Cryst. B41, 244-247, occupancy weighted as
+                    // Muller et al. (2003) Acta Cryst. D59, 32-37
+                    l.valence[c] = pa.occ * std::exp(((l.nitrogen ? R0_N[c] : R0_O[c]) - d) / BV_B);
+                    site.bvs[c] += l.valence[c];
+                    vec[c] += (nb.pos - q) * (l.valence[c] / d);
+                }
+            ligands.push_back(l);
+        }
+        // A site bound to a cation shares that cation's other ligands as neighbours: the cis waters of a
+        // Mg(H2O)6 are 2.9 A from each other, and are no hydrogen bonds of theirs. The scores leave them out.
+        std::erase_if(ligands, [&](const Ligand &l) {
+            for (const gemmi::Position &c : bound_cations)
+                if (l.pos.dist(c) < 2.6)
+                    return true;
+            return false;
                 });
                 for (int c = 0; c < BVS_CATIONS; c++)
                     site.nvecsum[c] = site.bvs[c] > 0 ? vec[c].length() / site.bvs[c] : 0;
@@ -803,11 +854,9 @@ std::vector<SolventSite> AnalyseSolventSites(const gemmi::Structure &st, double
                 if (std::isfinite(anomalous_scale))
                     site.anomalous_sigma = map_sigma(maps.anomalous, maps.anomalous_rms, q);
 
-                site_ligands.push_back(std::move(ligands));
-                site_expected.push_back(expected_electrons);
-                out.push_back(site);
-            }
-        }
+                site_ligands[index] = std::move(ligands);
+                site_expected[index] = expected_electrons;
+    });
 
     WaterElectrons waters;
     std::vector<double> w;

image_analysis/structure_refinement/SolventSites.h

@@ -115,7 +115,9 @@ constexpr double SOLVENT_PEAK_REACH = 0.6;
 // Every water and single-atom ion of the first model of `st` (in its cell and space group; hydrogens are
 // ignored), in model order. `d_min`: the resolution the model was refined at, for the coordinate error
 // (0 for a model with no data).
-std::vector<SolventSite> AnalyseSolventSites(const gemmi::Structure &st, double d_min, const SolventMaps &maps);
+// The sites are measured on `nthreads` threads; the result does not depend on their number.
+std::vector<SolventSite> AnalyseSolventSites(const gemmi::Structure &st, double d_min, const SolventMaps &maps,
+                                             size_t nthreads = 1);
 
 // Whether the site is labelled a water; whether its label is one of the species decided here (water, Mg, Na, K,
 // Ca, Cl); whether the best-supported species is its label.

rugnux/ResultReport.cpp

@@ -515,7 +515,7 @@ namespace {
             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, and waters and ions relabelled where their surroundings and density\n"
-                    " support another species, each round a refinement kept only where it lowered R-free: {}\n"
+                    " support another species, each round kept only where refining its solvent lowered R-free: {}\n"
                     " added, {} of the model's taken out, {} relabelled, in {} round{} ({}). {} peaks are not solvent\n"
                     " (bonded to the model, beside a displaced atom, a truncated side chain, a molecule's or elongated\n"
                     " density, denser than S and P): reported, not filled. Each, with its reason, in\n"

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, 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 << "   --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 its solvent rebuilt - waters added at difference-map peaks with a hydrogen-bond partner, removed from negative difference density, and waters and ions relabelled where their shell and density support another species (to 2.8 A, up to two rounds within the weight search, each kept where it lowers R-free; --no-refine-solvent keeps the solvent as given) - 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 << "   --no-refine-solvent             Refine the placed model with its solvent as given: no waters built or removed, no ion relabelled (the solvent audit, REFINE_SOLVENT_ keys and <prefix>_solvent.json, is still made). Also with --mode refine" << std::endl;

tests/RefinementTest.cpp

@@ -378,16 +378,6 @@ TEST_CASE("RefinementRestraints_metal_links_outer_sphere", "[Refinement]") {
             found += (r.atom[0] == mg && r.atom[1] == outer) || (r.atom[0] == outer && r.atom[1] == mg);
         CHECK(found == 1);
     }
-    // decided once: the record left as a contact stays one where the refinement has moved the water into the inner
-    // shell, and the six taken stay bonds where one has moved out of it
-    const std::set<std::string> outer_records = BuildRefinementTopology(st).outer_records;
-    REQUIRE(outer_records.size() == 1);
-    res[6].atoms[0].pos = site + gemmi::Position(1, 1, 1) * (2.1 / std::sqrt(3.0));
-    res[0].atoms[0].pos = site + gemmi::Position(1, 0, 0) * 2.9;
-    const RefinementTopology moved = BuildRefinementTopology(st, false, &outer_records);
-    CHECK(moved.metal_links == 6);
-    CHECK(moved.outer_records == outer_records);
-    CHECK(BuildRefinementTopology(st).metal_links == 6);   // re-gated: the other record now refused instead
 }
 
 // Deposited models cut down to a ligand and the residues within 4.2 A of it (test_data/ligands, from the PDB entries
@@ -831,7 +821,9 @@ TEST_CASE("Refinement_restrained_recovers_model", "[Refinement]") {
         std::vector<gemmi::Atom *> refined = atoms_of(r.model->models[0]);
         std::erase_if(refined, [](const gemmi::Atom *a) { return a->is_hydrogen(); });   // the riding hydrogens
         const std::vector<gemmi::Atom *> ref = atoms_of(truth.models[0]);
-        REQUIRE(refined.size() == ref.size());
+        REQUIRE(refined.size() == ref.size());   // the solvent build's placeholders not among them
+        for (const gemmi::Chain &ch : r.model->models[0].chains)
+            CHECK(ch.name != "~");
         for (size_t i = 0; i < ref.size(); i++)
             rms += refined[i]->pos.dist_sq(ref[i]->pos);
         rms = std::sqrt(rms / ref.size());