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]") {