Refinement restraints: contacts with lattice translates are kept, and an atom's own symmetry contact once

Contacts() dropped every pair found in image 0 with a nonzero lattice shift ("the identity term
covers the contact"), but that pair is a contact with a translated copy of j, not with j itself:
in P1 no intermolecular repulsion or B-contact restraint was made at all, and in every group the
contacts along a, b and c were missing (MeasureGeometry's clash count with them). Such a pair is
now kept like a symmetry mate's, from the lower atom of the two.

An atom against its own image x -> T x is also met through T^-1 at the same distance (another
operator, or the same one at another lattice shift for a screw), and both were kept: atoms near
a 3-, 4- or 6-fold or a screw axis were repelled twice. For i == j the image first in fractional
coordinates of T x and T^-1 x is kept; an involution is kept once as before.

Test: in P1 two waters 2.4 A apart across a cell face get one repulsion and one B contact (none
before); a water 1.4 A from a 3-fold axis gets one repulsion with its images (two before).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit 0d6e4d0491e655a49d59b709e1398843506d162a)

image_analysis/structure_refinement/RefinementRestraints.cpp

@@ -576,7 +576,7 @@ RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens)
 
 namespace {
 
-// The pairs of atoms of `model`, symmetry mates included, closer than `reach` and not excluded, where at least one of
+// The pairs of atoms of `model`, symmetry mates and lattice translates included, closer than `reach` and not excluded, where at least one of
 // the two is refinable, neither is a hydrogen and the two are not alternative conformations of each other: `make`
 // turns each into a record (the image of j as x -> rot x + tran, Cartesian, and the distance), or says to leave it out.
 // The search runs on `nthreads` threads; the list is in atom order whatever their number.
@@ -613,19 +613,34 @@ RestraintSet Contacts(const RefinementTopology &topo, gemmi::Model &model, const
                     return;
                 if (!topo.refinable[i] && !topo.refinable[j])
                     return;
-                if (m.image_idx == 0 && (static_cast<int>(i) >= j ||
-                                         std::binary_search(topo.excluded[i].begin(), topo.excluded[i].end(), j)))
-                    return;
-                if (m.image_idx != 0 && static_cast<int>(i) > j)
-                    return;
                 if (topo.altloc[i] && topo.altloc[j] && topo.altloc[i] != topo.altloc[j])
                     return;
                 // The image of j nearest to i, as a Cartesian x -> R x + t.
                 const gemmi::FTransform ft = ns.get_image_transformation(m.image_idx);
                 const gemmi::Fractional fj = ft.apply(cell.fractionalize(pos[j]));
                 const gemmi::Vec3 shift(std::round(fi.x - fj.x), std::round(fi.y - fj.y), std::round(fi.z - fj.z));
-                if (m.image_idx == 0 && (shift.x != 0 || shift.y != 0 || shift.z != 0))
-                    return;   // the same copy's own lattice translate: the identity term covers the contact
+                const bool same_copy = m.image_idx == 0 && shift.x == 0 && shift.y == 0 && shift.z == 0;
+                // Within the copy itself: each pair once, the excluded ones left out. With any other copy
+                // (a symmetry mate or a lattice translate): the pair is met again from j, so it is kept from
+                // the lower atom only.
+                if (same_copy && (static_cast<int>(i) >= j ||
+                                  std::binary_search(topo.excluded[i].begin(), topo.excluded[i].end(), j)))
+                    return;
+                if (!same_copy && static_cast<int>(i) > j)
+                    return;
+                // An atom against its own image x -> T x is met again through T^-1, at the same distance:
+                // the one of the two images that comes first in fractional coordinates is kept.
+                if (static_cast<int>(i) == j) {
+                    gemmi::Transform t = ft;
+                    t.vec += shift;
+                    const gemmi::Vec3 there = t.apply(fi), back = t.inverse().apply(fi);
+                    for (int a = 0; a < 3; a++) {
+                        if (back.at(a) < there.at(a) - 1e-6)
+                            return;
+                        if (back.at(a) > there.at(a) + 1e-6)
+                            break;
+                    }
+                }
                 const gemmi::Mat33 rot = cell.orth.mat.multiply(ft.mat).multiply(cell.frac.mat);
                 const gemmi::Vec3 tran = cell.orth.mat.multiply(ft.vec + shift);
                 const gemmi::Position image(rot.multiply(pos[j]) + tran);

tests/RefinementTest.cpp

@@ -507,6 +507,36 @@ TEST_CASE("RefinementRestraints_gradient", "[Refinement]") {
 #endif
 }
 
+// Contacts with other copies of the molecule: in P1 two waters 2.4 A apart across the a face (17.6 A apart within
+// the cell) repel each other once, and a water 1.4 A from a 3-fold axis, 2.4 A from its images x -> T x and
+// x -> T^-1 x (one contact per asymmetric unit), repels its image once.
+TEST_CASE("RefinementRestraints_lattice_and_symmetry_contacts", "[Refinement]") {
+    auto oxygens = [](gemmi::Structure &st, const std::vector<gemmi::Fractional> &where) {
+        std::vector<gemmi::Atom *> atoms = atoms_of(st.models[0]);
+        REQUIRE(atoms.size() == where.size());
+        for (size_t i = 0; i < where.size(); i++)
+            atoms[i]->pos = st.cell.orthogonalize(where[i]);
+    };
+    {
+        const gemmi::UnitCell cell(20, 21, 22, 90, 90, 90);
+        gemmi::Structure st = residue_crystal({"HOH", "HOH"}, cell, "P 1", 1);
+        oxygens(st, {gemmi::Fractional(0.01, 0.5, 0.5), gemmi::Fractional(0.89, 0.5, 0.5)});
+        const RefinementTopology topo = BuildRefinementTopology(st);
+        const RestraintSet rep = BuildRepulsions(topo, st.models[0], cell, *st.find_spacegroup(), 0.0);
+        REQUIRE(rep.records.size() == 1);
+        CHECK(rep.records[0].tran[0] == Catch::Approx(-20.0f));
+        CHECK(BuildBContacts(topo, st.models[0], cell, *st.find_spacegroup()).records.size() == 1);
+    }
+    {
+        const gemmi::UnitCell cell(20, 20, 22, 90, 90, 120);
+        gemmi::Structure st = residue_crystal({"HOH"}, cell, "P 3", 1);
+        oxygens(st, {cell.fractionalize(gemmi::Position(1.4, 0.0, 11.0))});
+        const RefinementTopology topo = BuildRefinementTopology(st);
+        const RestraintSet rep = BuildRepulsions(topo, st.models[0], cell, *st.find_spacegroup(), 0.0);
+        CHECK(rep.records.size() == 1);
+    }
+}
+
 // A complete model whose every atom is off by a Gaussian error of sigma_x per axis: D(s) = <cos(2 pi s.dx)> =
 // exp(-2 pi^2 sigma_x^2 s^2) (Luzzati). With a negligible measurement error D comes out at that; with an error of
 // half the rms amplitude of the outermost shell it falls below it there, where the error is a large part of the