diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md
index 781f80807..6b3a66af2 100644
--- a/THIRD_PARTY_NOTICES.md
+++ b/THIRD_PARTY_NOTICES.md
@@ -135,7 +135,8 @@ served frontend, so the shipped web UI carries its own attribution.
 * **CCP4 Monomer Library (LGPL-3.0)** is data, not code: the restraint dictionaries of the compounds
   rugnux's refinement restrains, copied unmodified per compound into
   `image_analysis/structure_refinement/monomers/` (`links_and_mods.cif` and `ener_lib.cif` cut down to
-  the entries those compounds use), with the library's own `COPYING`. The build embeds them in the binary
+  the entries those compounds use, and the library's metal-site distances of `metals.json` for its ions as
+  `metals.cif`), with the library's own `COPYING`. The build embeds them in the binary
   as text, and `monomers/update_monomers.py` regenerates the subset from the upstream library, so the
   copy can be replaced by a user with a version of their own.
 * **Qt (LGPL-3.0)** and **NVIDIA CUDA (EULA)** carry redistribution conditions beyond a copyright
diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md
index c58a481d4..82d2ac334 100644
--- a/docs/ACKNOWLEDGEMENT.md
+++ b/docs/ACKNOWLEDGEMENT.md
@@ -562,6 +562,15 @@ robust least squares - iteratively reweighted, with the pseudo-Huber weight of e
 so that a few reflections the model cannot explain do not set the scale of the rest: P. J. Huber, "Robust
 estimation of a location parameter" (1964), Ann. Math. Statist. 35, 73-101
 [doi:10.1214/aoms/1177703732](https://doi.org/10.1214/aoms/1177703732).
+The bonds of metal ions to their donor atoms that the model's LINK records name are restrained at the
+distances of the Monomer Library's own statistics of metal sites (its `metals.json`, embedded as
+`monomers/metals.cif`; per metal, donor element and coordination number), in line with the surveys of
+M. M. Harding, "Geometry of metal-ligand interactions in proteins" (2001), Acta Cryst. D57, 401-411
+[doi:10.1107/S0907444900019168](https://doi.org/10.1107/S0907444900019168) and "Small revisions to
+predicted distances around metal sites in proteins" (2006), Acta Cryst. D62, 678-682
+[doi:10.1107/S0907444906014594](https://doi.org/10.1107/S0907444906014594), and, for Mg²⁺, Na⁺ and K⁺ at
+nucleic acids, F. Leonarski, L. D'Ascenzo and P. Auffinger, "Mg2+ ions: do they bind to nucleobase
+nitrogens?" (2017), Nucleic Acids Res. 45, 987-1004 [doi:10.1093/nar/gkw1175](https://doi.org/10.1093/nar/gkw1175).
 
 **Unmodelled density** — the blobs `--refine` lists are the connected regions of the difference map above
 a cutoff, found by flood fill as Coot's "unmodelled blobs" finds them, in GEMMI's implementation
diff --git a/docs/RUGNUX_MODEL_AND_REFERENCE.md b/docs/RUGNUX_MODEL_AND_REFERENCE.md
index 9b6d58b11..adcb6e2be 100644
--- a/docs/RUGNUX_MODEL_AND_REFERENCE.md
+++ b/docs/RUGNUX_MODEL_AND_REFERENCE.md
@@ -242,7 +242,11 @@ nucleotides, ATP/ADP/AMP/GTP/GDP/GMP, DMSO, PEG fragments (PEG, PGE, PG4, 1PE),
 sulfate, phosphate, acetate, MPD, water and common ions (a subset of the CCP4 Monomer Library compiled into
 rugnux; nothing is read from disk). Any other residue is held fixed: it stays in the structure factors and
 the bulk-solvent mask, and the refined atoms are kept from clashing with it; `REFINE_RESIDUES_FIXED_NAMES`
-lists them. The dictionary's hydrogens are added and ride on their parents. The X-ray weight is chosen from a
+lists them. The dictionary's hydrogens are added and ride on their parents. Metal ions are restrained to the
+donor atoms (O, N, S) the model's LINK / `_struct_conn` records bind them to: by the dictionary's own Zn and Fe
+links to Cys and His where they apply, otherwise at the library's typical distance for the metal, the donor
+element and the number of such records of the ion (Mg-O 2.07 A, Na-O 2.42 A, K-O 2.82 A at their
+usual coordination). A contact without a record, or one to a symmetry mate, is not restrained as a bond. The X-ray weight is chosen from a
 ladder of complete refinements by R-free, among those whose bond and angle rmsd stay within 0.02 A and 2
 degrees; B-factor refinement is then tried at a few weights and kept only where it lowers R-free beyond
 0.001. The free set is the merge's own, so R-free after `--refine` judges the refinement against reflections it
diff --git a/image_analysis/structure_refinement/MonomerLibrary.cpp b/image_analysis/structure_refinement/MonomerLibrary.cpp
index df908fdd0..3c55b8fb3 100644
--- a/image_analysis/structure_refinement/MonomerLibrary.cpp
+++ b/image_analysis/structure_refinement/MonomerLibrary.cpp
@@ -3,8 +3,11 @@
 
 #include "MonomerLibrary.h"
 
+#include <algorithm>
 #include <cstring>
+#include <vector>
 
+#include "gemmi/numb.hpp"       // as_number
 #include "gemmi/read_cif.hpp"   // read_cif_from_memory
 
 // Generated by CMakeLists.txt from monomers/*.cif (MonomerLibraryData.cpp in the build tree).
@@ -13,27 +16,79 @@ extern const size_t embedded_monomer_file_count;
 
 namespace {
 
-gemmi::MonLib ReadEmbeddedMonomerLibrary() {
+struct MetalRow {
+    gemmi::El metal, ligand;
+    int coordination;
+    double distance, std;
+    int count;
+};
+
+struct EmbeddedLibrary {
     gemmi::MonLib monlib;
+    std::vector<MetalRow> metals;   // metals.cif
+};
+
+EmbeddedLibrary ReadEmbeddedMonomerLibrary() {
+    EmbeddedLibrary lib;
+    gemmi::MonLib &monlib = lib.monlib;
     for (size_t i = 0; i < embedded_monomer_file_count; i++) {
         const EmbeddedMonomerFile &f = embedded_monomer_files[i];
         std::string text;
         for (size_t k = 0; k < f.n_pieces; k++)
             text += f.pieces[k];
-        const gemmi::cif::Document doc = gemmi::read_cif_from_memory(text.data(), text.size(), f.name);
+        gemmi::cif::Document doc = gemmi::read_cif_from_memory(text.data(), text.size(), f.name);
         if (std::strcmp(f.name, "ener_lib.cif") == 0)
             monlib.ener_lib.read(doc);
-        else
+        else if (std::strcmp(f.name, "metals.cif") == 0) {
+            for (const auto &row : doc.blocks.at(0).find("_metal_coordination.",
+                                                         {"metal", "ligand", "coordination", "distance", "std", "count"}))
+                lib.metals.push_back({gemmi::find_element(row[0].c_str()), gemmi::find_element(row[1].c_str()),
+                                      gemmi::cif::as_int(row[2]), gemmi::cif::as_number(row[3]),
+                                      gemmi::cif::as_number(row[4]), gemmi::cif::as_int(row[5])});
+        } else
             monlib.read_monomer_doc(doc);
     }
-    return monlib;
+    return lib;
+}
+
+const EmbeddedLibrary &Library() {
+    static const EmbeddedLibrary lib = ReadEmbeddedMonomerLibrary();
+    return lib;
 }
 
 }  // namespace
 
 const gemmi::MonLib &EmbeddedMonomerLibrary() {
-    static const gemmi::MonLib monlib = ReadEmbeddedMonomerLibrary();
-    return monlib;
+    return Library().monlib;
+}
+
+bool EmbeddedMetalDistance(gemmi::El metal, gemmi::El ligand, int coordination, double &distance, double &sigma) {
+    // the rows of the pair resting on enough sites: the nearest at or below the coordination and above it
+    const MetalRow *lo = nullptr, *hi = nullptr;
+    for (const MetalRow &r : Library().metals) {
+        if (r.metal != metal || r.ligand != ligand || r.count < METAL_MIN_SITES)
+            continue;
+        if (r.coordination <= coordination && (lo == nullptr || r.coordination > lo->coordination))
+            lo = &r;
+        if (r.coordination > coordination && (hi == nullptr || r.coordination < hi->coordination))
+            hi = &r;
+    }
+    if (lo == nullptr && hi == nullptr)
+        return false;
+    if (lo != nullptr && (lo->coordination == coordination || hi == nullptr)) {
+        distance = lo->distance;
+        sigma = lo->std;
+    } else if (lo == nullptr) {
+        distance = hi->distance;
+        sigma = hi->std;
+    } else {
+        // between two coordination numbers: linear in the coordination
+        const double t = static_cast<double>(coordination - lo->coordination) / (hi->coordination - lo->coordination);
+        distance = lo->distance + t * (hi->distance - lo->distance);
+        sigma = lo->std + t * (hi->std - lo->std);
+    }
+    sigma = std::max(sigma, METAL_MIN_SIGMA);
+    return true;
 }
 
 bool InEmbeddedMonomerLibrary(const std::string &code) {
diff --git a/image_analysis/structure_refinement/MonomerLibrary.h b/image_analysis/structure_refinement/MonomerLibrary.h
index 1a21fe274..ec0e944d9 100644
--- a/image_analysis/structure_refinement/MonomerLibrary.h
+++ b/image_analysis/structure_refinement/MonomerLibrary.h
@@ -14,6 +14,15 @@
 // the list. Parsed from memory once per process; no file is read.
 const gemmi::MonLib &EmbeddedMonomerLibrary();
 
+// The ideal distance and its sigma of a bond from a metal ion to a donor atom (O, N, S) at a site of `coordination`
+// donors, from the library's statistics of metal sites (metals.json, embedded as monomers/metals.cif). Only rows
+// resting on at least METAL_MIN_SITES sites count: the row of that coordination number, else the two nearest
+// below and above it interpolated linearly, else the nearest one; the sigma is the row's standard deviation, at
+// least METAL_MIN_SIGMA. False where the pair has no such row.
+constexpr int METAL_MIN_SITES = 20;
+constexpr double METAL_MIN_SIGMA = 0.02;   // A
+bool EmbeddedMetalDistance(gemmi::El metal, gemmi::El ligand, int coordination, double &distance, double &sigma);
+
 // Whether `code` (a residue name) has a description in it.
 bool InEmbeddedMonomerLibrary(const std::string &code);
 
diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp
index 9772c78e0..b055ddc82 100644
--- a/image_analysis/structure_refinement/Refinement.cpp
+++ b/image_analysis/structure_refinement/Refinement.cpp
@@ -751,6 +751,7 @@ void PlaceRiders(Context &ctx, RefinementTopology &topo) {
 RestraintSet RestraintsNow(Context &ctx, const RefinementTopology &topo) {
     const RefinementStage stage("repulsions");
     RestraintSet set = topo.geometry;
+    set.Append(topo.metal_bonds);
     set.Append(BuildRepulsions(topo, ctx.Model(), ctx.cell, *ctx.sg, REPULSION_MARGIN, ctx.nthreads));
     set.Finalize(ctx.atoms.size());
     return set;
@@ -874,13 +875,16 @@ TrialOutcome RunTrial(TrialSlot slot, const std::vector<int> &refinable, const s
     {
         const RefinementStage s1("measure");
         t.geometry = MeasureGeometry(topo.geometry,
-                                     BuildRepulsions(topo, ctx.Model(), ctx.cell, *ctx.sg, 0.0, ctx.nthreads), ctx.XyzB());
+                                     BuildRepulsions(topo, ctx.Model(), ctx.cell, *ctx.sg, 0.0, ctx.nthreads), ctx.XyzB(),
+                                     &topo.metal_bonds);
     }
     t.within_bounds = t.geometry.rmsd_bonds <= MAX_RMSD_BONDS && t.geometry.rmsd_angles <= MAX_RMSD_ANGLES;
     out.line = fmt::format("Refinement, X-ray weight {:.3g}, B {}: R-work {:.4f}, R-free {:.4f}, rmsd bonds {:.4f} A, "
-                           "angles {:.2f} deg{}",
+                           "angles {:.2f} deg{}{}",
                            w, adp_w > 0 ? fmt::format("refined at weight {:.3g}", adp_w) : std::string("held"),
                            rf.r_work, rf.r_free, t.geometry.rmsd_bonds, t.geometry.rmsd_angles,
+                           t.geometry.metal_bonds > 0 ? fmt::format(", metal bonds {:.3f} A", t.geometry.rmsd_metal)
+                                                      : std::string(),
                            t.within_bounds ? "" : " (geometry out of bounds)");
     out.pos = ctx.Positions();
     out.b = ctx.B();
@@ -1020,11 +1024,13 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
         result.fixed_residues = topo.fixed_residues;
         result.fixed_names = topo.fixed_names;
         result.geometry_start = MeasureGeometry(topo.geometry, BuildRepulsions(topo, ctx.Model(), cell, sg, 0.0),
-                                                ctx.XyzB());
+                                                ctx.XyzB(), &topo.metal_bonds);
         finish();
-        logger.Info("Refinement: the model scored only; geometry rmsd bonds {:.4f} A (rmsZ {:.2f}), angles {:.2f} deg "
-                    "(rmsZ {:.2f})", result.geometry_start.rmsd_bonds, result.geometry_start.rmsz_bonds,
-                    result.geometry_start.rmsd_angles, result.geometry_start.rmsz_angles);
+        const GeometryStatistics &g = result.geometry_start;
+        logger.Info("Refinement: the model scored only; geometry rmsd bonds {:.4f} A (rmsZ {:.2f}, {} bonds), metal "
+                    "bonds {:.4f} A (rmsZ {:.2f}, {}; {} links made), both {:.4f} A; angles {:.2f} deg (rmsZ {:.2f})",
+                    g.rmsd_bonds, g.rmsz_bonds, g.bonds, g.rmsd_metal, g.rmsz_metal, g.metal_bonds, topo.metal_links,
+                    g.RmsdAllBonds(), g.rmsd_angles, g.rmsz_angles);
         return result;
     }
 
@@ -1047,18 +1053,18 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
             if (topo.refinable[i])
                 refinable.push_back(static_cast<int>(i));
         result.refined_atoms = static_cast<int>(refinable.size());
-        logger.Info("Refinement: {} residues restrained, {} held fixed{}; {} bonds, {} angles, {} torsions, {} chiral "
-                    "volumes, {} planes",
+        logger.Info("Refinement: {} residues restrained, {} held fixed{}; {} bonds, {} metal bonds ({} links made), {} "
+                    "angles, {} torsions, {} chiral volumes, {} planes",
                     topo.restrained_residues, topo.fixed_residues,
                     topo.fixed_names.empty() ? std::string() : " (" + Joined(topo.fixed_names) + ")",
-                    topo.bonds, topo.angles, topo.torsions, topo.chirals, topo.planes);
+                    topo.bonds, topo.metal_bonds.records.size(), topo.metal_links, topo.angles, topo.torsions, topo.chirals, topo.planes);
 
         ctx.engine->SetRestraints(1, topo.adp);
         {
             const RestraintSet set = RestraintsNow(ctx, topo);
             ctx.engine->SetRestraints(0, set);
             RestraintSet repulsions = BuildRepulsions(topo, ctx.Model(), cell, sg, 0.0, ctx.nthreads);
-            result.geometry_start = MeasureGeometry(topo.geometry, repulsions, ctx.XyzB());
+            result.geometry_start = MeasureGeometry(topo.geometry, repulsions, ctx.XyzB(), &topo.metal_bonds);
         }
         double w0 = 0;
         {
diff --git a/image_analysis/structure_refinement/RefinementRestraints.cpp b/image_analysis/structure_refinement/RefinementRestraints.cpp
index 7b3c0173c..ed5f22145 100644
--- a/image_analysis/structure_refinement/RefinementRestraints.cpp
+++ b/image_analysis/structure_refinement/RefinementRestraints.cpp
@@ -94,6 +94,89 @@ V3<double> FromImage(const RestraintRecord &r, V3<double> g) {
             r.rot[2] * g.x + r.rot[5] * g.y + r.rot[8] * g.z};
 }
 
+// A metal ion: a residue of one atom, of a metal element.
+bool IsMetalIon(const gemmi::CRA &cra) {
+    return cra.residue->atoms.size() == 1 && IsMetal(cra.atom->element);
+}
+
+bool IsMetalDonor(const gemmi::Element &el) {
+    return el == gemmi::El::O || el == gemmi::El::N || el == gemmi::El::S;
+}
+
+// The modification a donor's residue takes when it binds a metal: the proton the bond replaces, as the
+// library's own ZN-HISND, ZN-HISNE and ZN-CYS links have it.
+std::string MetalDonorModification(const std::string &resname, const std::string &atom) {
+    if (resname == "HIS" && atom == "ND1")
+        return "DEL-HD1";
+    if (resname == "HIS" && atom == "NE2")
+        return "DEL-HE2";
+    if (resname == "CYS" && atom == "SG")
+        return "CYSmod1";
+    return "";
+}
+
+// The model's LINK / struct_conn records between a metal ion and an O, N or S that no link of the dictionary
+// describes (its ZN-CYS, ZN-HISND, ZN-HISNE, FE-CYS, FE-HISND and FE-HISNE are taken as they are): a link of
+// one bond is made for each in `monlib`, at the distance of the library's statistics of metal sites
+// (EmbeddedMetalDistance) for the two elements and the site's coordination - the number of records between the
+// ion and an O, N or S - and the record is pointed at it, so that prepare_topology() restrains it like any other link.
+// Only the bond: no angles at the metal. Sites are taken from the records alone, as REFMAC takes links from
+// the coordinate file; a contact without a record stays a repulsion (BuildRepulsions, METAL_CONTACT).
+// Records to a symmetry mate are left out, as every bond to a symmetry mate is (BuildRefinementTopology).
+// Returns the number of links made.
+// 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) {
+    gemmi::Model &model = st.models.at(0);
+    struct Site {
+        gemmi::Connection *conn;
+        gemmi::CRA metal, donor;
+        bool metal_first;   // the metal is partner1
+    };
+    std::vector<Site> sites;
+    std::unordered_map<const gemmi::Atom *, int> coordination;
+    for (gemmi::Connection &conn : st.connections) {
+        gemmi::CRA a = model.find_cra(conn.partner1, true), b = model.find_cra(conn.partner2, true);
+        if (!a.atom || !b.atom)
+            continue;
+        const bool a_metal = IsMetalIon(a), b_metal = IsMetalIon(b);
+        if (a_metal == b_metal)
+            continue;
+        const gemmi::CRA &metal = a_metal ? a : b, &donor = a_metal ? b : a;
+        if (!IsMetalDonor(donor.atom->element))
+            continue;
+        // the site's coordination: every record of the ion, a symmetry mate's or one the dictionary describes too
+        coordination[metal.atom]++;
+        if (!conn.link_id.empty() || conn.asu == gemmi::Asu::Different ||
+            std::get<0>(monlib.match_link(*a.residue, a.atom->name, conn.partner1.altloc, *b.residue, b.atom->name,
+                                          conn.partner2.altloc)) != nullptr)
+            continue;
+        sites.push_back({&conn, metal, donor, a_metal});
+    }
+    int made = 0;
+    for (const Site &site : sites) {
+        double distance = 0, sigma = 0;
+        if (!EmbeddedMetalDistance(site.metal.atom->element.elem, site.donor.atom->element.elem,
+                                   coordination[site.metal.atom], distance, sigma))
+            continue;
+        const gemmi::Connection &conn = *site.conn;
+        gemmi::ChemLink link;
+        link.side1.comp = conn.partner1.res_id.name;
+        link.side2.comp = conn.partner2.res_id.name;
+        (site.metal_first ? link.side2 : link.side1).mod =
+            MetalDonorModification(site.donor.residue->name, site.donor.atom->name);
+        link.id = "metal-" + std::to_string(made);
+        link.name = link.id;
+        link.rt.bonds.push_back({gemmi::Restraints::AtomId{1, conn.partner1.atom_name},
+                                 gemmi::Restraints::AtomId{2, conn.partner2.atom_name},
+                                 gemmi::BondType::Metal, false, distance, sigma, distance, sigma});
+        monlib.links.emplace(link.id, link);
+        site.conn->link_id = link.id;
+        made++;
+    }
+    return made;
+}
+
 }  // namespace
 
 void RestraintSet::Add(RestraintRecord r, int n_atoms_of_r) {
@@ -247,6 +330,7 @@ RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens)
     out.monlib = std::make_shared<gemmi::MonLib>(EmbeddedMonomerLibrary());
     gemmi::MonLib &monlib = *out.monlib;
     gemmi::setup_entities(st);
+    out.metal_links = AddMetalLinks(st, monlib);
     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,
@@ -349,8 +433,12 @@ RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens)
         r.atom[1] = j;
         r.ideal = static_cast<float>(b.restr->value);
         r.sigma = static_cast<float>(b.restr->esd);
-        out.geometry.Add(r, 2);
-        out.bonds++;
+        if (out.metal[i] || out.metal[j])
+            out.metal_bonds.Add(r, 2);   // a metal's coordination, kept apart from the covalent bonds
+        else {
+            out.geometry.Add(r, 2);
+            out.bonds++;
+        }
         RestraintRecord rb;
         rb.kind = static_cast<int>(RestraintKind::BSimilarity);
         rb.atom[0] = i;
@@ -454,6 +542,7 @@ RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens)
     for (size_t i = 0; i < n; i++)
         out.excluded[i].assign(within3[i].begin(), within3[i].end());
     out.geometry.Finalize(n);
