Refinement: common cofactors and modified RNA nucleotides in the embedded dictionary; heme iron restrained from LINK records

The embedded CCP4 Monomer Library subset grows from 62 to 84 compounds:
- the modified tRNA/rRNA nucleotides PSU, 5MC, 1MA, OMG, OMC, OMU, H2U, M2G, 7MG, 2MG, 5MU (RNA group,
  joined by the library's backbone link "p");
- the cofactors NAD, NAI, NAP, NDP, FAD, FMN, HEM, HEC, SAM, SAH, and dioxygen (OXY). CMO is left out:
  its dictionary has C-O at 1.25 A, where bound CO is near 1.15 A.
links_and_mods.cif gains the links these compounds use: HEC-CYS1/HEC-CYS2 (heme C thioethers),
CYS-HEMB/CYS-HEMC, CYS-FAD and HIS-FAD1/2, with their modifications.

AddMetalLinks now takes a metal atom of any residue, not only a one-atom ion, so a heme iron's LINK records
to His NE2 / Met SD / water become bonds at the library's metal-site distances; the site's coordination counts
the metal's bonds in its own dictionary (a heme's four pyrrole N) as well as the records. A metal within a
ligand is also a metal for the repulsions (1.9 A to N, O, S) and its dictionary bonds are kept with the metal
bonds. Note HEC.cif's own Fe-N(pyrrole) ideal is 1.90 A (HEM: 2.04 A), so a cytochrome c reports a metal-bond
rmsd near 0.1 A.

Test RefinementRestraints_deposited_ligands builds the topology of eleven deposited models cut down to the
ligand and its surroundings (tests/test_data/ligands, 66 KB, from PDB 1EHZ, 8BKH, 6K9I, 2W3V, 5VML, 6F3P,
8CON, 3TK0, 7TCD, 6W4H, 5EPE): every ligand restrained, bond rmsd < 0.03 A, the RNA backbone bonded to each
modified nucleotide, the heme Fe-His bonds, and the HEC thioethers. rugnux --mode refine on the full entries
against their deposited structure factors: R-free within +-0.003 of the build where these were held fixed,
ligand bond rmsd vs the dictionary from 0.004-0.058 A down to 0.002-0.042 A. rugnux grows by 425 KiB.

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

tests/RefinementTest.cpp

@@ -10,6 +10,7 @@
 #include <random>
 #include <vector>
 
+#include <gemmi/mmread_gz.hpp>
 #include <gemmi/symmetry.hpp>
 
 #include "../common/CUDAWrapper.h"
@@ -332,6 +333,103 @@ TEST_CASE("RefinementRestraints_metal_links", "[Refinement]") {
         CHECK(RestraintValue(metal_bonds, r, xyzb) == Catch::Approx(r.ideal).margin(0.005));
 }
 
