From 2a0e91bd6aac5b3a3b5bead3e8aadc07eca9f40a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 27 Oct 2026 18:39:47 +0100 Subject: [PATCH] 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 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/CPU_DATA_ANALYSIS_REFINEMENT.md | 22 +- docs/RUGNUX_REPORT.md | 4 +- .../structure_refinement/Refinement.cpp | 546 +++++++++++++----- .../structure_refinement/Refinement.h | 15 +- .../structure_refinement/RefinementGPU.cpp | 13 +- .../RefinementRestraints.cpp | 100 +++- .../RefinementRestraints.h | 15 +- .../structure_refinement/SolventBuild.cpp | 108 ++-- .../structure_refinement/SolventBuild.h | 24 +- .../structure_refinement/SolventSites.cpp | 437 +++++++------- .../structure_refinement/SolventSites.h | 4 +- rugnux/ResultReport.cpp | 2 +- rugnux/rugnux_cli.cpp | 2 +- tests/RefinementTest.cpp | 14 +- 14 files changed, 868 insertions(+), 438 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index a79a7d8d0..20ddb30cc 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/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 `_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 diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 719f7bf62..12fa07cf6 100644 --- a/docs/RUGNUX_REPORT.md +++ b/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 `_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 `_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 `_blobs.json`, strongest first, one object per blob: `kind`, diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index b23612f77..ce3d4e757 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/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 &merged, double fwt_data.v.push_back({ctx.refl.hkl[m], std::polar(static_cast(fwt), static_cast(phi))}); } result.maps = maps; - // The difference map beside the anomalous one below: neither reads the other. - std::future difference_map = std::async(std::launch::async, [&diff, &fwt_data, &result] { - const std::array 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 size = gemmi::get_size_for_hkl(diff, {{0, 0, 0}}, 3.0); + std::future difference_map = std::async(std::launch::async, [&diff, &size, &result] { result.difference_map = std::make_shared>(MapFromFPhi(gemmi::get_f_phi_on_grid(diff, size, true))); result.difference_map_rms = gemmi::calculate_data_statistics(result.difference_map->data).rms; + }); + std::future map = std::async(std::launch::async, [&fwt_data, &size, &result] { result.map = std::make_shared>(MapFromFPhi(gemmi::get_f_phi_on_grid(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 &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 MakeOwnSlot(const gemmi::Structure &before_topology, const Context &main, int device, - size_t nthreads, const RefinementSettings &settings, - const std::set *outer_records) { + size_t nthreads, const RefinementSettings &settings) { auto slot = std::make_unique(); Context &ctx = slot->ctx; ctx.st = before_topology; @@ -1066,7 +1070,7 @@ std::unique_ptr MakeOwnSlot(const gemmi::Structure &before_topology, co ctx.cancel = main.cancel; ctx.b_ceiling = main.b_ceiling; ctx.mask = std::make_unique(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 RunTrials(const std::vector &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 chain; // empty: not put in + std::vector seqid; + std::map 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 moved; + std::vector 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 &merged, const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double wavelength_A, - Logger &logger, const RefinementSettings &settings, - const std::set *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(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(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 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(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 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(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(std::move(model)); result.ok = true; result.seconds = std::chrono::duration(std::chrono::steady_clock::now() - start).count(); }; @@ -1410,7 +1552,7 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector refinable; @@ -1438,7 +1580,7 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector 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> 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 &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> 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 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 added; // the waters put in and still there, by name + std::set relabelled; // the sites relabelled, by their new names + std::vector not_solvent; + // per slot, the residue of each atom (the atoms stay where they are in memory) + std::vector> 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 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 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 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> undo(trial_slots.size()); + std::vector 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 &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 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(e.occupancy); + a.b_iso = static_cast(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(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 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(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 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 changed_pos = ctx.Positions(); + const std::vector 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 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(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(added.size()); + b.edits.insert(b.edits.end(), not_solvent.begin(), not_solvent.end()); + b.seconds = std::chrono::duration(std::chrono::steady_clock::now() - build_start).count(); + 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 &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 outer_records = result.outer_records; - std::set added; // the waters added and still there - std::set relabelled; // the sites relabelled, by their new names - std::vector not_solvent; // the not-solvent peaks of the model last picked on - for (int round = 1;; round++) { - const SolventMaps maps{result.map.get(), result.map_rms, result.difference_map.get(), - result.difference_map_rms, nullptr, 0}; - SolventRound picked = PickSolvent(*result.model, maps, result.solvent, added, 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(added.size()); - b.edits.insert(b.edits.end(), not_solvent.begin(), not_solvent.end()); - b.seconds = std::chrono::duration(std::chrono::steady_clock::now() - build_start).count(); - result.seconds = std::chrono::duration(std::chrono::steady_clock::now() - start).count(); - logger.Info("Solvent build: {} waters added, {} removed, {} 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; } diff --git a/image_analysis/structure_refinement/Refinement.h b/image_analysis/structure_refinement/Refinement.h index 77304fc91..d5c245401 100644 --- a/image_analysis/structure_refinement/Refinement.h +++ b/image_analysis/structure_refinement/Refinement.h @@ -60,9 +60,9 @@ struct RefinementSettings { std::vector final_ladder = {0.5, 1.0, 2.0}; double max_weight = 1.0 / 8; // the positions' X-ray weight at most, times the balanced one bool riding_hydrogens = true; // the dictionary's hydrogens added, riding on their parents - // Waters added at the difference map's peaks and removed where the maps no longer hold them, between refinements - // (SolventBuild.h), to SOLVENT_BUILD_D_MIN; 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 solvent; double b_ceiling = 0; // the highest B the refinement allowed - std::set 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 &merged, const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double wavelength_A, Logger &logger, const RefinementSettings &settings); diff --git a/image_analysis/structure_refinement/RefinementGPU.cpp b/image_analysis/structure_refinement/RefinementGPU.cpp index 753894ad3..1c1643c54 100644 --- a/image_analysis/structure_refinement/RefinementGPU.cpp +++ b/image_analysis/structure_refinement/RefinementGPU.cpp @@ -303,6 +303,7 @@ private: const int n = static_cast(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(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 engine_; std::vector atom_ptrs_; // the model's atoms at the last SetAtoms() std::vector 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 last_key_; PinnedVector atoms_; // pinned: a B rung sends them every evaluation (37 MB on a large model) PinnedVector> pos_; std::vector changed_; diff --git a/image_analysis/structure_refinement/RefinementRestraints.cpp b/image_analysis/structure_refinement/RefinementRestraints.cpp index cad4e2160..5397c4ca4 100644 --- a/image_analysis/structure_refinement/RefinementRestraints.cpp +++ b/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 *outer_given, - std::set &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 sites; std::unordered_map 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::setname, 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 &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(donors.size()), distance, sigma)) + continue; + RestraintRecord r; + r.kind = static_cast(RestraintKind::Bond); + r.atom[0] = j; + r.atom[1] = i; + r.ideal = static_cast(distance); + r.sigma = static_cast(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 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 *outer_records) { +RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens) { RefinementTopology out; gemmi::Model &model = st.models.at(0); out.monlib = std::make_shared(EmbeddedMonomerLibrary()); @@ -453,8 +503,7 @@ RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens, std::vector 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(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> first(model.chains.size()); std::vector pos; + std::vector 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(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; diff --git a/image_analysis/structure_refinement/RefinementRestraints.h b/image_analysis/structure_refinement/RefinementRestraints.h index 0a2613342..d51d4dcb7 100644 --- a/image_analysis/structure_refinement/RefinementRestraints.h +++ b/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 outer_records; // and which (MetalRecordKey) // The atoms on special positions (FindSpecialPositions), placed exactly on them; per atom its entry, -1 none. std::vector special; std::vector 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 *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 &donors); + // The riding hydrogens placed again on their parents as they are now (gemmi's place_hydrogens_on_all_atoms). void PlaceRidingHydrogens(RefinementTopology &topo); diff --git a/image_analysis/structure_refinement/SolventBuild.cpp b/image_analysis/structure_refinement/SolventBuild.cpp index e8d02e8f6..93970c93d 100644 --- a/image_analysis/structure_refinement/SolventBuild.cpp +++ b/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 InnerSphereDonors(const gemmi::Model &model, const gemmi::Residue &own, const gemmi::Atom &metal) { + std::vector 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(c), static_cast(r), static_cast(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 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 diff --git a/image_analysis/structure_refinement/SolventBuild.h b/image_analysis/structure_refinement/SolventBuild.h index 818ab642b..36fe3c4d4 100644 --- a/image_analysis/structure_refinement/SolventBuild.h +++ b/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 &built, const std::set &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 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 diff --git a/image_analysis/structure_refinement/SolventSites.cpp b/image_analysis/structure_refinement/SolventSites.cpp index 7cb279b0d..1d1b19c12 100644 --- a/image_analysis/structure_refinement/SolventSites.cpp +++ b/image_analysis/structure_refinement/SolventSites.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -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> &wanted, size_t nthreads) { + std::map> made; + for (const auto &[sc, b] : wanted) + made[KeyOf(sc, b)]; + std::vector> *> todo; + for (auto &entry : made) + if (!tables_.count(entry.first)) + todo.push_back(&entry); + ParallelFor(static_cast(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 &Table(const Scatterer &sc, double b) { - const int bi = static_cast(std::lround(std::clamp(b, 0.0, 500.0))); - const auto key = std::make_tuple(static_cast(sc.el), static_cast(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 t(static_cast(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 &t, double r) { @@ -166,9 +177,23 @@ public: } private: + using Key = std::tuple; // element, charge, hydrogens, B to the nearest A^2 + static Key KeyOf(const Scatterer &sc, double b) { + return {static_cast(sc.el), static_cast(sc.charge), sc.hydrogens, + static_cast(std::lround(std::clamp(b, 0.0, 500.0)))}; + } + std::vector Make(const Key &key) const { + const Scatterer sc{static_cast(std::get<0>(key)), static_cast(std::get<1>(key)), + std::get<2>(key)}; + std::vector t(static_cast(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::vector> tables_; + std::mutex m_; + std::map> 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(s.best); } -std::vector AnalyseSolventSites(const gemmi::Structure &st, double d_min, const SolventMaps &maps) { +std::vector AnalyseSolventSites(const gemmi::Structure &st, double d_min, const SolventMaps &maps, + size_t nthreads) { std::vector out; if (st.models.empty()) return out; @@ -523,6 +549,14 @@ std::vector AnalyseSolventSites(const gemmi::Structure &st, double }; const bool density = maps.map != nullptr && d_min > 0; RadialDensity radial(d_min); + if (density) { + std::vector> 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 AnalyseSolventSites(const gemmi::Structure &st, double std::map atom_calibration; std::pair global_calibration{NAN, NAN}; if (density) { - std::vector c; + std::vector> 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 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 v(atoms.size()); + ParallelFor(static_cast(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 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 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> 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(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 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 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(ci) && nb.residue == static_cast(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 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 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(ci) && nb.residue == static_cast(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(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(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 ligands; - std::vector bound_cations; // other cations closer than 2.6 A - std::vector 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(ci) && nb.residue == static_cast(ri) && - (nb.atom != static_cast(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 ligands; + std::vector bound_cations; // other cations closer than 2.6 A + std::vector 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(ci) && nb.residue == static_cast(ri) && + (nb.atom != static_cast(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 ring; - const gemmi::Atom *o6n6 = nullptr; - for (const gemmi::Atom &ra : pr.atoms) { - static const std::set 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 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 ring; + const gemmi::Atom *o6n6 = nullptr; + for (const gemmi::Atom &ra : pr.atoms) { + static const std::set 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 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 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 w; diff --git a/image_analysis/structure_refinement/SolventSites.h b/image_analysis/structure_refinement/SolventSites.h index 2306670c6..b2cd024cb 100644 --- a/image_analysis/structure_refinement/SolventSites.h +++ b/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 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 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. diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 23d976f89..6719a6d32 100644 --- a/rugnux/ResultReport.cpp +++ b/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" diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index b9679f5c9..7ec595a8e 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -256,7 +256,7 @@ void print_usage() { std::cout << " --scaling-iterations Cap on the per-frame scaling iterations; the loop stops earlier once the scales settle (default: 100)" << std::endl; std::cout << " -z, --reference Reference data of the same crystal form: an MTZ, or a PDB structure-factor mmCIF (e.g. 1abc-sf.cif, first merged block with a usable column; R-free set from _refln.status f/o, else pdbx_r_free_flag). Gzipped or not, recognised by content. --reference-mtz is the same option" << std::endl; std::cout << " --reference-column