Merge branch 'af-mr3' into rc175: --af-model searches several copies per asymmetric unit (Matthews probability, sequential translation function, fixed-copy rigid body, tNCS)

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>

docs/ACKNOWLEDGEMENT.md

@@ -465,7 +465,15 @@ crystallographic software" (2007), J. Appl. Cryst. 40, 658-674
 [doi:10.1107/S0021889807021206](https://doi.org/10.1107/S0021889807021206). Whether one copy fits
 the cell is read off the Matthews coefficient: B. W. Matthews, "Solvent content of protein
 crystals" (1968), J. Mol. Biol. 33, 491-497
-[doi:10.1016/0022-2836(68)90205-2](https://doi.org/10.1016/0022-2836(68)90205-2).
+[doi:10.1016/0022-2836(68)90205-2](https://doi.org/10.1016/0022-2836(68)90205-2). How many copies
+to search for is weighed by a distribution of the Matthews coefficient over protein crystals made from
+the mode and mean of Kantardjieff and Rupp's: K. A. Kantardjieff and B. Rupp, "Matthews coefficient
+probabilities: Improved estimates for unit cell contents of proteins, DNA, and protein-nucleic acid
+complex crystals" (2003), Protein Sci. 12, 1865-1871
+[doi:10.1110/ps.0350503](https://doi.org/10.1110/ps.0350503). A copy refined beside the others has the
+others' placements read off by least-squares superposition: W. Kabsch, "A solution for the best rotation
+to relate two sets of vectors" (1976), Acta Cryst. A32, 922-923
+[doi:10.1107/S0567739476001873](https://doi.org/10.1107/S0567739476001873).
 
 ## Software and computing methods
 

docs/RUGNUX_REPORT.md

@@ -974,22 +974,35 @@ nothing else.
 
 `--af-model <file>` takes a predicted model - an AlphaFold mmCIF, with the pLDDT in the B-factor
 column - and places it by molecular replacement: residues below pLDDT 70 removed, B-factors from the
-pLDDT, then a rotation function over the point group (6 deg sampling, 4-12 A), a translation function
-for the 50 best orientations in every space group the merge cannot tell apart, and the best three
-placements per group through the rigid body, bulk solvent and R-free (to `SCORE_D_MIN`, 3 A or the
-data's limit). One copy in the asymmetric unit only, for now.
+pLDDT, then a rotation function over the point group (6 deg sampling, 4-12 A, each of the 50 best
+peaks refined off the sampling), a translation function for those orientations in every space group the
+merge cannot tell apart, and the best three placements per group that pack through the rigid body, bulk
+solvent and R-free (to `SCORE_D_MIN`, 3 A or the data's limit). Further copies follow one at a time, up
+to the largest copy number with a Matthews probability of 1% (`COPIES_LIKELY`: the most probable): the
+translation function again with the copies already placed as a known partial structure, its peaks
+ranked by the correlation of all the copies with the data, each new copy refined beside the others held
+fixed, and kept where its `TF_Z` reaches 8, it packs and R-free drops by 0.01 (`COPIES`, and
+`MODEL_PRIOR_COPIES=` for the group chosen). The groups are searched side by side; every group the best
+one leads is searched once more with the best one's copies' orientations tried first, and the better
+search kept.
 
-The table has one row per group and merge. A row is `PLACED` where its translation-function Z-score
-(`TF_Z`) reaches 8, at most 5% of its C-alpha atoms clash with symmetry mates (`CLASH`), one copy fits
-the cell (`MATTHEWS`: Vm below 1.5 A^3/Da; `MULTI_COPY`: Vm of 4 or more, which is reported but not
-supported yet). The null - the same model refitted from random orientations, whose R-work the
-placement's has to beat by 3 sd (`NULL_SIGMA`) - is built for the candidate the choice falls on; one
-that does not beat it is marked `NULL` and the choice falls on the next. A placement clashing with its
-symmetry mates over twice the limit is not refined. The rotation and translation functions run on the
+The table has one row per group and merge. A row is `PLACED` where one of its copies' translation-
+function Z-scores (`TF_Z`, the highest) reaches 8, at most 5% of the C-alpha atoms clash with another
+copy or a symmetry mate (`CLASH`), and one copy fits the cell (`MATTHEWS`: Vm below 1.5 A^3/Da). The
+null - the same model refitted from random orientations, each placed by the translation function as the
+search placed the real one (all copies turned together about their centroid where there are several) -
+is built for the candidate the choice falls on; its R-work has to beat the null's by 4.51 sd
+(`NULL_SIGMA`: 3 sd allowing for the null being nine draws, as `MODEL_FIT_SIGMA`), and one that does not
+is marked `NULL` and the choice falls on the next. A placement clashing over twice the limit is not
+refined. Translational NCS is looked for in the native Patterson (a peak off the lattice of a fifth of
+the origin's): with it, a group whose screw axes the tNCS translation turns into the chosen group's and
+whose search placed fewer copies leaves the choice `UNDECIDED`, since one copy beside a translated twin
+of itself fits the data no better than at random. The rotation and translation functions, the rigid
+body and the structure factors the R-factors are taken on run on the
 GPU where there is one (`MODEL_PRIOR_SEARCH_DEVICE=`). `MODEL_PRIOR_DECISION=` is then:
 
 - `AGREES` - the placed group with the lowest R-free is the answer's, ahead of every other group of
-  the merge that packs by 0.03 in R-free;
+  the merge that packs by 0.03 in R-free (with tNCS, see above);
 - `ENANTIOMORPH` - it is the answer's enantiomorph, and the answer's files take that label (a label
   only, as with `--model`; no reflection moves). It is taken only by the rule a `--model`'s
   enantiomorph is taken by: the model has to fit against its null (the choice already required that)

image_analysis/structure_refinement/MolecularReplacement.h

@@ -3,7 +3,7 @@
 
 #pragma once
 
-// Molecular replacement of ONE copy of a predicted model (rugnux --af-model, experimental): a global
+// Molecular replacement of a predicted model (rugnux --af-model, experimental): a global
 // search for the model's orientation and position in the cell, so that a model with no placement of
 // its own - an AlphaFold prediction - can be scored against the data, and the space groups the data
 // cannot tell apart can be scored against it.
@@ -13,11 +13,13 @@
 //    correlation of the observed E^2 with the model's own symmetry-summed |F|^2 - read off the
 //    model's molecular transform, computed once by FFT and interpolated per orientation.
 //  - Translation: the full-symmetry translation function on one FFT per orientation and space group.
+//  - Further copies, one at a time, as many as the cell can hold: the translation function again with
+//    the copies already placed as a known partial structure, each new copy refined against the others.
 //  - Rescoring: the best placements through the rigid body, bulk solvent, R-free, a packing check
 //    against the symmetry mates and a null of the same model in random orientations.
 //
-// Everything up to the rescoring runs on the CPU, in both builds; the rigid body runs on the GPU
-// where there is one.
+// The rotation and translation functions, the rigid body and the structure factors the R-factors are
+// taken on run on the GPU where there is one, and on the CPU otherwise; the same functions either way.
 
 #include <complex>
 #include <memory>
@@ -114,14 +116,32 @@ struct TranslationSetup {
     std::vector<gemmi::Vec3> s;                     // terms x ops
     std::vector<std::complex<double>> phase;        // terms x ops
     std::vector<double> weight;                     // terms
+    std::vector<double> e2, multiplicity;           // terms
     std::vector<int> bin;                           // terms
     int bins = 0;
     std::vector<int> pair_j, pair_k;                // the ordered pairs
     std::vector<int> target_u, target_v, target_l;  // terms x pairs
+    // With a known partial structure (MakeTranslationSetup's `known_part`): where the cross terms
+    // conj(F_known) g_j go, at h R_j - or, as their conjugate F_known conj(g_j), at -h R_j where that
+    // has l >= 0; both where l = 0. One entry per coefficient: its term, operator, index and whether it
+    // is the conjugate.
+    std::vector<int> cross_term, cross_op, cross_u, cross_v, cross_l;
+    std::vector<char> cross_conj;
     size_t Terms() const { return weight.size(); }
 };
+// `known_part`: the grid is sized for the cross terms with a known partial structure as well, which
+// reach out to h R_j where the copies' own terms reach h R_j - h R_k (in a polar group no further than
+// the plane the polar axis is free in).
 TranslationSetup MakeTranslationSetup(const SearchReflections &refl, const gemmi::UnitCell &cell,
-                                      const gemmi::SpaceGroup &sg);
+                                      const gemmi::SpaceGroup &sg, bool known_part = false);
+
+// F of copies of the model already placed - copy c turned by rotations[c] about its centroid, the
+// centroid at fractional positions[c] - at each term of `setup`, summed over the copies and their
+// symmetry mates, in the units the translation function's g_j are in.
+std::vector<std::complex<double>> KnownStructureFactors(const MolecularTransform &ft, const TranslationSetup &setup,
+                                                        const gemmi::UnitCell &cell,
+                                                        const std::vector<gemmi::Mat33> &rotations,
+                                                        const std::vector<gemmi::Fractional> &positions);
 
 // The translation function of the model in orientation `rotation`: its highest point as a Z-score over
 // every point of the function, and where it is (the fractional position of the model's centroid).
@@ -129,20 +149,60 @@ struct TranslationPeak {
     double z = 0;
     gemmi::Fractional t;
 };
+// With `known` (KnownStructureFactors on a setup made with known_part), the function of one more copy
+// beside the known ones: the overlap of E^2 - 1 with |F_known + F(t)|^2, whose cross term
+// 2 Re(conj(F_known) F(t)) is a Fourier series in t too and goes on the same grid.
 TranslationPeak TranslationFunction(const MolecularTransform &ft, const TranslationSetup &setup,
-                                    const gemmi::Mat33 &rotation);
+                                    const gemmi::Mat33 &rotation,
+                                    const std::vector<std::complex<double>> *known = nullptr);
+
+// The correlation of E^2 with |F|^2 of the known copies (KnownStructureFactors) and one more copy turned by
+// `rotation` with its centroid at `t`, over the terms of `setup`, each weighted by its multiplicity and each
+// bin's |F|^2 over its mean: the score the translation function's peaks for a further copy are ranked by.
+// Unlike the translation function, which is linear in |F|^2, it does not reward a copy that only repeats
+// the known ones - their own copy moved by one of the group's origin shifts doubles F on half the
+// reflections and cancels it on the other half, which the overlap with E^2 counts as a gain.
+double AssemblyCorrelation(const MolecularTransform &ft, const TranslationSetup &setup, const gemmi::UnitCell &cell,
+                           const std::vector<std::complex<double>> &known, const gemmi::Mat33 &rotation,
+                           const gemmi::Fractional &t);
 
 // The model placed: rotated about its centroid by `rotation`, the centroid at `t`, in `cell` and `sg`.
 gemmi::Structure PlaceSearchModel(const SearchModel &model, const gemmi::Mat33 &rotation, const gemmi::Fractional &t,
                                   const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg);
 
-// The fraction of the model's C-alpha atoms within 3 A of a C-alpha of a symmetry mate.
-double ClashFraction(const gemmi::Structure &placed);
+// The fraction of the C-alpha atoms within 3 A of a C-alpha of another copy: a symmetry mate, or one of
+// the other `copies` the structure's chains are split into evenly (in order), or their mates.
+double ClashFraction(const gemmi::Structure &placed, size_t copies = 1);
+
+// Copies of the model as one structure: their chains in order, renamed so that no two share a name.
+gemmi::Structure AssembleCopies(const std::vector<gemmi::Structure> &copies);
+
+// The probabilities of 1, 2, ... copies of a model of `mass_da` in the asymmetric unit of a cell of
+// `asu_volume` A^3, from the distribution of the Matthews coefficient over protein crystals; one entry
+// per copy number whose Vm is at least 1.5 A^3/Da.
+std::vector<double> MatthewsProbabilities(double asu_volume, double mass_da);
 
 // The space groups with the rotations and the centring of `sg`, in its setting: those merged
 // intensities cannot tell apart, apart from systematic absences. `sg` itself first.
 std::vector<gemmi::SpaceGroup> SameRotationGroups(const gemmi::SpaceGroup &sg);
 
+// Translational non-crystallographic symmetry: copies in the same orientation, a translation apart, make
+// the native Patterson peak off its origin where no lattice vector is. The highest peak of the Patterson
+// of E^2 (5-10 A), at least 15 A from the origin and from every centring vector, as a fraction of the
+// origin's height, and where it is; tNCS where that is a fifth or more.
+struct TranslationalNcs {
+    bool present = false;
+    double height = 0;              // over the origin's
+    gemmi::Fractional vector;
+};
+TranslationalNcs DetectTranslationalNcs(const std::vector<MergedReflection> &merged, const gemmi::UnitCell &cell,
+                                        const gemmi::SpaceGroup &sg);
+
+// Whether groups `a` and `b` (the same rotations) differ only in screw translations that the tNCS vector
+// `v` makes up: then copies a translation v apart put one group's operators onto the other's, and neither
+// is told from the other by a search that has not placed all of them.
+bool TncsRelatesGroups(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b, const gemmi::Fractional &v);
+
 // One merge the model is searched against: the answer's, or the twin hypothesis' in a subgroup.
 struct MRMerge {
     std::string name;                                // "MERGE" or "TWIN_SUBGROUP_MERGE"
@@ -151,20 +211,23 @@ struct MRMerge {
     gemmi::SpaceGroup space_group;                   // the group it was merged in
 };
 
-// One space group's best single-copy placement.
+// One space group's best placement: its first copy and every further one the search added.
 struct MRCandidate {
     int merge = 0;                                   // index into the merges searched
     gemmi::SpaceGroup space_group;
     double vm = 0;                                   // Matthews coefficient of one copy
-    int copies_likely = 1;
-    double rf_z = 0;                                 // of the orientation, over the rotation function
-    double tf_z = 0;
-    double clash_fraction = 0;
+    int copies_likely = 1;                           // the most probable copy number, by Matthews
+    int copies_max = 1;                              // the most the search tried (Matthews probability >= 1%)
+    int copies = 0;                                  // placed
+    double rf_z = 0;                                 // of the first copy's orientation, over the rotation function
+    double tf_z = 0;                                 // the highest of the copies'
+    std::vector<double> copy_tf_z;                   // each copy's, the later ones beside the earlier
+    double clash_fraction = 0;                       // of all the copies
     double r_work = NAN, r_free = NAN;               // after the rigid body, to d_min_score
     double k_sol = NAN, b_sol = NAN;                 // the bulk solvent R was computed with
     double null_sigma = NAN;                         // R-work above its random-orientation null, in sd
     std::string verdict;                             // PLACED, or why not
-    std::shared_ptr<gemmi::Structure> placed;
+    std::shared_ptr<gemmi::Structure> placed;        // every copy
 };
 
 struct MolecularReplacementResult {
@@ -180,6 +243,7 @@ struct MolecularReplacementResult {
     int chosen = -1;
     std::string decision;
     bool gpu_search = false;   // the rotation and translation functions ran on the GPU
+    TranslationalNcs tncs;     // of the answer's merge
     // An ENANTIOMORPH decision is taken into the answer's label only as model validation takes a model's
     // enantiomorph (ValidateAgainstModel): the placed model has to fit against its null and the anomalous
     // map at its atoms must not be inverted. Filled by the caller, which runs that validation.
@@ -199,11 +263,20 @@ struct MolecularReplacementResult {
     double t_model = 0, t_transform = 0, t_rotation = 0, t_translation = 0, t_rescore = 0, t_null = 0;
 };
 
-// Search `merges` (the answer's first) for one copy of the predicted model at `model_path`, in every
-// space group with the rotations of each merge's group, and decide. A candidate is PLACED where its
-// translation function stands out (TF Z >= 8), it packs, one copy fits the cell, and its R-work beats
-// its own null. The model chooses among the answer merge's groups where one placed candidate leads
+// Search `merges` (the answer's first) for the predicted model at `model_path`, in every space group with
+// the rotations of each merge's group - one copy, then further copies one at a time while each new one
+// stands out, packs and lowers R-free - and decide. A candidate is PLACED where a copy's translation
+// function stands out (TF Z >= 8), the copies pack, one copy fits the cell, and its R-work beats its own
+// null. The model chooses among the answer merge's groups where one placed candidate leads
 // every other by an R-free margin; it chooses the twin hypothesis' merge only where none of the
 // answer's groups can hold the model at all (Matthews, packing) and a group of the subgroup can.
 MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector<MRMerge> &merges,
                                                    Logger &logger, size_t nthreads);
+
+// The hand check of an ENANTIOMORPH decision, filled into result.hand_check: taken into the answer by
+// --model's rule for a model's enantiomorph - the model has to fit against its own null, which the choice
+// already required of it, and the anomalous map at its atoms (`merged`, the answer's merge in
+// `data_space_group`) must not come out inverted, read by model validation's own code on the placement
+// and the bulk solvent MR scored. Where it is not taken, the decision becomes ENANTIOMORPH_VETOED.
+void CheckEnantiomorphHand(MolecularReplacementResult &result, const std::vector<MergedReflection> &merged,
+                           const gemmi::SpaceGroup &data_space_group, Logger &logger, size_t nthreads);

image_analysis/structure_refinement/MolecularReplacementGPU.cu

@@ -194,6 +194,28 @@ __global__ void translation_gather(const int *offset, const int *entry, const fl
     coef[c] = sum;
 }
 
+// The cross terms with a known partial structure, added to the coefficients: w conj(F_known) g_j, or its
+// conjugate, for every entry landing on a grid point, in the order they were listed.
+__global__ void translation_gather_cross(const int *offset, const int *entry, const float2 *copies,
+                                         const float *scaled, const float2 *known, const int *cross_term,
+                                         const int *cross_op, const char *cross_conj, int nops, size_t cells,
+                                         float2 *coef) {
+    const size_t c = static_cast<size_t>(blockIdx.x) * blockDim.x + threadIdx.x;
+    if (c >= cells)
+        return;
+    float2 sum = coef[c];
+    for (int e = offset[c]; e < offset[c + 1]; e++) {
+        const int x = entry[e], t = cross_term[x];
+        const float2 k = known[t];
+        const float2 g = copies[static_cast<size_t>(t) * nops + cross_op[x]];
+        const float w = scaled[t];
+        const float re = k.x * g.x + k.y * g.y, im = k.x * g.y - k.y * g.x;   // conj(k) g
+        sum.x += w * re;
+        sum.y += w * (cross_conj[x] ? -im : im);
+    }
+    coef[c] = sum;
+}
+
 // Per block: the sum and the sum of squares of the function, and its highest point (the lowest index
 // on a tie).
 __global__ void translation_stats(const float *map, size_t n, double *sum, double *sum2, float *best,
@@ -329,12 +351,13 @@ std::vector<double> MolecularReplacementGPU::RotationFunction(const SearchReflec
     return out;
 }
 
-std::vector<TranslationPeak> MolecularReplacementGPU::TranslationFunction(const TranslationSetup &setup,
-                                                                          const std::vector<gemmi::Mat33> &rotations) {
+std::vector<TranslationPeak> MolecularReplacementGPU::TranslationFunction(
+        const TranslationSetup &setup, const std::vector<gemmi::Mat33> &rotations,
+        const std::vector<std::complex<double>> *known) {
     std::vector<TranslationPeak> out(rotations.size());
     const int terms = static_cast<int>(setup.Terms()), nops = setup.ops;
     const int npairs = static_cast<int>(setup.pair_j.size());
-    if (nops < 2 || terms == 0 || rotations.empty())
+    if (terms == 0 || rotations.empty() || (nops < 2 && known == nullptr))
         return out;
     cudaStream_t stream = state_->stream;
     const int nu = setup.nu, nv = setup.nv, nw = setup.nw, nwh = nw / 2 + 1;
@@ -356,6 +379,28 @@ std::vector<TranslationPeak> MolecularReplacementGPU::TranslationFunction(const
                 entry[fill[(static_cast<size_t>(setup.target_u[e]) * nv + setup.target_v[e]) * nwh
                            + setup.target_l[e]]++] = static_cast<int>(e);
     }
+    // The same for the cross terms.
+    std::vector<int> cross_offset(cells + 1, 0), cross_entry;
+    std::vector<float2> known_f;
+    if (known != nullptr) {
+        for (size_t e = 0; e < setup.cross_term.size(); e++)
+            cross_offset[(static_cast<size_t>(setup.cross_u[e]) * nv + setup.cross_v[e]) * nwh + setup.cross_l[e] + 1]++;
+        for (size_t c = 0; c < cells; c++)
+            cross_offset[c + 1] += cross_offset[c];
+        cross_entry.resize(cross_offset[cells]);
+        std::vector<int> fill(cross_offset.begin(), cross_offset.end() - 1);
+        for (size_t e = 0; e < setup.cross_term.size(); e++)
+            cross_entry[fill[(static_cast<size_t>(setup.cross_u[e]) * nv + setup.cross_v[e]) * nwh
+                             + setup.cross_l[e]]++] = static_cast<int>(e);
+        for (const std::complex<double> &k : *known)
+            known_f.push_back(make_float2(static_cast<float>(k.real()), static_cast<float>(k.imag())));
+    }
+    auto d_cross_offset = upload(cross_offset, stream);
+    auto d_cross_entry = upload(cross_entry, stream);
+    auto d_known = upload(known_f, stream);
+    auto d_cross_term = upload(setup.cross_term, stream);
+    auto d_cross_op = upload(setup.cross_op, stream);
+    auto d_cross_conj = upload(setup.cross_conj, stream);
     std::vector<float2> phase(setup.phase.size());
     for (size_t i = 0; i < phase.size(); i++)
         phase[i] = make_float2(static_cast<float>(setup.phase[i].real()), static_cast<float>(setup.phase[i].imag()));
@@ -391,6 +436,10 @@ std::vector<TranslationPeak> MolecularReplacementGPU::TranslationFunction(const
                                                            scaled);
             translation_gather<<<static_cast<unsigned>((cells + THREADS - 1) / THREADS), THREADS, 0, stream>>>(
                     d_offset, d_entry, copies, scaled, d_pair_j, d_pair_k, nops, npairs, cells, coef);
+            if (known != nullptr)
+                translation_gather_cross<<<static_cast<unsigned>((cells + THREADS - 1) / THREADS), THREADS, 0, stream>>>(
+                        d_cross_offset, d_cross_entry, copies, scaled, d_known, d_cross_term, d_cross_op, d_cross_conj,
+                        nops, cells, coef);
             cufft_err(cufftExecC2R(plan, reinterpret_cast<cufftComplex *>(coef.get()), map.get()));
             translation_stats<<<stat_blocks, THREADS, 0, stream>>>(
                     map, points, d_sum.get() + r * stat_blocks, d_sum2.get() + r * stat_blocks,

image_analysis/structure_refinement/MolecularReplacementGPU.h

@@ -28,9 +28,11 @@ public:
     std::vector<double> RotationFunction(const SearchReflections &refl, const gemmi::UnitCell &cell,
                                          const gemmi::SpaceGroup &sg, const std::vector<gemmi::Mat33> &rotations);
 
-    // The translation function of one group at each of `rotations`.
+    // The translation function of one group at each of `rotations`; with `known`, beside a known partial
+    // structure (KnownStructureFactors).
     std::vector<TranslationPeak> TranslationFunction(const TranslationSetup &setup,
-                                                     const std::vector<gemmi::Mat33> &rotations);
+                                                     const std::vector<gemmi::Mat33> &rotations,
+                                                     const std::vector<std::complex<double>> *known = nullptr);
 
 private:
     std::unique_ptr<MolecularReplacementGPUState> state_;

image_analysis/structure_refinement/RigidBodyGPU.cpp

@@ -280,6 +280,10 @@ void RigidBodyTargetGPU::Placement(const double q[6], double rotation[9], double
 }
 
 void RigidBodyTargetGPU::SetZone(const gemmi::AsuData<gemmi::ValueSigma<float>> &fobs, double d_min) {
+    // The body where this zone starts, for the mask it shares with fixed atoms: the last evaluation's.
+    std::array<double, 6> q_start{};
+    if (have_point_)
+        q_start = q_;
     d_min_ = d_min;
     solvent_fitted_ = false;
     host_scale_ = false;
@@ -287,6 +291,42 @@ void RigidBodyTargetGPU::SetZone(const gemmi::AsuData<gemmi::ValueSigma<float>>
     zone_ = pool_.ObservedZone(model_, cell_, sg_, fobs, d_min);
     try {
         engine_.SetZone(*zone_);
+        if (fixed != nullptr) {
+            // The fixed atoms are copies of the body (RefineRigidBody): each one's placement of the body's
+            // atoms, by the least-squares superposition of the two.
+            // Following Kabsch (1976) Acta Cryst. A32, 922-923
+            const std::vector<gemmi::Position> all = ModelPositions(*fixed);
+            const size_t n = base_.size();
+            std::vector<std::array<double, 12>> copies;
+            for (size_t c = 0; c + n <= all.size(); c += n) {
+                Eigen::Vector3d mean = Eigen::Vector3d::Zero();
+                for (size_t i = 0; i < n; i++)
+                    mean += Eigen::Vector3d(all[c + i].x, all[c + i].y, all[c + i].z);
+                mean /= static_cast<double>(n);
+                Eigen::Matrix3d h = Eigen::Matrix3d::Zero();
+                for (size_t i = 0; i < n; i++) {
+                    const Eigen::Vector3d p(base_[i].x - centre_.x, base_[i].y - centre_.y, base_[i].z - centre_.z);
+                    const Eigen::Vector3d q(all[c + i].x - mean.x(), all[c + i].y - mean.y(), all[c + i].z - mean.z());
+                    h += p * q.transpose();
+                }
+                const Eigen::JacobiSVD<Eigen::Matrix3d> svd(h, Eigen::ComputeFullU | Eigen::ComputeFullV);
+                Eigen::Matrix3d d = Eigen::Matrix3d::Identity();
+                d(2, 2) = (svd.matrixV() * svd.matrixU().transpose()).determinant() < 0 ? -1.0 : 1.0;
+                const Eigen::Matrix3d r = svd.matrixV() * d * svd.matrixU().transpose();
+                // The body's relative coordinates have their mean at zero, so the copy's centroid is t.
+                std::array<double, 12> x{};
+                for (int i = 0; i < 3; i++)
+                    for (int j = 0; j < 3; j++)
+                        x[3 * i + j] = r(i, j);
+                x[9] = mean.x();
+                x[10] = mean.y();
+                x[11] = mean.z();
+                copies.push_back(x);
+            }
+            std::array<double, 12> moving{};
+            Placement(q_start.data(), moving.data(), moving.data() + 9);
+            engine_.SetFixedCopies(copies, moving);
+        }
     } catch (const JFJochException &e) {
         throw RigidBodyGPUFailure(e.what());
     }
@@ -331,7 +371,7 @@ bool RigidBodyTargetGPU::Residuals(const double q[6], double *residuals) {
         if (zone_->row_hkl.empty())
             return false;
         engine_.Fcalc(rotation, translation);
-        if (!(hold_mask && have_point_))
+        if (fixed == nullptr && !(hold_mask && have_point_))
             engine_.Fmask();
 
         if (zone_->point_obs.empty())

image_analysis/structure_refinement/RigidBodyGPU.cu

@@ -106,7 +106,7 @@ __global__ void place_kernel(const double *rel, int n, Placement p, int nu, int
 // dF/dt where `df_dt` is not null. `f1` is cuFFT's forward transform, which carries exp(-2 pi i h.x):
 // conjugated and scaled by V/N it is gemmi's F1.
 __global__ void compose_kernel(const float2 *f1, float norm, const RigidBodyGPUTerm *terms, const double *row_scale,
-                               int rows, int ops, double2 *f, double2 *df_dt) {
+                               int rows, int ops, const double2 *fixed, double2 *f, double2 *df_dt) {
     const int m = blockIdx.x * blockDim.x + threadIdx.x;
     if (m >= rows)
         return;
@@ -126,7 +126,8 @@ __global__ void compose_kernel(const float2 *f1, float norm, const RigidBodyGPUT
         }
     }
     const double scale = row_scale[m];
-    f[m] = make_double2(scale * sr, scale * si);
+    f[m] = fixed != nullptr ? make_double2(scale * sr + fixed[m].x, scale * si + fixed[m].y)
+                            : make_double2(scale * sr, scale * si);
     if (df_dt != nullptr) {
         // times 2 pi i
         const double two_pi = 2 * 3.14159265358979323846;
@@ -135,6 +136,12 @@ __global__ void compose_kernel(const float2 *f1, float norm, const RigidBodyGPUT
     }
 }
 
+__global__ void add_kernel(const double2 *a, int n, double2 *sum) {
+    const int m = blockIdx.x * blockDim.x + threadIdx.x;
+    if (m < n)
+        sum[m] = make_double2(sum[m].x + a[m].x, sum[m].y + a[m].y);
+}
+
 // Fmask at each row, read at h directly (the mask is symmetric): gemmi's get_value_by_hkl() on the
 // transform of the mask, taking the Friedel mate where h is outside the stored half.
 __global__ void fmask_kernel(const float2 *fm, float norm, const int *row_hkl, int rows, int nu, int nv, int nw,
@@ -324,6 +331,8 @@ struct RigidBodyGPUEngineImpl {
     CudaDevicePtr<int> mask_slot;
     CudaDevicePtr<float> mask_radius;
     CudaDevicePtr<double> mask_atoms;   // 4 per mask atom: fractional x, y, z in [0,1), radius
+    CudaDevicePtr<double> copies_mask_atoms;   // the same for the body and fixed copies of it, grown as needed
+    size_t copies_mask_capacity = 0;
 
     CudaDevicePtr<float> grid;          // 3 grids back to back
     CudaDevicePtr<float2> spectrum;     // their 3 transforms
@@ -332,7 +341,8 @@ struct RigidBodyGPUEngineImpl {
     CudaDevicePtr<RigidBodyGPUTerm> terms;
     CudaDevicePtr<int> row_hkl;
     CudaDevicePtr<double> row_scale, row_stol2;
-    CudaDevicePtr<double2> f, df_dt, shifted;
+    CudaDevicePtr<double2> f, df_dt, shifted, fixed;
+    bool has_fixed = false;   // the zone's fixed Fcalc is in `fixed`
     CudaDevicePtr<float2> fmask;
 
     CudaDevicePtr<int> obs_row, point_obs;
@@ -375,6 +385,7 @@ struct RigidBodyGPUEngineImpl {
         f = CudaDevicePtr<double2>(nr, sync);
         df_dt = CudaDevicePtr<double2>(3 * nr, sync);
         shifted = CudaDevicePtr<double2>(3 * nr, sync);
+        fixed = CudaDevicePtr<double2>(nr, sync);
         fmask = CudaDevicePtr<float2>(nr, sync);
         obs_row = CudaDevicePtr<int>(no, sync);
         point_obs = CudaDevicePtr<int>(no, sync);
@@ -443,7 +454,7 @@ struct RigidBodyGPUEngineImpl {
     void Compose(int slot, double2 *out, double2 *derivative) {
         compose_kernel<<<blocks(rows), THREADS, 0, stream>>>(spectrum.get() + slot * ComplexPoints(), Norm(), terms,
                                                             row_scale, static_cast<int>(rows), static_cast<int>(ops),
-                                                            out, derivative);
+                                                            has_fixed ? fixed.get() : nullptr, out, derivative);
         cuda_err(cudaGetLastError());
     }
 
@@ -500,7 +511,7 @@ size_t RigidBodyGPUEngine::DeviceBytes(const RigidBodyGPUCapacity &c) {
     bytes += 3 * c.grid_points * sizeof(float) + 3 * c.complex_points * sizeof(float2) + c.fft_work_bytes;
     bytes += ModelDensityGPU::DeviceBytes(c.atoms, c.pairs, c.bricks);
     bytes += c.terms * sizeof(RigidBodyGPUTerm);
-    bytes += nr * (3 * sizeof(int) + 2 * sizeof(double) + 7 * sizeof(double2) + sizeof(float2));
+    bytes += nr * (3 * sizeof(int) + 2 * sizeof(double) + 8 * sizeof(double2) + sizeof(float2));
     bytes += no * (2 * sizeof(int) + sizeof(float) + (1 + 6 + MAX_SCALE_PARAMS) * sizeof(double) + 2 * sizeof(float2));
     bytes += ModelMaskGPU::DeviceBytes(c.grid_points) + ModelScaleGPU::DeviceBytes(no);
     return bytes;
@@ -559,6 +570,7 @@ void RigidBodyGPUEngine::SetZone(const RigidBodyGPUZone &zone) {
     e.n_constraints = zone.constraints.size();
     e.constraints = zone.constraints;
     e.f_mean = zone.f_mean;
+    e.has_fixed = false;
 
     std::vector<int> slot(zone.atoms.size(), -1);
     std::vector<float> radius(zone.atoms.size(), 0.0f);
@@ -607,6 +619,48 @@ void RigidBodyGPUEngine::Fcalc(const double rotation[9], const double translatio
     e.Compose(0, e.f, e.df_dt);
 }
 
+void RigidBodyGPUEngine::SetFixedCopies(const std::vector<std::array<double, 12>> &fixed,
+                                        const std::array<double, 12> &moving) {
+    RigidBodyGPUEngineImpl &e = *impl_;
+    cuda_err(cudaSetDevice(e.device));
+    const size_t copies = fixed.size() + 1;
+    if (e.copies_mask_capacity < copies * e.n_mask_atoms) {
+        e.copies_mask_capacity = copies * e.n_mask_atoms;
+        e.copies_mask_atoms = CudaDevicePtr<double>(4 * std::max<size_t>(e.copies_mask_capacity, 1), CudaAlloc::Synchronous);
+    }
+    cuda_err(cudaMemsetAsync(e.fixed, 0, e.rows * sizeof(double2), e.stream));
+    for (size_t c = 0; c < copies; c++) {
+        const std::array<double, 12> &x = c < fixed.size() ? fixed[c] : moving;
+        Placement p{};
+        std::copy(x.begin(), x.begin() + 9, p.r);
+        std::copy(x.begin() + 9, x.end(), p.t);
+        std::copy(e.frac, e.frac + 9, p.frac);
+        place_kernel<<<blocks(e.n_atoms), THREADS, 0, e.stream>>>(e.rel, e.n_atoms, p, e.geom.nu, e.geom.nv, e.geom.nw,
+                                                                e.mask_slot, e.mask_radius, e.pos,
+                                                                e.copies_mask_atoms.get() + 4 * c * e.n_mask_atoms);
+        cuda_err(cudaGetLastError());
+        if (c < fixed.size()) {
+            e.density->Compute(e.pos, e.grid.get());
+            cufft_err(cufftExecR2C(e.Plan(1), e.grid, reinterpret_cast<cufftComplex *>(e.spectrum.get())));
+            e.has_fixed = false;
+            e.Compose(0, e.shifted, nullptr);
+            add_kernel<<<blocks(e.rows), THREADS, 0, e.stream>>>(e.shifted, static_cast<int>(e.rows), e.fixed);
+            cuda_err(cudaGetLastError());
+        }
+    }
+    float *mask_grid = e.grid.get() + e.GridPoints();
+    cufftComplex *mask_spectrum = reinterpret_cast<cufftComplex *>(e.spectrum.get() + e.ComplexPoints());
+    e.mask->Compute(reinterpret_cast<const ModelMaskAtom *>(e.copies_mask_atoms.get()),
+                    static_cast<int>(copies * e.n_mask_atoms), mask_grid);
+    cufft_err(cufftExecR2C(e.Plan(1), mask_grid, mask_spectrum));
+    fmask_kernel<<<blocks(e.rows), THREADS, 0, e.stream>>>(reinterpret_cast<float2 *>(mask_spectrum), e.Norm(),
+                                                          e.row_hkl, static_cast<int>(e.rows), e.geom.nu, e.geom.nv,
+                                                          e.geom.nw, e.fmask);
+    cuda_err(cudaGetLastError());
+    cuda_err(cudaStreamSynchronize(e.stream));
+    e.has_fixed = true;
+}
+
 void RigidBodyGPUEngine::Fmask() {
     RigidBodyGPUEngineImpl &e = *impl_;
     cuda_err(cudaSetDevice(e.device));

image_analysis/structure_refinement/RigidBodyGPUEngine.h

@@ -109,6 +109,11 @@ public:
     // solvent mask of the same placement (ModelMaskGPU) and Fmask.
     void Fcalc(const double rotation[9], const double translation[3]);
     void Fmask();
+    // Per zone, after SetZone(): fixed atoms beside the body that are copies of it, each placed at (R, t) -
+    // rotation row-major, then translation, as Fcalc() takes them. Their Fcalc, composed one copy at a time,
+    // is added to the body's in every composition from now on, and the bulk-solvent mask of them and of the
+    // body at `moving` is held at the rows in place of Fmask()'s.
+    void SetFixedCopies(const std::vector<std::array<double, 12>> &fixed, const std::array<double, 12> &moving);
     // The scale at this evaluation's Fcalc and Fmask (ModelScaleGPU): FitModelScale()'s solvent grid, and
     // the overall scale and anisotropic B at a fixed solvent. FitSolvent() throws ModelScaleGPUTooFewReflections where
     // the grid's fits would chain (too few reflections for the isotropic start), which the host then does.

image_analysis/structure_refinement/RigidBodyRefine.cpp

@@ -274,6 +274,31 @@ void RigidBodyTarget::SetZone(const gemmi::AsuData<gemmi::ValueSigma<float>> &fo
     for (const auto &hv : fobs_.v)
         hkl.push_back(hv.hkl);
     composition_.emplace(grid, d_min, hkl);
+    if (fixed != nullptr)
+        fixed_part_ = RigidBodyFixed(model_, *fixed, cell_, sg_, d_min_, *composition_, nthreads_);
+}
+
+RigidBodyFixedPart RigidBodyFixed(const gemmi::Model &moving, const gemmi::Model &fixed, const gemmi::UnitCell &cell,
+                                  const gemmi::SpaceGroup &sg, double d_min, const SymmetryComposition &composition,
+                                  size_t nthreads) {
+    RigidBodyFixedPart out;
+    gemmi::DensityCalculator<Table, float> dc;
+    dc.d_min = d_min;
+    dc.rate = GRID_RATE;
+    dc.grid.unit_cell = cell;
+    dc.grid.spacegroup = &sg;
+    dc.set_refmac_compatible_blur(fixed);
+    PutModelDensityOnGrid(dc, fixed, {}, nthreads);
+    composition.Compose(MapToFPhi(dc.grid), dc.blur, out.fcalc, nullptr, nthreads);
+
+    gemmi::Model all = moving;
+    all.chains.insert(all.chains.end(), fixed.chains.begin(), fixed.chains.end());
+    gemmi::Grid<float> mask_grid = RigidBodyZoneGrid(cell, sg, d_min);
+    PutMaskOnGrid(mask_grid, all, OrbitLeaders(mask_grid, nthreads), nthreads);
+    const gemmi::FPhiGrid<float> fm = MapToFPhi(mask_grid);
+    for (const gemmi::Miller &h : composition.Hkl())
+        out.fmask.push_back(fm.get_value_by_hkl(h));
+    return out;
 }
 
 // The model's Fcalc at the zone's composed indices, from one copy of it on the zone's grid.
@@ -300,7 +325,11 @@ bool RigidBodyTarget::Residuals(const double q[6], double *residuals) {
 
     const std::vector<gemmi::Miller> &hkl = composition_->Hkl();
     std::vector<std::complex<float>> fmask;
-    if (hold_mask && have_point_) {
+    if (fixed != nullptr) {
+        for (size_t m = 0; m < fc.size(); m++)
+            fc[m] += fixed_part_.fcalc[m];
+        fmask = fixed_part_.fmask;
+    } else if (hold_mask && have_point_) {
         fmask = fmask_;
     } else {
         // The mask is symmetric as gridded, so it is read at h directly.
@@ -400,6 +429,9 @@ bool RigidBodyTarget::Jacobian(const double q[6], double *jacobian) {
             qj[j] += step;
             Place(qj, column_models_[j]);
             CopyFcalc(column_models_[j], shifted[j], nullptr);
+            if (fixed != nullptr)
+                for (size_t m = 0; m < shifted[j].size(); m++)
+                    shifted[j][m] += fixed_part_.fcalc[m];
         }));
     for (std::future<void> &f : running)
         f.get();
@@ -498,7 +530,8 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model,
                                       double d_min,
                                       Logger &logger,
                                       size_t nthreads,
-                                      RigidBodyGPUPool *gpu) {
+                                      RigidBodyGPUPool *gpu,
+                                      const gemmi::Model *fixed) {
     const auto t0 = std::chrono::steady_clock::now();
     RigidBodyRefineResult result;
     if (fobs.v.empty())
@@ -513,9 +546,12 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model,
     RigidBodyTargetBase &target = *target_backend;
     if (!target.Usable())
         return result;
+    target.fixed = fixed;
 
     const std::vector<double> ladder = RigidBodyLadder(d_min);
-    const gemmi::Mat33 gauge = RigidBodyGaugeProjector(sg, cell);
+    // With fixed atoms beside the body the origin is theirs, not free.
+    const gemmi::Mat33 gauge = fixed != nullptr ? gemmi::Mat33(0, 0, 0, 0, 0, 0, 0, 0, 0)
+                                                : RigidBodyGaugeProjector(sg, cell);
     double q[6] = {0, 0, 0, 0, 0, 0};
     bool any_zone_solved = false;
     for (double zone : ladder) {

image_analysis/structure_refinement/RigidBodyRefine.h

@@ -53,6 +53,11 @@ double RigidBodyReachDeg(const gemmi::Model &model, double d_min);
 //
 // `gpu`: evaluate the target on an engine of this pool (RigidBodyGPU.h, CUDA builds only) rather than
 // on the CPU. The same target either way; null is the CPU.
+//
+// `fixed`: other atoms of the crystal that do not move - the other copies of a molecular-replacement
+// search - whose Fcalc is added to the body's at every evaluation, and which the bulk-solvent mask covers
+// as well; that mask is computed once per zone, with the body where the zone starts, and held. With them
+// the origin is no longer free, and the translation is refined in all three directions.
 class RigidBodyGPUPool;
 RigidBodyRefineResult RefineRigidBody(gemmi::Model &model,
                                       const gemmi::UnitCell &cell,
@@ -61,7 +66,8 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model,
                                       double d_min,
                                       Logger &logger,
                                       size_t nthreads = 1,
-                                      RigidBodyGPUPool *gpu = nullptr);
+                                      RigidBodyGPUPool *gpu = nullptr,
+                                      const gemmi::Model *fixed = nullptr);
 
 // --- The pieces RefineRigidBody is built from, exposed for the tests and the GPU backend ---
 
@@ -123,6 +129,16 @@ private:
     std::vector<Term> terms_;
 };
 
+// What RefineRigidBody's fixed atoms add to a zone: their Fcalc at the zone's composed indices, and the
+// bulk-solvent mask of them and the body together (`moving`, where the zone starts) at the same indices.
+struct RigidBodyFixedPart {
+    std::vector<std::complex<double>> fcalc;
+    std::vector<std::complex<float>> fmask;
+};
+RigidBodyFixedPart RigidBodyFixed(const gemmi::Model &moving, const gemmi::Model &fixed, const gemmi::UnitCell &cell,
+                                  const gemmi::SpaceGroup &sg, double d_min, const SymmetryComposition &composition,
+                                  size_t nthreads);
+
 // The target RefineRigidBody minimises over one resolution zone: the amplitude residuals of the model
 // placed at q (the angle-axis rotation about the centroid times the rms radius, then the translation,
 // all in angstroms), with the scale re-fitted at every evaluation. The placement and the counters are
@@ -146,6 +162,7 @@ public:
     virtual bool Jacobian(const double q[6], double *jacobian) = 0;
 
     bool hold_mask = false;        // for the tests: keep the last evaluation's bulk-solvent mask
+    const gemmi::Model *fixed = nullptr;   // RefineRigidBody's fixed atoms; set before the first zone
     int evaluations = 0;
     int jacobians = 0;
     int unmatched = 0;             // zone observations with no calculated amplitude to compare against
@@ -186,6 +203,7 @@ private:
     double f_mean_ = 1;
     std::vector<size_t> orbit_leaders_;
     std::optional<SymmetryComposition> composition_;   // the zone's
+    RigidBodyFixedPart fixed_part_;                    // the zone's, with fixed atoms
     bool solvent_fitted_ = false;
 
     // The last evaluation, which the Jacobian is taken at.

rugnux/ResultReport.cpp

@@ -2591,7 +2591,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
             Add(s, KeyReal("SCORE_D_MIN", mp.d_min_score, "{:.2f}"));
             Add(s, KeyText("MODEL_PRIOR_DECISION", mp.decision));
             if (mp.chosen >= 0)
+            {
                 Add(s, KeyText("MODEL_PRIOR_SPACE_GROUP", mp.candidates[mp.chosen].space_group.xhm()));
+                Add(s, KeyInt("MODEL_PRIOR_COPIES", mp.candidates[mp.chosen].copies));
+            }
             if (mp.hand_check.performed) {
                 const auto &h = mp.hand_check;
                 Add(s, KeyEnum("MODEL_PRIOR_HAND_CHECK", h.adopted ? "ADOPTED" : "VETOED", {"ADOPTED", "VETOED"}));
@@ -2613,22 +2616,25 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
             Add(s, Blank());
             ReportEntry e;
             e.kind = ReportEntry::Kind::Table;
-            e.table.columns = {"MERGE", "SPACE_GROUP", "VM", "RF_Z", "TF_Z", "CLASH", "R_WORK", "R_FREE",
-                               "NULL_SIGMA", "VERDICT"};
-            e.table.text_header = "  MERGE                SPACE_GROUP     VM   RF_Z   TF_Z  CLASH  R_WORK  R_FREE  NULL_SIGMA  VERDICT\n"
-                                  "  -------------------  -----------  -----  -----  -----  -----  ------  ------  ----------  ----------";
+            e.table.columns = {"MERGE", "SPACE_GROUP", "VM", "COPIES", "COPIES_LIKELY", "RF_Z", "TF_Z", "CLASH",
+                               "R_WORK", "R_FREE", "NULL_SIGMA", "VERDICT"};
+            e.table.text_header = "  MERGE                SPACE_GROUP     VM  COPIES  LIKELY   RF_Z   TF_Z  CLASH  R_WORK  R_FREE  NULL_SIGMA  VERDICT\n"
+                                  "  -------------------  -----------  -----  ------  ------  -----  -----  -----  ------  ------  ----------  ----------";
             for (const auto &c : mp.candidates) {
                 const std::string merge = c.merge == 0 ? "MERGE" : "TWIN_SUBGROUP_MERGE";
                 const auto r = [](double v) { return std::isfinite(v) ? fmt::format("{:.3f}", v) : std::string("-"); };
                 e.table.text_rows.push_back(fmt::format(
-                        "  {:<19s}  {:<11s}  {:5.2f}  {:5.1f}  {:5.1f}  {:4.1f}%  {:>6s}  {:>6s}  {:>10s}  {}",
-                        merge, c.space_group.xhm(), c.vm, c.rf_z, c.tf_z, 100 * c.clash_fraction, r(c.r_work),
+                        "  {:<19s}  {:<11s}  {:5.2f}  {:6d}  {:6d}  {:5.1f}  {:5.1f}  {:4.1f}%  {:>6s}  {:>6s}  {:>10s}  {}",
+                        merge, c.space_group.xhm(), c.vm, c.copies, c.copies_likely, c.rf_z, c.tf_z,
+                        100 * c.clash_fraction, r(c.r_work),
                         r(c.r_free), std::isfinite(c.null_sigma) ? fmt::format("{:+.1f}", c.null_sigma) : std::string("-"),
                         c.verdict));
                 std::vector<ReportValue> row;
                 row.push_back(KeyText("", merge).value);
                 row.push_back(KeyText("", c.space_group.xhm()).value);
                 row.push_back(KeyReal("", c.vm, "{:.2f}").value);
+                row.push_back(KeyInt("", c.copies).value);
+                row.push_back(KeyInt("", c.copies_likely).value);
                 row.push_back(KeyReal("", c.rf_z, "{:.1f}").value);
                 row.push_back(KeyReal("", c.tf_z, "{:.1f}").value);
                 row.push_back(KeyReal("", c.clash_fraction, "{:.3f}").value);
@@ -2641,17 +2647,20 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
             Add(s, std::move(e));
             Add(s, Blank());
             Add(s, Prose(fmt::format(
-                    " One copy of the predicted model, trimmed to pLDDT >= 70, was searched for in every space group\n"
-                    " each merge cannot tell apart (rotation function on a 6 deg sampling, translation function,\n"
-                    " {:.1f}-{:.1f} A), and the best placements taken through the rigid body; R_WORK and R_FREE are\n"
-                    " to SCORE_D_MIN. PLACED needs TF_Z >= 8,\n"
-                    " at most 5% of the C-alpha atoms clashing with symmetry mates, one copy fitting the cell\n"
-                    " (MATTHEWS: Vm below 1.5; MULTI_COPY: Vm of 4 or more, not supported yet) and, for the\n"
-                    " candidate the choice falls on, an R-work at least 3 sd better than the same model refitted\n"
-                    " from random orientations (NULL where it is not; the next candidate is then taken). The model chooses\n"
-                    " within a merge where its placed group leads every other that packs by 0.03 in R-free, and\n"
-                    " the twin hypothesis' merge - searched only then - where no group of the answer's merge can\n"
-                    " hold it at all. A group one copy cannot fill is checked for packing only, not refined.\n"
+                    " The predicted model, trimmed to pLDDT >= 70, was searched for in every space group each merge\n"
+                    " cannot tell apart (rotation function on a 6 deg sampling, translation function, {:.1f}-{:.1f} A),\n"
+                    " and the best placements taken through the rigid body; then further copies, one at a time, up to\n"
+                    " the largest copy number with a Matthews probability of 1% (COPIES_LIKELY: the most probable),\n"
+                    " each with the translation function beside the copies already placed, kept where its TF_Z is 8\n"
+                    " or more, it packs and R-free drops by 0.01. R_WORK and R_FREE are of all COPIES, to\n"
+                    " SCORE_D_MIN; TF_Z is the highest of the copies'. PLACED needs TF_Z >= 8, at most 5% of the\n"
+                    " C-alpha atoms clashing with another copy or a symmetry mate, one copy fitting the cell\n"
+                    " (MATTHEWS: Vm below 1.5) and, for the candidate the choice falls on, an R-work at least 3 sd\n"
+                    " better than the same copies refitted from random orientations (NULL where it is not; the next\n"
+                    " candidate is then taken). The model chooses within a merge where its placed group leads every\n"
+                    " other that packs by 0.03 in R-free, and the twin hypothesis' merge - searched only then - where\n"
+                    " no group of the answer's merge can hold it at all. A group one copy cannot fit is checked for\n"
+                    " packing only, not refined.\n"
                     " The answer's files take the enantiomorph from it - only where model validation of the placed\n"
                     " model accepts it as it accepts a --model's: a fit over its null and an anomalous map at its\n"
                     " atoms that is not inverted (ENANTIOMORPH_VETOED otherwise) - and nothing else; another group it\n"

rugnux/RugnuxScaleMerge.cpp

@@ -4000,47 +4000,9 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
                                       "The predicted model was not placed: " + mr.failure_reason});
         if (mr.chosen >= 0) {
             const MRCandidate &c = mr.candidates[mr.chosen];
-            if (mr.decision == "ENANTIOMORPH") {
-                // Taken into the answer by --model's rule for a model's enantiomorph: the model has to fit
-                // against its own null - which the MR choice already required of it - and the anomalous map
-                // at its atoms must not come out inverted, read by model validation's own code on the
-                // placement and the bulk solvent MR scored.
-                const auto t_hand = std::chrono::steady_clock::now();
-                Logger quiet = Logger::Buffered();
-                const auto data_sg = experiment_.GetSpaceGroupOrP1();
-                ModelValidationResult v;
-                try {
-                    v = ReadModelAnomalousMap(sm.merged, *c.placed, data_sg, c.k_sol, c.b_sol, quiet,
-                                              static_cast<size_t>(config_.nthreads));
-                } catch (const std::runtime_error &e) {   // a GPU failure: no reading, so no relabelling
-                    logger.Warning("Model prior: the anomalous map could not be read ({})", e.what());
-                }
-                v.r_work_sigma = c.null_sigma;
-                v.model_fits = c.null_sigma >= 3.0;
-                v.adopted_model_enantiomorph = v.ok && v.model_fits && !v.anomalous_hands_disagree;
-                auto &h = mr.hand_check;
-                h.performed = true;
-                h.adopted = v.adopted_model_enantiomorph;
-                h.model_fits = v.model_fits;
-                h.fit_sigma = v.r_work_sigma;
-                h.anomalous_inverted = v.anomalous_hands_disagree;
-                h.anomalous_pairs = v.anomalous_pairs;
-                h.anomalous_scatterer_mean_sigma = v.anomalous_scatterer_mean_sigma;
-                h.deepest_site = v.anomalous_deepest_site;
-                h.deepest_sigma = v.anomalous_deepest_sigma;
-                h.seconds = std::chrono::duration<double>(std::chrono::steady_clock::now() - t_hand).count();
-                if (!h.adopted) {
-                    mr.decision = "ENANTIOMORPH_VETOED";
-                    logger.Info("Model prior: the enantiomorph {} is not taken into the answer - {} ({:.3f} s)",
-                                c.space_group.xhm(),
-                                h.anomalous_inverted ? fmt::format("the anomalous density at the placed model's "
-                                                                   "atoms is inverted ({} at {:.1f} sigma)",
-                                                                   h.deepest_site, h.deepest_sigma)
-                                                     : fmt::format("the placed model does not fit against its "
-                                                                   "null ({:+.1f} sd)", h.fit_sigma),
-                                h.seconds);
-                }
-            }
+            if (mr.decision == "ENANTIOMORPH")
+                CheckEnantiomorphHand(mr, sm.merged, experiment_.GetSpaceGroupOrP1(), logger,
+                                      static_cast<size_t>(config_.nthreads));
             if (mr.decision == "ENANTIOMORPH") {
                 // A change of label only (see ValidateAgainstModel): the two groups index identically.
                 experiment_.SetSpaceGroup(c.space_group);

tests/MolecularReplacementTest.cpp

@@ -11,9 +11,12 @@
 
 #include "../image_analysis/structure_refinement/ModelFFT.h"
 #include "../image_analysis/structure_refinement/MolecularReplacement.h"
+#include "../image_analysis/structure_refinement/RigidBodyRefine.h"
+#include "../common/Logger.h"
 #ifdef JFJOCH_USE_CUDA
 #include "../common/CUDAWrapper.h"
 #include "../image_analysis/structure_refinement/MolecularReplacementGPU.h"
+#include "../image_analysis/structure_refinement/RigidBodyGPU.h"
 #endif
 
 namespace {
@@ -39,7 +42,7 @@ namespace {
         c *= 1.0 / pos.size();
         for (const auto &p : pos) {
             gemmi::Atom a;
-            a.name = "C";
+            a.name = "CA";   // every atom counts in the clash check
             a.element = gemmi::El::C;
             a.pos = gemmi::Position(p - c);
             a.b_iso = 15.0f;
@@ -173,6 +176,114 @@ TEST_CASE("MolecularReplacement_SyntheticPlacement", "[MolecularReplacement]") {
     }
 }
 
+// Two copies in the asymmetric unit: with the first one known, the translation function of the second
+// at its true orientation peaks at its true position - not merely at one of the group's origins, which
+// only the cross term with the known copy tells apart.
+TEST_CASE("MolecularReplacement_SecondCopyBesideKnownFirst", "[MolecularReplacement]") {
+    const SearchModel model = synthetic_model();
+    const gemmi::UnitCell cell(52, 60, 70, 90, 90, 90);
+    const gemmi::SpaceGroup &sg = gemmi::get_spacegroup_by_name("P 21 21 21");
+    const gemmi::Mat33 r1 = axis_angle({1, 2, 3}, 50), r2 = axis_angle({-2, 1, 1}, 110);
+    const gemmi::Fractional t1(0.17, 0.31, 0.09), t2(0.41, 0.12, 0.33);
+    const gemmi::Structure crystal = AssembleCopies({PlaceSearchModel(model, r1, t1, cell, sg),
+                                                     PlaceSearchModel(model, r2, t2, cell, sg)});
+    const auto merged = observed(crystal);
+    const SearchReflections refl = MakeSearchReflections(merged, cell, sg, 4.0, 12.0);
+    const MolecularTransform ft(model, 4.0, 4);
+    const TranslationSetup setup = MakeTranslationSetup(refl, cell, sg, true);
+    const auto known = KnownStructureFactors(ft, setup, cell, {r1}, {t1});
+    const TranslationPeak peak = TranslationFunction(ft, setup, r2, &known);
+    CHECK(peak.z > 10.0);
+    for (double d : {peak.t.x - t2.x, peak.t.y - t2.y, peak.t.z - t2.z})
+        CHECK(std::fabs(d - std::round(d)) < 0.04);
+    // Two copies of the same body in one cell do not clash where they are placed apart.
+    CHECK(ClashFraction(crystal, 2) < 0.05);
+    CHECK(ClashFraction(AssembleCopies({PlaceSearchModel(model, r1, t1, cell, sg),
+                                        PlaceSearchModel(model, r2, t1, cell, sg)}), 2) > 0.5);
+}
+
+TEST_CASE("MolecularReplacement_MatthewsProbabilities", "[MolecularReplacement]") {
+    // Vm 2.5 A^3/Da for one copy: one copy is by far the most probable, four impossible (Vm 0.6).
+    const auto p = MatthewsProbabilities(2.5 * 20000, 20000);
+    REQUIRE(p.size() == 1);
+    CHECK(p[0] == Catch::Approx(1.0));
+    // Vm 9 for one copy: three or four copies (Vm 3, 2.25) are the likely ones, one copy is not.
+    const auto q = MatthewsProbabilities(9.0 * 20000, 20000);
+    REQUIRE(q.size() == 6);
+    CHECK(q[0] < 0.01);
+    CHECK(std::max_element(q.begin(), q.end()) - q.begin() == 3);
+}
+
+TEST_CASE("MolecularReplacement_TranslationalNcs", "[MolecularReplacement]") {
+    const SearchModel model = synthetic_model();
+    const gemmi::UnitCell cell(52, 60, 90, 90, 90, 90);
+    const gemmi::SpaceGroup &sg = gemmi::get_spacegroup_by_name("P 21 21 21");
+    const gemmi::Mat33 r = axis_angle({1, 2, 3}, 50);
+    const gemmi::Fractional t(0.17, 0.31, 0.09), v(0.0, 0.0, 0.5);
+    const auto one = observed(PlaceSearchModel(model, r, t, cell, sg));
+    CHECK_FALSE(DetectTranslationalNcs(one, cell, sg).present);
+    const auto two = observed(AssembleCopies({PlaceSearchModel(model, r, t, cell, sg),
+                                              PlaceSearchModel(model, r, gemmi::Fractional(t.x, t.y, t.z + v.z), cell, sg)}));
+    const TranslationalNcs tncs = DetectTranslationalNcs(two, cell, sg);
+    REQUIRE(tncs.present);
+    CHECK(std::fabs(tncs.vector.x - std::round(tncs.vector.x)) < 0.03);
+    CHECK(std::fabs(tncs.vector.y - std::round(tncs.vector.y)) < 0.03);
+    CHECK(std::fabs(tncs.vector.z - 0.5) < 0.03);
+}
+
+TEST_CASE("MolecularReplacement_TncsRelatesGroups", "[MolecularReplacement]") {
+    const auto &p3121 = gemmi::get_spacegroup_by_name("P 31 2 1"), &p3221 = gemmi::get_spacegroup_by_name("P 32 2 1");
+    // A translation of c/3 turns a 3(1) axis into a 3(2) one, wherever it lies in the plane.
+    CHECK(TncsRelatesGroups(p3121, p3221, gemmi::Fractional(0.3, 0.1, 1.0 / 3)));
+    CHECK_FALSE(TncsRelatesGroups(p3121, p3221, gemmi::Fractional(0.3, 0.1, 0.1)));
+    // A translation of a/2 turns 2(1) along a into 2, and nothing else.
+    const auto &p212121 = gemmi::get_spacegroup_by_name("P 21 21 21");
+    CHECK(TncsRelatesGroups(gemmi::get_spacegroup_by_name("P 2 21 21"), p212121, gemmi::Fractional(0.5, 0, 0)));
+    CHECK_FALSE(TncsRelatesGroups(gemmi::get_spacegroup_by_name("P 21 2 2"), p212121, gemmi::Fractional(0.5, 0, 0)));
+}
+
+// A copy refined beside another one held fixed goes back where it belongs, on the CPU and on the GPU alike.
+TEST_CASE("MolecularReplacement_RigidBodyBesideFixedCopy", "[MolecularReplacement]") {
+    const SearchModel model = synthetic_model();
+    const gemmi::UnitCell cell(52, 60, 70, 90, 90, 90);
+    const gemmi::SpaceGroup &sg = gemmi::get_spacegroup_by_name("P 21 21 21");
+    const gemmi::Mat33 r1 = axis_angle({1, 2, 3}, 50), r2 = axis_angle({-2, 1, 1}, 110);
+    const gemmi::Fractional t1(0.17, 0.31, 0.09), t2(0.41, 0.12, 0.33);
+    const gemmi::Structure first = PlaceSearchModel(model, r1, t1, cell, sg);
+    const gemmi::Structure truth = PlaceSearchModel(model, r2, t2, cell, sg);
+    const auto merged = observed(AssembleCopies({first, truth}));
+    gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
+    fobs.unit_cell_ = cell;
+    fobs.spacegroup_ = &sg;
+    for (const MergedReflection &m : merged)
+        fobs.v.push_back({{{m.h, m.k, m.l}}, {m.F, 1.0f}});
+    fobs.ensure_asu();
+    fobs.ensure_sorted();
+    // 3 deg and 0.5 A out.
+    const gemmi::Structure start = PlaceSearchModel(model, axis_angle({0, 0, 1}, 3).multiply(r2),
+                                                    gemmi::Fractional(t2.x + 0.5 / 52, t2.y, t2.z), cell, sg);
+    const auto worst = [&](const gemmi::Structure &st) {
+        const auto a = ModelPositions(st.models[0]), b = ModelPositions(truth.models[0]);
+        double w = 0;
+        for (size_t i = 0; i < a.size(); i++)
+            w = std::max(w, a[i].dist(b[i]));
+        return w;
+    };
+    Logger logger("MolecularReplacement_RigidBodyBesideFixedCopy");
+    gemmi::Structure cpu = start;
+    RefineRigidBody(cpu.models[0], cell, sg, fobs, 3.0, logger, 4, nullptr, &first.models[0]);
+    CHECK(worst(start) > 1.0);
+    CHECK(worst(cpu) < 0.1);
+#ifdef JFJOCH_USE_CUDA
+    if (get_gpu_count() > 0) {
+        gemmi::Structure gpu = start;
+        auto pool = RigidBodyGPUPool::Create(gpu.models[0], cell, sg, 3.0, fobs.v.size(), 1, logger);
+        REQUIRE(pool);
+        RefineRigidBody(gpu.models[0], cell, sg, fobs, 3.0, logger, 4, pool.get(), &first.models[0]);
+        CHECK(worst(gpu) < 0.1);
+    }
+#endif
+}
 
 #ifdef JFJOCH_USE_CUDA
 TEST_CASE("MolecularReplacement_GPUMatchesCPU", "[MolecularReplacement]") {
@@ -207,5 +318,13 @@ TEST_CASE("MolecularReplacement_GPUMatchesCPU", "[MolecularReplacement]") {
             CHECK(tf_gpu[i].t.z == Catch::Approx(cpu.t.z));
         }
     }
+    // Beside a known copy.
+    const TranslationSetup setup2 = MakeTranslationSetup(refl, cell, sg, true);
+    const auto known = KnownStructureFactors(ft, setup2, cell, {rotations[5]}, {gemmi::Fractional(0.3, 0.1, 0.4)});
+    const auto tf2_gpu = gpu.TranslationFunction(setup2, orientations, &known);
+    for (size_t i = 0; i < orientations.size(); i++) {
+        const TranslationPeak cpu = TranslationFunction(ft, setup2, orientations[i], &known);
+        CHECK(std::fabs(tf2_gpu[i].z - cpu.z) < 1e-3);
+    }
 }
 #endif

tools/CMakeLists.txt

@@ -12,6 +12,11 @@ TARGET_LINK_LIBRARIES(jfjoch_extract_hkl JFJochReader)
 
 # rugnux, the single offline analysis CLI, lives in rugnux/ next to its library.
 
+# rugnux --af-model's molecular replacement re-run on a run's merged MTZ files (developer utility, not
+# installed): a survey of the search over many datasets in minutes instead of a reprocessing each.
+ADD_EXECUTABLE(jfjoch_mr_replay jfjoch_mr_replay.cpp)
+TARGET_LINK_LIBRARIES(jfjoch_mr_replay JFJochStructureRefinement)
+
 # In-place re-compress /entry/data/data of a _data file from bitshuffle/LZ4 to bitshuffle/zstd.
 ADD_EXECUTABLE(jfjoch_recompress jfjoch_recompress.cpp)
 TARGET_LINK_LIBRARIES(jfjoch_recompress JFJochHDF5Wrappers)

tools/jfjoch_mr_replay.cpp

@@ -0,0 +1,112 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+// Developer utility: rugnux --af-model's molecular replacement re-run on the merged reflections a rugnux
+// run wrote, without reprocessing the images - the answer's p.mtz and, where the run made one, the twin
+// hypothesis' p_twin_<group>.mtz. The search, the decision and the hand check are rugnux's own code; the
+// output is the MODEL PRIOR keys and table of the rugnux report.
+//
+//   jfjoch_mr_replay <model> <answer group> <p.mtz> [<twin group> <p_twin.mtz>]
+//
+// The answer group is the one the run searched in: where it took the model's enantiomorph into its
+// label, the file carries the relabelled group and the original has to be given here.
+
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include <thread>
+
+#include <gemmi/mtz.hpp>
+
+#include "../common/Logger.h"
+#include "../image_analysis/structure_refinement/MolecularReplacement.h"
+
+namespace {
+
+std::vector<MergedReflection> read_merged(const std::string &path, UnitCell &cell) {
+    gemmi::Mtz mtz;
+    mtz.read_file(path);
+    const gemmi::UnitCell &c = mtz.cell;
+    cell = UnitCell(c.a, c.b, c.c, c.alpha, c.beta, c.gamma);
+    const auto col = [&](const char *label) -> int {
+        const gemmi::Mtz::Column *x = mtz.column_with_label(label);
+        return x ? static_cast<int>(x->idx) : -1;
+    };
+    const int i_mean = col("IMEAN"), s_mean = col("SIGIMEAN"), ip = col("I(+)"), sip = col("SIGI(+)"),
+              im = col("I(-)"), sim = col("SIGI(-)"), f = col("F"), sf = col("SIGF"), fp = col("F(+)"),
+              sfp = col("SIGF(+)"), fm = col("F(-)"), sfm = col("SIGF(-)"), free = col("FreeR_flag");
+    const size_t nc = mtz.columns.size();
+    const auto get = [&](size_t row, int i) { return i < 0 ? NAN : mtz.data[row * nc + i]; };
+    std::vector<MergedReflection> out;
+    for (size_t row = 0; row < static_cast<size_t>(mtz.nreflections); row++) {
+        MergedReflection r;
+        r.h = static_cast<int>(get(row, 0));
+        r.k = static_cast<int>(get(row, 1));
+        r.l = static_cast<int>(get(row, 2));
+        r.I = get(row, i_mean);
+        r.sigma = get(row, s_mean);
+        r.I_plus = get(row, ip);
+        r.sigma_plus = get(row, sip);
+        r.I_minus = get(row, im);
+        r.sigma_minus = get(row, sim);
+        r.F = get(row, f);
+        r.sigmaF = get(row, sf);
+        r.F_plus = get(row, fp);
+        r.sigmaF_plus = get(row, sfp);
+        r.F_minus = get(row, fm);
+        r.sigmaF_minus = get(row, sfm);
+        r.rfree_flag = get(row, free) == 0.0f;
+        r.d = static_cast<float>(1.0 / std::sqrt(c.calculate_1_d2(gemmi::Miller{{r.h, r.k, r.l}})));
+        out.push_back(r);
+    }
+    return out;
+}
+
+}  // namespace
+
+int main(int argc, char **argv) {
+    if (argc != 4 && argc != 6) {
+        std::cerr << "Usage: jfjoch_mr_replay <model> <answer group> <p.mtz> [<twin group> <p_twin.mtz>]\n";
+        return 1;
+    }
+    Logger logger("jfjoch_mr_replay");
+    const gemmi::SpaceGroup *answer_sg = gemmi::find_spacegroup_by_name(argv[2]);
+    UnitCell cell{}, twin_cell{};
+    const std::vector<MergedReflection> merged = read_merged(argv[3], cell);
+    std::vector<MergedReflection> twin;
+    std::vector<MRMerge> merges{{"MERGE", &merged, cell, *answer_sg}};
+    if (argc == 6) {
+        twin = read_merged(argv[5], twin_cell);
+        merges.push_back({"TWIN_SUBGROUP_MERGE", &twin, twin_cell, *gemmi::find_spacegroup_by_name(argv[4])});
+    }
+    const size_t nthreads = std::max(1u, std::thread::hardware_concurrency());
+    MolecularReplacementResult mr = RunMolecularReplacement(argv[1], merges, logger, nthreads);
+    if (mr.chosen >= 0 && mr.decision == "ENANTIOMORPH")
+        CheckEnantiomorphHand(mr, merged, *answer_sg, logger, nthreads);
+
+    std::cout << "MODEL_PRIOR= " << (mr.ok ? "PERFORMED" : "NOT_PERFORMED") << "\n";
+    if (!mr.ok) {
+        std::cout << "MODEL_PRIOR_REASON= " << mr.failure_reason << "\n";
+        return 0;
+    }
+    std::cout << "MODEL_PRIOR_DECISION= " << mr.decision << "\n";
+    if (mr.chosen >= 0) {
+        std::cout << "MODEL_PRIOR_SPACE_GROUP= " << mr.candidates[mr.chosen].space_group.xhm() << "\n";
+        std::cout << "MODEL_PRIOR_COPIES= " << mr.candidates[mr.chosen].copies << "\n";
+    }
+    if (mr.hand_check.performed)
+        std::cout << "MODEL_PRIOR_HAND_CHECK= " << (mr.hand_check.adopted ? "ADOPTED" : "VETOED") << "\n";
+    std::cout << fmt::format("MODEL_PRIOR_SECONDS= {:.2f}\n", mr.t_model + mr.t_transform + mr.t_rotation
+                                                                   + mr.t_translation + mr.t_rescore + mr.t_null
+                                                                   + mr.hand_check.seconds);
+    for (const MRCandidate &c : mr.candidates) {
+        const auto r = [](double v) { return std::isfinite(v) ? fmt::format("{:.3f}", v) : std::string("-"); };
+        std::cout << fmt::format(
+                "#  {:<19s}  {:<11s}  {:5.2f}  {:6d}  {:6d}  {:5.1f}  {:5.1f}  {:4.1f}%  {:>6s}  {:>6s}  {:>10s}  {}\n",
+                c.merge == 0 ? "MERGE" : "TWIN_SUBGROUP_MERGE", c.space_group.xhm(), c.vm, c.copies, c.copies_likely,
+                c.rf_z, c.tf_z,
+                100 * c.clash_fraction, r(c.r_work), r(c.r_free),
+                std::isfinite(c.null_sigma) ? fmt::format("{:+.1f}", c.null_sigma) : std::string("-"), c.verdict);
+    }
+    return 0;
+}