+    out.metal_bonds.Finalize(n);
     out.adp.Finalize(n);
     return out;
 }
@@ -549,8 +638,19 @@ void PlaceRidingHydrogens(RefinementTopology &topo) {
 }
 
 GeometryStatistics MeasureGeometry(const RestraintSet &geometry, const RestraintSet &repulsions,
-                                   const std::vector<std::array<double, 4>> &xyzb) {
+                                   const std::vector<std::array<double, 4>> &xyzb, const RestraintSet *metal_bonds) {
     GeometryStatistics s;
+    if (metal_bonds != nullptr) {
+        double sm = 0, zm = 0;
+        for (const RestraintRecord &r : metal_bonds->records) {
+            const double d = RestraintValue(*metal_bonds, r, xyzb) - r.ideal;
+            sm += d * d;
+            zm += d * d / (double(r.sigma) * r.sigma);
+            s.metal_bonds++;
+        }
+        s.rmsd_metal = s.metal_bonds > 0 ? std::sqrt(sm / s.metal_bonds) : 0;
+        s.rmsz_metal = s.metal_bonds > 0 ? std::sqrt(zm / s.metal_bonds) : 0;
+    }
     double sb = 0, sa = 0, sc = 0, zb = 0, za = 0;
     int nc = 0;
     for (const RestraintRecord &r : geometry.records) {
diff --git a/image_analysis/structure_refinement/RefinementRestraints.h b/image_analysis/structure_refinement/RefinementRestraints.h
index 969fc3a01..b4717c44b 100644
--- a/image_analysis/structure_refinement/RefinementRestraints.h
+++ b/image_analysis/structure_refinement/RefinementRestraints.h
@@ -4,6 +4,7 @@
 #pragma once
 
 #include <array>
+#include <cmath>
 #include <string>
 #include <vector>
 
@@ -48,6 +49,7 @@ double EvaluateRestraints(const RestraintSet &set, const std::vector<std::array<
 // What refinement knows of a model's chemistry.
 struct RefinementTopology {
     RestraintSet geometry;   // bonds, angles, torsions, chiral volumes, planes
+    RestraintSet metal_bonds;   // the bonds of metal ions to their donors (LINK records), apart from the rest
     RestraintSet adp;        // B of bonded and angle-related atoms alike
     std::vector<char> refinable;           // per atom: its residue is described by the dictionary
     // Per atom: the atoms the repulsion leaves out - itself and those up to three bonds away (sorted).
@@ -66,12 +68,14 @@ struct RefinementTopology {
     int fixed_residues = 0;
     std::vector<std::string> fixed_names;  // the residue names held fixed, each once
     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)
 };
 
 // The topology of the first model of `st` (hydrogens already removed) from the embedded monomer library
 // (MonomerLibrary.h), through gemmi's prepare_topology(): links between residues (peptide, cis-peptide,
 // nucleic-acid backbone, disulfides, the model's LINK records the library describes) and their
-// modifications included. A residue the library does not describe, or with an atom it does not know, is
+// modifications included; the LINK records of metal ions the library has no link for are restrained at the
+// library's metal-site distances (RefinementRestraints.cpp, AddMetalLinks). A residue the library does not describe, or with an atom it does not know, is
 // held fixed: it stays in the structure factors and the bulk-solvent mask, and it is restrained against.
 // Torsions: those with an esd, except the main-chain phi/psi and the sugar puckers given as alternatives.
 // With `hydrogens`, the dictionary's hydrogens are added to the residues it describes (waters excepted) and
@@ -95,6 +99,17 @@ struct GeometryStatistics {
     double rmsz_bonds = 0, rmsz_angles = 0;   // rms of (deviation / the dictionary's sigma)
     int bonds = 0, angles = 0;
     int clashes = 0;     // repulsions overlapping by more than 0.4 A
+    // The metal-donor bonds, apart: their sigma is the spread of real sites (0.04-0.1 A), so they sit further
+    // from ideal than covalent bonds do and are not what the bond rmsd's bound is about.
+    double rmsd_metal = 0, rmsz_metal = 0;
+    int metal_bonds = 0;
+    // Bonds and metal bonds together, as one rmsd (what REFMAC lists as its bond distances).
+    double RmsdAllBonds() const {
+        const int n = bonds + metal_bonds;
+        return n > 0 ? std::sqrt((rmsd_bonds * rmsd_bonds * bonds + rmsd_metal * rmsd_metal * metal_bonds) / n) : 0;
+    }
 };
+// `metal_bonds`: where given, its bonds are measured too (rmsd_metal), apart from `geometry`'s.
 GeometryStatistics MeasureGeometry(const RestraintSet &geometry, const RestraintSet &repulsions,
-                                   const std::vector<std::array<double, 4>> &xyzb);
+                                   const std::vector<std::array<double, 4>> &xyzb,
+                                   const RestraintSet *metal_bonds = nullptr);
diff --git a/image_analysis/structure_refinement/monomers/metals.cif b/image_analysis/structure_refinement/monomers/metals.cif
new file mode 100644
index 000000000..ef6358c3c
--- /dev/null
+++ b/image_analysis/structure_refinement/monomers/metals.cif
@@ -0,0 +1,197 @@
+data_metal_coordination
+loop_
+_metal_coordination.metal
+_metal_coordination.ligand
+_metal_coordination.coordination
+_metal_coordination.distance
+_metal_coordination.std
+_metal_coordination.count
+Ca O 4 2.3410 0.0508 11
+Ca O 5 2.3739 0.0757 42
+Ca O 6 2.3394 0.0530 1553
+Ca O 8 2.4677 0.0959 1383
+Ca O 10 2.6507 0.0849 20
+Ca O 11 2.2922 0.0000 1
+Ca O 12 2.3827 0.0231 12
+Ca O 13 2.4900 0.0362 9
+Ca N 4 2.3575 0.0404 16
+Ca N 5 2.3958 0.0865 27
+Ca N 6 2.4784 0.0664 88
+Ca N 8 2.5801 0.0620 140
+Ca N 9 2.2771 0.0055 2
+Ca N 12 2.5812 0.0464 4
+Ca S 8 2.9364 0.0301 64
+Cd O 4 2.2696 0.0541 1231
+Cd O 5 2.3075 0.0823 1290
+Cd O 6 2.3098 0.0733 9252
+Cd O 8 2.4498 0.1269 995
+Cd N 2 2.0928 0.0679 4
+Cd N 4 2.2704 0.0584 1034
+Cd N 5 2.3182 0.0729 978
+Cd N 6 2.3420 0.0567 7948
+Cd N 8 2.4168 0.0730 251
+Cd S 4 2.5344 0.0403 865
+Cd S 5 2.6051 0.0822 139
+Cd S 6 2.7112 0.0790 651
+Cd S 8 2.7723 0.0830 22
+Co O 2 1.8382 0.0023 4
+Co O 3 1.9027 0.0645 25
+Co O 4 1.9674 0.0748 2524
+Co O 5 2.0540 0.1288 2442
+Co O 6 2.0794 0.0873 31908
+Co O 8 2.0525 0.1404 23
+Co O 12 2.0981 0.0038 36
+Co N 2 1.8597 0.0228 19
+Co N 3 1.9920 0.0856 235
+Co N 4 1.9857 0.0829 3126
+Co N 5 2.0660 0.1097 2890
+Co N 6 2.0556 0.1079 23780
+Co N 7 1.8744 0.0681 30
+Co N 8 1.9732 0.0457 57
+Co N 9 1.9558 0.0000 1
+Co N 10 2.0748 0.0032 2
+Co N 12 1.9696 0.0309 12
+Co S 3 2.5588 0.1295 8
+Co S 4 2.2766 0.0681 327
+Co S 5 2.2863 0.1133 230
+Co S 6 2.3645 0.1482 797
+Co S 7 2.1350 0.0240 62
+Co S 8 2.2357 0.0388 110
+Co S 9 2.2874 0.0020 2
+Cu O 2 1.8209 0.0270 17
+Cu O 3 2.1634 0.2038 278
+Cu O 4 1.9484 0.0691 9401
+Cu O 5 2.0597 0.1858 18517
+Cu O 6 2.1733 0.2472 18687
+Cu N 2 1.8821 0.0563 144
+Cu N 3 1.9652 0.0713 2191
+Cu N 4 2.0018 0.0671 13836
+Cu N 5 2.0314 0.0969 20631
+Cu N 6 2.0861 0.1672 13385
+Cu N 7 2.0480 0.0464 10
+Cu S 2 2.1481 0.0143 21
+Cu S 3 2.2507 0.0481 651
+Cu S 4 2.3296 0.0931 2190
+Cu S 5 2.4463 0.2184 425
+Cu S 6 2.5901 0.2628 147
+Fe O 2 1.8469 0.0002 2
+Fe O 3 1.9645 0.1260 120
+Fe O 4 1.9111 0.0994 1017
+Fe O 5 1.9840 0.1333 1160
+Fe O 6 2.0322 0.0892 18340
+Fe O 7 2.0055 0.0302 14
+Fe O 8 2.2664 0.1337 37
+Fe N 2 1.8946 0.0223 22
+Fe N 3 1.9782 0.0914 396
+Fe N 4 1.9621 0.1558 1500
+Fe N 5 2.0733 0.1352 2715
+Fe N 6 2.0641 0.1154 22783
+Fe N 7 2.0293 0.1030 17
+Fe N 8 2.0743 0.1851 201
+Fe N 10 2.0233 0.0311 38
+Fe S 3 2.3024 0.0660 15
+Fe S 4 2.3269 0.0592 294
+Fe S 5 2.2888 0.0897 323
+Fe S 6 2.3543 0.1142 752
+Fe S 7 2.3237 0.0000 1
+Fe S 8 2.2729 0.0686 31
+Fe S 11 2.1977 0.0000 1
+K O 2 2.9059 0.5750 8
+K O 3 2.8089 0.1893 36
+K O 4 2.6188 0.0582 34
+K O 6 2.7546 0.1389 722
+K O 8 2.8184 0.0605 4680
+K O 9 2.8434 0.0817 18
+K O 10 2.8835 0.0759 87
+K O 11 2.6983 0.0465 4
+K O 12 2.8222 0.1108 16
+K O 16 2.7728 0.0428 16
+K N 3 2.8173 0.0237 3
+K N 4 2.7418 0.0241 3
+K N 6 2.8734 0.1215 408
+K N 8 2.9957 0.0534 1141
+K N 12 3.0862 0.0013 2
+K S 3 3.1901 0.0274 9
+K S 4 3.3220 0.1171 8
+K S 6 3.2685 0.0823 18
+K S 8 3.2406 0.1292 41
+Mg O 2 2.0554 0.0097 4
+Mg O 3 1.9937 0.0918 12
+Mg O 4 1.9685 0.0894 352
+Mg O 5 2.0575 0.0846 285
+Mg O 6 2.0741 0.0443 5665
+Mg O 8 2.0342 0.0034 2
+Mg O 12 2.0810 0.0003 12
+Mg N 2 2.0210 0.0000 2
+Mg N 3 2.0244 0.0434 50
+Mg N 4 2.0668 0.0526 431
+Mg N 5 2.1362 0.0752 213
+Mg N 6 2.1893 0.0689 532
+Mg N 7 2.0256 0.0253 12
+Mg S 4 2.4449 0.0066 4
+Mg S 6 2.6233 0.0847 18
+Mn O 3 2.1300 0.0000 1
+Mn O 4 1.9886 0.1364 712
+Mn O 5 2.0454 0.1325 2352
+Mn O 6 2.1235 0.1389 24192
+Mn O 8 2.3379 0.1153 187
+Mn N 2 1.9518 0.0006 2
+Mn N 3 2.1064 0.0836 30
+Mn N 4 2.1247 0.0998 452
+Mn N 5 2.1128 0.1560 1294
+Mn N 6 2.1915 0.1254 12485
+Mn N 7 2.0334 0.0117 8
+Mn N 8 2.1638 0.2290 63
+Mn S 4 2.4160 0.0516 82
+Mn S 5 2.3974 0.1473 65
+Mn S 6 2.5308 0.1543 759
+Mn S 8 2.7530 0.4496 4
+Na O 2 2.7451 0.2229 9
+Na O 3 2.5148 0.2189 356
+Na O 4 2.3601 0.1126 936
+Na O 5 2.3785 0.0889 1964
+Na O 6 2.4167 0.0890 10682
+Na O 8 2.5803 0.1583 1174
+Na O 11 2.4729 0.0682 11
+Na O 12 2.4485 0.0245 12
+Na O 16 2.4544 0.0868 64
+Na N 3 2.3817 0.0468 19
+Na N 4 2.4397 0.0763 110
+Na N 5 2.4662 0.0916 162
+Na N 6 2.5158 0.0858 535
+Na N 8 2.6107 0.1369 108
+Na S 3 3.0905 0.0704 9
+Na S 4 2.8085 0.0742 22
+Na S 5 2.9437 0.0994 24
+Na S 6 2.9463 0.0887 217
+Na S 8 3.0014 0.0301 8
+Ni O 2 1.7925 0.0333 3
+Ni O 3 2.0494 0.0771 147
+Ni O 4 1.8650 0.0474 1934
+Ni O 5 2.0181 0.1088 684
+Ni O 6 2.0680 0.0569 22549
+Ni N 2 1.8330 0.1294 18
+Ni N 3 1.9148 0.0996 263
+Ni N 4 1.9086 0.0578 7110
+Ni N 5 2.0509 0.0767 1323
+Ni N 6 2.0921 0.0531 18876
+Ni N 7 1.8852 0.1013 24
+Ni S 3 2.1662 0.0390 33
+Ni S 4 2.1727 0.0306 4582
+Ni S 5 2.2907 0.0973 198
+Ni S 6 2.4548 0.0691 1081
+Ni S 7 2.1761 0.0533 15
+Zn O 2 2.1082 0.0000 1
+Zn O 3 2.0096 0.0923 171
+Zn O 4 1.9578 0.0451 15632
+Zn O 5 2.0414 0.0799 10506
+Zn O 6 2.1095 0.0766 13008
+Zn N 2 1.8575 0.0462 21
+Zn N 3 2.0281 0.0767 383
+Zn N 4 2.0304 0.0529 10112
+Zn N 5 2.1003 0.0802 7734
+Zn N 6 2.1602 0.0689 6819
+Zn S 3 2.4083 0.0793 16
+Zn S 4 2.3391 0.0545 1189
+Zn S 5 2.4363 0.1220 308
+Zn S 6 2.5526 0.0669 103
diff --git a/image_analysis/structure_refinement/monomers/update_monomers.py b/image_analysis/structure_refinement/monomers/update_monomers.py
index 5b0b8d871..1e13afc30 100644
--- a/image_analysis/structure_refinement/monomers/update_monomers.py
+++ b/image_analysis/structure_refinement/monomers/update_monomers.py
@@ -13,8 +13,12 @@ is held fixed. To add one, add its code and run this again; CMake embeds every .
 The compound files are copied verbatim. links_and_mods.cif is cut down to the links and modifications
 that can apply to these compounds (the polymer links, the termini, the disulfide, ...), and ener_lib.cif
 to its atom-type table (_lib_atom: the van der Waals and ion radii and the hydrogen-bond types).
+metals.json, the library's statistics of metal-ligand distances, becomes metals.cif: for each metal ion of
+COMPOUNDS and each donor element of METAL_LIGANDS, the distance per coordination number (the restraints of
+metal sites the links of links_and_mods.cif do not describe).
 Requires the gemmi Python module."""
 
+import json
 import os
 import shutil
 import sys
@@ -29,6 +33,8 @@ ADDITIVES = ["DMS", "PEG", "PGE", "PG4", "1PE", "GOL", "EDO", "SO4", "PO4", "ACT
 WATER_AND_IONS = ["HOH", "NA", "K", "MG", "CA", "ZN", "MN", "FE", "FE2", "NI", "CU", "CO", "CD",
                   "CL", "BR", "IOD"]
 COMPOUNDS = AMINO_ACIDS + NUCLEOTIDES + NUCLEOTIDE_LIGANDS + ADDITIVES + WATER_AND_IONS
+# The donor elements of metal sites whose distances metals.cif carries.
+METAL_LIGANDS = ["O", "N", "S"]
 
 # Groups of the polymers above; links and modifications of other groups (sugars) are left out.
 GROUPS = {".", "peptide", "P-peptide", "M-peptide", "DNA/RNA", "NON-POLYMER"}
@@ -76,6 +82,30 @@ def trim_ener_lib(lib, out):
     new.write_file(out)
 
 
+def trim_metals(lib, out):
+    """metals.json -> metals.cif: one row per metal ion of COMPOUNDS, ligand element and coordination number,
+    with the mean distance, its standard deviation and the number of sites behind it."""
+    with open(os.path.join(lib, "metals.json")) as f:
+        table = json.load(f)["metal_coordination"]
+    metals = set()
+    for code in WATER_AND_IONS:
+        for atom in gemmi.make_chemcomp_from_block(gemmi.cif.read(compound_path(lib, code))[-1]).atoms:
+            if atom.el.is_metal:
+                metals.add(atom.el.name)
+    new = gemmi.cif.Document()
+    block = new.add_new_block("metal_coordination")
+    loop = block.init_loop("_metal_coordination.", ["metal", "ligand", "coordination", "distance", "std", "count"])
+    n = 0
+    for metal in sorted(metals):
+        for ligand in METAL_LIGANDS:
+            for row in table.get(metal, {}).get(ligand, []):
+                loop.add_row([metal, ligand, str(row["coord"]), "%.4f" % row["mean"], "%.4f" % row["std"],
+                              str(row["count"])])
+                n += 1
+    new.write_file(out)
+    return n
+
+
 def main():
     if len(sys.argv) != 2:
         sys.exit(__doc__)
@@ -88,8 +118,9 @@ def main():
         shutil.copyfile(compound_path(lib, code), os.path.join(here, code + ".cif"))
     n_links, n_mods = trim_links_and_mods(lib, os.path.join(here, "links_and_mods.cif"))
     trim_ener_lib(lib, os.path.join(here, "ener_lib.cif"))
+    n_metal_rows = trim_metals(lib, os.path.join(here, "metals.cif"))
     shutil.copyfile(os.path.join(lib, "COPYING"), os.path.join(here, "COPYING"))
-    print(f"{len(COMPOUNDS)} compounds, {n_links} links, {n_mods} modifications")
+    print(f"{len(COMPOUNDS)} compounds, {n_links} links, {n_mods} modifications, {n_metal_rows} metal distances")
 
 
 if __name__ == "__main__":
diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp
index c8bdd3636..b0ec872f1 100644
--- a/rugnux/ResultReport.cpp
+++ b/rugnux/ResultReport.cpp
@@ -476,6 +476,8 @@ namespace {
                 Add(s, KeyReal("REFINE_R_WORK", rf.r_work, "{:.4f}"));
                 Add(s, KeyReal("REFINE_R_FREE", rf.r_free, "{:.4f}"));
                 Add(s, KeyReal("REFINE_RMSD_BONDS", rf.geometry.rmsd_bonds, "{:.4f}"));
+                if (rf.geometry.metal_bonds > 0)
+                    Add(s, KeyReal("REFINE_RMSD_METAL_BONDS", rf.geometry.rmsd_metal, "{:.4f}"));
                 Add(s, KeyReal("REFINE_RMSD_ANGLES", rf.geometry.rmsd_angles, "{:.2f}"));
                 Add(s, KeyReal("REFINE_RMSD_BONDS_START", rf.geometry_start.rmsd_bonds, "{:.4f}", true));
                 Add(s, KeyReal("REFINE_RMSD_ANGLES_START", rf.geometry_start.rmsd_angles, "{:.2f}", true));
@@ -537,6 +539,8 @@ namespace {
             Add(s, KeyReal("REFINE_SCORE_FOM", sc.fom, "{:.4f}", true));
             Add(s, KeyReal("REFINE_SCORE_RMSD_BONDS", rf.geometry_start.rmsd_bonds, "{:.4f}", true));
             Add(s, KeyReal("REFINE_SCORE_RMSZ_BONDS", rf.geometry_start.rmsz_bonds, "{:.3f}", true));
+            Add(s, KeyReal("REFINE_SCORE_RMSD_METAL_BONDS", rf.geometry_start.rmsd_metal, "{:.4f}", true));
+            Add(s, KeyReal("REFINE_SCORE_RMSD_BONDS_ALL", rf.geometry_start.RmsdAllBonds(), "{:.4f}", true));
             Add(s, KeyReal("REFINE_SCORE_RMSD_ANGLES", rf.geometry_start.rmsd_angles, "{:.3f}", true));
             Add(s, KeyReal("REFINE_SCORE_RMSZ_ANGLES", rf.geometry_start.rmsz_angles, "{:.3f}", true));
             if (!sc.shells.empty()) {
diff --git a/tests/RefinementTest.cpp b/tests/RefinementTest.cpp
index 78c741632..7c6cae710 100644
--- a/tests/RefinementTest.cpp
+++ b/tests/RefinementTest.cpp
@@ -3,6 +3,8 @@
 
 #include <catch2/catch_all.hpp>
 
+#include <algorithm>
+#include <atomic>
 #include <cmath>
 #include <map>
 #include <random>
@@ -247,6 +249,88 @@ TEST_CASE("RefinementRestraints_ideal_residues", "[Refinement]") {
     CHECK(g.rmsd_chirals < 0.2);
 }
 
+// A zinc bound by four cysteines (the dictionary's own ZN-CYS link) and a magnesium by six waters (links made from
+// the library's metal-site statistics), both from LINK records: the bonds come out at the expected distances,
+// and a metal set off its site returns there under them alone.
+TEST_CASE("RefinementRestraints_metal_links", "[Refinement]") {
+    const gemmi::UnitCell cell(40, 40, 40, 90, 90, 90);
+    gemmi::Structure st = residue_crystal({"CYS", "CYS", "CYS", "CYS", "ZN", "HOH", "HOH", "HOH", "HOH", "HOH", "HOH",
+                                           "MG"}, cell, "P 1", 11);
+    double mg_o = 0, mg_sigma = 0;
+    REQUIRE(EmbeddedMetalDistance(gemmi::El::Mg, gemmi::El::O, 6, mg_o, mg_sigma));
+    CHECK(mg_o == Catch::Approx(2.07).margin(0.02));   // Mg-O of an octahedral site
+    constexpr double ZN_S = 2.34;                      // the ZN-CYS link
+    std::vector<gemmi::Residue> &res = st.models[0].chains[0].residues;
+    // the donors on a tetrahedron and an octahedron around their sites, each residue moved as a whole
+    const gemmi::Position zn_site(10, 10, 10), mg_site(25, 25, 25);
+    const double t = 1.0 / std::sqrt(3.0);
+    const gemmi::Position tetra[4] = {{t, t, t}, {t, -t, -t}, {-t, t, -t}, {-t, -t, t}};
+    const gemmi::Position octa[6] = {{1, 0, 0}, {-1, 0, 0}, {0, 1, 0}, {0, -1, 0}, {0, 0, 1}, {0, 0, -1}};
+    auto donor_to = [](gemmi::Residue &r, const std::string &atom, const gemmi::Position &target) {
+        const gemmi::Position shift = target - r.find_atom(atom, '*')->pos;
+        for (gemmi::Atom &a : r.atoms)
+            a.pos += shift;
+    };
+    for (int k = 0; k < 4; k++)
+        donor_to(res[k], "SG", zn_site + tetra[k] * ZN_S);
+    for (int k = 0; k < 6; k++)
+        donor_to(res[5 + k], "O", mg_site + octa[k] * mg_o);
+    // the metals off their sites by 0.3 A
+    res[4].atoms[0].pos = zn_site + gemmi::Position(0.2, -0.15, 0.17);
+    res[11].atoms[0].pos = mg_site + gemmi::Position(-0.1, 0.25, 0.12);
+    auto link = [&](int metal, int donor, const char *atom) {
+        gemmi::Connection c;
+        c.type = gemmi::Connection::MetalC;
+        c.asu = gemmi::Asu::Same;
+        c.partner1 = gemmi::AtomAddress("A", res[donor].seqid, res[donor].name, atom);
+        c.partner2 = gemmi::AtomAddress("A", res[metal].seqid, res[metal].name, res[metal].atoms[0].name);
+        st.connections.push_back(c);
+    };
+    for (int k = 0; k < 4; k++)
+        link(4, k, "SG");
+    for (int k = 0; k < 6; k++)
+        link(11, 5 + k, "O");
+
+    const RefinementTopology topo = BuildRefinementTopology(st);
+    CHECK(topo.metal_links == 6);   // the Mg-O links; Zn-S is the dictionary's
+    std::vector<gemmi::Atom *> atoms = atoms_of(st.models[0]);
+    const int zn = static_cast<int>(std::find(atoms.begin(), atoms.end(), &res[4].atoms[0]) - atoms.begin());
+    const int mg = static_cast<int>(std::find(atoms.begin(), atoms.end(), &res[11].atoms[0]) - atoms.begin());
+    // the metals' bonds kept apart from the covalent ones, none left among them
+    for (const RestraintRecord &r : topo.geometry.records)
+        CHECK_FALSE((r.kind == static_cast<int>(RestraintKind::Bond) &&
+                     (r.atom[0] == zn || r.atom[1] == zn || r.atom[0] == mg || r.atom[1] == mg)));
+    const RestraintSet &metal_bonds = topo.metal_bonds;
+    REQUIRE(metal_bonds.records.size() == 10);
+    for (const RestraintRecord &r : metal_bonds.records) {
+        const bool is_zn = r.atom[0] == zn || r.atom[1] == zn;
+        CHECK((is_zn || r.atom[0] == mg || r.atom[1] == mg));
+        CHECK(r.ideal == Catch::Approx(is_zn ? ZN_S : mg_o).margin(1e-4));
+    }
+
+    // the two metals moved down the gradient of their bonds, everything else held
+    std::vector<std::array<double, 4>> xyzb;
+    for (gemmi::Atom *a : atoms)
+        xyzb.push_back({a->pos.x, a->pos.y, a->pos.z, a->b_iso});
+    double e = EvaluateRestraints(metal_bonds, xyzb, nullptr), step = 1e-4;
+    for (int it = 0; it < 2000 && step > 1e-12; it++) {
+        std::vector<std::array<double, 4>> grad(xyzb.size(), {0, 0, 0, 0}), trial = xyzb;
+        EvaluateRestraints(metal_bonds, xyzb, &grad);
+        for (int a : {zn, mg})
+            for (int k = 0; k < 3; k++)
+                trial[a][k] -= step * grad[a][k];
+        const double e_trial = EvaluateRestraints(metal_bonds, trial, nullptr);
+        if (e_trial < e) {
+            xyzb = trial;
+            e = e_trial;
+            step *= 1.5;
+        } else
+            step *= 0.5;
+    }
+    for (const RestraintRecord &r : metal_bonds.records)
+        CHECK(RestraintValue(metal_bonds, r, xyzb) == Catch::Approx(r.ideal).margin(0.005));
+}
+
 TEST_CASE("RefinementRestraints_gradient", "[Refinement]") {
     const gemmi::UnitCell cell(14, 15, 16, 90, 90, 90);
     gemmi::Structure st = residue_crystal({"TYR", "DG"}, cell, "P 21 21 21", 5);
diff --git a/tools/battery/refine_parity.py b/tools/battery/refine_parity.py
index d7efed0fb..a00786b10 100644
--- a/tools/battery/refine_parity.py
+++ b/tools/battery/refine_parity.py
@@ -18,7 +18,7 @@ is the input F times REFMAC's overall scale, so FC_ALL_LS * F_in / F_out is the
 with. With --per-reflection the R of each data subset (resolution bins) is recomputed from both programs'
 model amplitudes on the same reflections, rugnux's being FC of its _refined_maps.mtz.
 
-  refine_parity.py --rugnux build/rugnux/rugnux --workdir DIR [--d-min D] [--refmac-extra 'make hydr no'] model data.mtz
+  refine_parity.py --rugnux build/rugnux/rugnux --workdir DIR [--d-min D] [--refmac-extra 'make hydr no'] [--gpu] model data.mtz
 """
 import argparse
 import json
@@ -76,10 +76,10 @@ def parse_refmac(text):
     return out
 
 
-def run_rugnux(rugnux, workdir, model, mtz, d_min):
+def run_rugnux(rugnux, workdir, model, mtz, d_min, gpu=False):
     prefix = os.path.join(workdir, "rugnux")
-    cmd = ["nice", rugnux, "--mode", "refine", "--refine-score-only", "--developer", "--no-gpu", "-N", "8", "-o", prefix,
-           "--model", model]
+    cmd = ["nice", rugnux, "--mode", "refine", "--refine-score-only", "--developer", "-N", "8", "-o", prefix,
+           "--model", model] + ([] if gpu else ["--no-gpu"])
     if d_min:
         cmd += ["--refine-d-min", f"{d_min:.3f}"]
     subprocess.run(cmd + [mtz], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
@@ -90,6 +90,8 @@ def run_rugnux(rugnux, workdir, model, mtz, d_min):
             "rmsd_bonds", "rmsz_bonds", "rmsd_angles", "rmsz_angles")}
     out["b_aniso"] = [float(x) for x in key("REFINE_SCORE_B_ANISO").split()]
     out["mask_spacing"] = float(key("REFINE_SCORE_MASK_SPACING"))
+    # REFMAC's bond rmsd counts the metal links among its bonds; rugnux reports them apart
+    out["rmsd_bonds"] = float(key("REFINE_SCORE_RMSD_BONDS_ALL"))
     return out
 
 
@@ -134,6 +136,7 @@ def main():
     ap.add_argument("--d-min", type=float, default=None)
     ap.add_argument("--refmac-extra", action="append", default=[])
     ap.add_argument("--no-refmac", action="store_true", help="reuse refmac.log in the workdir")
+    ap.add_argument("--gpu", action="store_true", help="rugnux on the GPU (the default is its CPU engine)")
     args = ap.parse_args()
     os.makedirs(args.workdir, exist_ok=True)
     model, data = os.path.abspath(args.model), os.path.abspath(args.data)
@@ -141,7 +144,7 @@ def main():
         ref = parse_refmac(open(os.path.join(args.workdir, "refmac.log"), errors="replace").read())
     else:
         ref = run_refmac(args.workdir, model, data, args.d_min, args.refmac_extra)
-    rug = run_rugnux(args.rugnux, args.workdir, model, data, args.d_min)
+    rug = run_rugnux(args.rugnux, args.workdir, model, data, args.d_min, args.gpu)
     print(f"{'term':<14} {'REFMAC':>10} {'rugnux':>10} {'diff':>10}")
     for k in ("r_work", "r_free", "k_overall", "b_iso", "k_sol", "b_sol", "fom", "nll_work", "nll_free",
               "rmsd_bonds", "rmsz_bonds", "rmsd_angles", "rmsz_angles"):