+// Deposited models cut down to a ligand and the residues within 4.2 A of it (test_data/ligands, from the PDB entries
+// their names start with; the tRNA as its modified nucleotides and their neighbours in the chain): every ligand is
+// restrained, its own restraints included, the deposited geometry agrees with the dictionary, and the links the
+// ligands need are made - the RNA backbone to and from each modified nucleotide, the iron of a heme to its histidine
+// (and a heme C's to its methionine) from the LINK records, and a heme C's two thioethers to its cysteines.
+TEST_CASE("RefinementRestraints_deposited_ligands", "[Refinement]") {
+    const std::vector<std::string> modified_nucleotides = {"PSU", "5MC", "1MA", "OMG", "OMC", "H2U", "M2G", "7MG",
+                                                           "2MG", "5MU"};
+    const std::vector<std::pair<std::string, std::vector<std::string>>> cases = {
+        {"1ehz_psu", modified_nucleotides}, {"8bkh_hem", {"HEM"}}, {"6k9i_hec", {"HEC"}}, {"2w3v_ndp", {"NDP"}},
+        {"5vml_nap", {"NAP"}}, {"6f3p_nad", {"NAD"}}, {"8con_nai", {"NAI"}}, {"3tk0_fad", {"FAD"}},
+        {"7tcd_fmn", {"FMN"}}, {"6w4h_sam", {"SAM"}}, {"5epe_sah", {"SAH"}}};
+    for (const auto &[name, ligands] : cases) {
+        DYNAMIC_SECTION(name) {
+            gemmi::Structure st = gemmi::read_structure_gz("../../tests/test_data/ligands/" + name + ".pdb.gz");
+            const RefinementTopology topo = BuildRefinementTopology(st);
+            std::vector<gemmi::Atom *> atoms = atoms_of(st.models[0]);
+            std::vector<const gemmi::Residue *> residue_of;
+            for (const gemmi::Chain &ch : st.models[0].chains)
+                for (const gemmi::Residue &r : ch.residues)
+                    for (size_t k = 0; k < r.atoms.size(); k++)
+                        residue_of.push_back(&r);
+            auto index_of = [&](const gemmi::Residue &r, const char *atom) {
+                const gemmi::Atom *a = r.find_atom(atom, '*');
+                return a == nullptr ? -1 : static_cast<int>(std::find(atoms.begin(), atoms.end(), a) - atoms.begin());
+            };
+            auto bonded = [&](const RestraintSet &set, int i, int j) {
+                for (const RestraintRecord &r : set.records)
+                    if (r.kind == static_cast<int>(RestraintKind::Bond) &&
+                        ((r.atom[0] == i && r.atom[1] == j) || (r.atom[0] == j && r.atom[1] == i)))
+                        return true;
+                return false;
+            };
+
+            std::map<std::string, int> restraints, found;
+            for (const RestraintRecord &r : topo.geometry.records)
+                if (r.kind != static_cast<int>(RestraintKind::Plane))
+                    restraints[residue_of[r.atom[0]]->name]++;
+            for (size_t i = 0; i < atoms.size(); i++)
+                if (std::find(ligands.begin(), ligands.end(), residue_of[i]->name) != ligands.end()) {
+                    found[residue_of[i]->name]++;
+                    CHECK(topo.refinable[i]);
+                }
+            for (const std::string &code : ligands) {
+                INFO(code);
+                CHECK(found[code] > 0);
+                CHECK(restraints[code] > 0);
+                CHECK(std::find(topo.fixed_names.begin(), topo.fixed_names.end(), code) == topo.fixed_names.end());
+            }
+
+            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});
+            const GeometryStatistics g = MeasureGeometry(topo.geometry, {}, xyzb, &topo.metal_bonds);
+            CHECK(g.rmsd_bonds < 0.03);
+            CHECK(g.rmsd_angles < 3.0);
+            CHECK(g.rmsd_metal < 0.1);
+
+            if (name == "1ehz_psu") {
+                // each modified nucleotide bonded to the residue before it in the chain, where that one is there
+                int links = 0;
+                for (const gemmi::Residue &r : st.models[0].chains[0].residues) {
+                    const gemmi::Residue *prev = &r == &st.models[0].chains[0].residues.front() ? nullptr : &r - 1;
+                    if (prev == nullptr || prev->seqid.num != r.seqid.num - 1 ||
+                        std::find(ligands.begin(), ligands.end(), r.name) == ligands.end())
+                        continue;
+                    INFO(r.name << " " << *r.seqid.num);
+                    CHECK(bonded(topo.geometry, index_of(*prev, "O3'"), index_of(r, "P")));
+                    links++;
+                }
+                CHECK(links >= 10);
+            }
+            if (name == "8bkh_hem" || name == "6k9i_hec") {
+                const gemmi::Residue *heme = nullptr, *his = nullptr;
+                for (const gemmi::Residue &r : st.models[0].chains[0].residues) {
+                    if (r.name == ligands[0])
+                        heme = &r;
+                    if (r.name == "HIS" && r.seqid.num == (name == "8bkh_hem" ? 93 : 18))
+                        his = &r;
+                }
+                REQUIRE(heme != nullptr);
+                REQUIRE(his != nullptr);
+                const int fe = index_of(*heme, "FE");
+                CHECK(topo.metal[fe]);
+                CHECK(topo.metal_links == (name == "8bkh_hem" ? 1 : 2));   // the histidine; and the methionine
+                CHECK(bonded(topo.metal_bonds, fe, index_of(*his, "NE2")));
+                CHECK(bonded(topo.metal_bonds, fe, index_of(*heme, "NA")));   // the porphyrin's own, among the metal's
+                if (name == "6k9i_hec")
+                    for (auto [cys, carbon] : {std::pair{14, "CAB"}, std::pair{17, "CAC"}})
+                        for (const gemmi::Residue &r : st.models[0].chains[0].residues)
+                            if (r.name == "CYS" && r.seqid.num == cys)
+                                CHECK(bonded(topo.geometry, index_of(r, "SG"), index_of(*heme, carbon)));
+            }
+        }
+    }
+}
+
 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);