From c4df332618baeedcf66f040759f3d07e223bf8f0 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 5 Nov 2026 15:51:05 +0100 Subject: [PATCH] Solvent build: write each new water beside the model, not a symmetry mate A difference peak is taken once per symmetry-equivalent set (the map's asymmetric unit) and was to be moved to its copy nearest the model's own atoms (NearestCopy) before a water is built there. That step took the nearest atom's NeighborSearch mark, whose position is already the symmetry image, and applied the image's operator to it a second time: outside the identity image the water was written beside a symmetry mate or a unit cell away (by a lattice translation for a two-fold, beside another mate for a three-, four- or six-fold). The copy is now taken against the atom as the model has it: the symmetry operator and lattice translation that put the peak nearest to that atom (UnitCell::find_nearest_image / find_nearest_pbc_position). The water's environment, and so the refinement, is the same; only the coordinates written move. New test: a peak beside a four-fold mate in P4 is written 2.8 A from the model's water (the old step wrote it a lattice vector away). 5reo (open arm, C2), --refine-solvent, share of the added waters whose nearest non-water, non-hydrogen model atom over all copies is in the written copy: from the model without waters 67/155 (43%) before, 154/155 after; from the deposited model 27/43 before, 42/43 after (the one left a near tie moved by refinement). R-work/R-free 0.2046/0.2335 -> 0.2046/0.2334 and 0.1979/0.2172 unchanged; 156 and 43 waters added either way. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UqXYdL4qBGYRiSXSGXfCLu --- docs/CPU_DATA_ANALYSIS_REFINEMENT.md | 4 +++- .../structure_refinement/SolventBuild.cpp | 13 +++++-------- tests/SolventBuildTest.cpp | 19 +++++++++++++++++++ 3 files changed, 27 insertions(+), 9 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index da42e52fc..ec43655ae 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -476,7 +476,9 @@ higher, one 0.0016 lower). `--refine-score-only` builds nothing, and the build a restrains the model's own (§15.5). - *Candidates* are the local maxima of mFo-DFc from 3 σ (phenix.refine's peak cutoff), each symmetry-equivalent - set once, moved to the copy nearest the model. A peak closer than 1.2 Å to a model atom is that atom's own + set once, moved to the copy nearest the model - the symmetry operator and lattice translation that put it + closest to an atom of the model as written, so a water built there is written beside the model, not beside a + symmetry mate or in another unit cell. A peak closer than 1.2 Å to a model atom is that atom's own residual and is passed over. - *Not solvent* - reported in `_solvent.json` with the reason, never filled: a peak 1.2-1.9 Å from a model atom (a missing or displaced atom); one beside a model atom (within 2.6 Å) that sits in mFo-DFc of -3 σ or diff --git a/image_analysis/structure_refinement/SolventBuild.cpp b/image_analysis/structure_refinement/SolventBuild.cpp index 8e4e17178..407c89ca9 100644 --- a/image_analysis/structure_refinement/SolventBuild.cpp +++ b/image_analysis/structure_refinement/SolventBuild.cpp @@ -132,18 +132,15 @@ std::vector DifferencePeaks(const gemmi::Grid &g, double rms, doubl return out; } -// The copy of `p` nearest to the model's own atoms (not their symmetry mates); false where no atom is within -// `reach`. +// The copy of `p` - the symmetry operator and lattice translation - nearest to the model's own atoms (not their +// symmetry mates); false where no atom is within `reach`. The mark's position is already the atom's image, so the +// copy is taken by the atom as the model has it: p's copy nearest to that atom. bool NearestCopy(gemmi::NeighborSearch &ns, const gemmi::UnitCell &cell, double reach, gemmi::Position &p) { gemmi::NeighborSearch::Mark *m = ns.find_nearest_atom(p, reach); if (m == nullptr) return false; - const gemmi::FTransform t = ns.get_image_transformation(m->image_idx); - const gemmi::Fractional fc = cell.fractionalize(p); - const gemmi::Fractional fi = t.apply(cell.fractionalize(m->pos)); - const gemmi::Fractional shifted(fc.x - std::round(fc.x - fi.x), fc.y - std::round(fc.y - fi.y), - fc.z - std::round(fc.z - fi.z)); - p = cell.orthogonalize(gemmi::FTransform(t.inverse()).apply(shifted)); + const gemmi::Position atom = m->to_cra(*ns.model).atom->pos; + p = cell.find_nearest_pbc_position(atom, p, cell.find_nearest_image(atom, p, gemmi::Asu::Any).sym_idx); return true; } diff --git a/tests/SolventBuildTest.cpp b/tests/SolventBuildTest.cpp index 4103145ea..ee6fae5ee 100644 --- a/tests/SolventBuildTest.cpp +++ b/tests/SolventBuildTest.cpp @@ -101,6 +101,25 @@ TEST_CASE("SolventBuild_peak_with_and_without_partner", "[Solvent]") { CHECK(r.edits.at(0).site == "HOH W1"); } +// A peak beside a symmetry mate of the model's water, in another unit cell, is written beside the water itself: the +// copy of the peak - symmetry operator and lattice translation - nearest to the model's own atoms. +TEST_CASE("SolventBuild_peak_written_beside_the_model", "[Solvent]") { + gemmi::Structure st = box(); + st.spacegroup_hm = "P 4"; + st.setup_cell_images(); + add_atom(st, "HOH", "O", gemmi::El::O, {6, 4, 6}); + // the water's mate by -y,x,z is at 20 6 6; the peak 2.8 A from it, only there (the map drawn without symmetry) + const gemmi::Position peak(20, 6, 8.8); + const gemmi::Grid diff = map_with(gemmi::UnitCell(EDGE, EDGE, EDGE, 90, 90, 90), {peak}, 6.0f); + const gemmi::Grid map = map_with(st.cell, {}, 0.0f, 1.5f); + const SolventMaps maps{&map, 1.0, &diff, 1.0, nullptr, 0}; + + const SolventRound r = PickSolvent(st, maps, {}, {}, {}, 300, 1); + REQUIRE(r.added == 1); + const SolventEdit &e = r.edits.at(0); + CHECK(e.pos.dist(gemmi::Position(6, 4, 6)) == Catch::Approx(2.8).margin(0.25)); +} + // A sulfate-like cluster - a peak with four more at 1.5 A on a tetrahedron - beside a water is no water: the centre // is reported as a tetrahedral group, its corners as atoms of one molecule; nothing is added. TEST_CASE("SolventBuild_tetrahedral_cluster_refused", "[Solvent]") {