Merge refine-correct into rc181: free-set class from the structure factors, the as-given model in the final choice, the balanced-weight U trial on its own criterion

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37

# Conflicts:
#	image_analysis/structure_refinement/Refinement.h
#	rugnux/RugnuxScaleMerge.cpp

docs/CPU_DATA_ANALYSIS_REFINEMENT.md

@@ -313,8 +313,10 @@ are:
    as TLSMD cuts them (by dynamic programming over the least-squares fit of an isotropic TLS model, a quadratic B(r),
    to the refined B factors), each cut where it gives more groups than the one before; and a U per atom, at the U
    weight balanced as $w_{B,0}$ is (on the U shaken by 20% of its B) times $\{1/4, 1\}$. The cuts between one group
-   per chain and the finest run only where the finest fits the free set better (lower $-\ln L$) than one per chain;
-   on 49 refinements of 33 structures neither this nor leaving out the U weight $\times 4$ changed a choice taken. A
+   per chain and the finest run only where the finest fits the free set better (lower $-\ln L$) than one per chain,
+   and the U at its balanced weight only where the U at a quarter of it lowers R-free or $-\ln L$ of the free
+   reflections from where the weights ended; on 49 refinements of 33 structures neither the cuts' rule nor leaving
+   out the U weight $\times 4$ changed a choice taken. A
    trial displaces the one before it (the first: where the weights ended) only where its geometry is within bounds
    and both its R-free (beyond the tie) and the $-\ln L$ of the free reflections are lower: the free set tests the
    parameters a model adds, so the choice is made on the evidence of each run, not at a resolution (Merritt 2012).
@@ -322,6 +324,12 @@ are:
    `REFINE_ANISO_ATOMS`, `REFINE_TLS_GROUPS`). It is decided after the solvent build (§15.10), on the model with its
    solvent: a water a round put in refines a U of its own in the per-atom trial as every refinable atom does, and no
    TLS group holds a water.
+5. **The model as given** is the first candidate in order of trust: where the refinement does not displace it by the
+   rule every trial is held to (its geometry within bounds and R-free lower beyond the tie), the model is written as
+   it came, scaled to the data, and the solvent build's edits are undone (`REFINE_AS_GIVEN`). What comes out is never
+   worse on the free set than what went in. This holds only where the free set is free to the model as given: a model
+   placed as given on a free set drawn by the run (`MODEL_FREE_SET`; in `--mode refine`, a file without free flags) has
+   most likely been refined against those reflections, and its R-free as given does not compete.
 
 Where the model has copies (§15.5), the NCS restraints are one more weight: every rung of the first ladder is run
 without them, and then with them at $\sigma = 0.1$ Å at the B factors' weight the rungs without chose (NCS moves the

docs/RUGNUX_MODEL_AND_REFERENCE.md

@@ -339,7 +339,13 @@ columns are found, not named: amplitudes with their sigmas (F/SIGF, FP/SIGFP, FO
 map coefficients such as FC, FWT, DELFWT are passed over) where the file has them, else intensities (IMEAN/SIGIMEAN,
 I/SIGI, ...) made into amplitudes by French-Wilson, else the mean of the anomalous pairs. The free set is read from
 a FreeR_flag-like column (flag 0 is the test set; where 0 marks the majority, as phenix and CNS write it, the
-minority value is); a file without one gets rugnux's own 5 % set. The refinement is the one run after a merge on a placed model, and
+minority value is; in a structure-factor mmCIF, `_refln.status` 'f' first); a file without one gets rugnux's own 5 % set.
+Where a structure-factor mmCIF's flags are classes (`_refln.pdbx_r_free_flag` 0-19), the class the file marks is not
+always the one the depositors held out, so the model as given is scored on each class: the mark is kept where the model
+fits that class worst or the test-set count the model file records (REMARK 3, `_refine.ls_number_reflns_R_free`) is
+nearest its count; otherwise the class the model fits worst is the test set where that count is nearest it, or where
+it is above the marked class by more than twice the scatter of the other classes' R (`REFINE_FREE_SET_SOURCE=
+STATUS|REMARK3_COUNT|WORST_FIT`, `REFINE_FREE_SET_CLASS`). The refinement is the one run after a merge on a placed model, and
 writes the same files - `<prefix>_refined_model.cif`/`.pdb`, `<prefix>_refined_maps.mtz`, `<prefix>_blobs.json` - with
 the `REFINE_` keys in `<prefix>_report.txt`, beside `REFINE_FRENCH_WILSON=` and `REFINE_FREE_FLAGS= FILE|GENERATED`
 saying what had to be made of the input.

docs/RUGNUX_REPORT.md

@@ -1235,6 +1235,7 @@ Where a model was placed and refined ([Refinement](RUGNUX_MODEL_AND_REFERENCE.md
 | `REFINE_RIGID_BODIES` | The chains moved as rigid bodies, where that step was kept (0: not kept) |
 | `REFINE_R_WORK_PLACEMENT`, `REFINE_R_FREE_PLACEMENT` | After the placement stage (rigid bodies, then one B per residue, each kept only where it lowers R-free) |
 | `REFINE_R_WORK`, `REFINE_R_FREE` | After the restrained stage: the refined model, the one written as `<prefix>_refined_model.cif` |
+| `REFINE_AS_GIVEN` | `TRUE` where no refinement lowered R-free beyond the tie of the model as given, which is then the model written (as given, scaled to the data) |
 | `REFINE_RMSD_BONDS`, `REFINE_RMSD_ANGLES` | The refined model's rmsd from the dictionary's bonds (Å) and angles (°); the weight is chosen within 0.02 Å and 2° |
 | `REFINE_RMSD_METAL_BONDS` | The rmsd of the metal-donor bonds made from the model's LINK records (those within the metal's inner sphere), apart from the covalent bonds. Present where there are any |
 | `REFINE_RMSD_BONDS_START`, `REFINE_RMSD_ANGLES_START` | *dev*. The same for the model as given |
@@ -1307,4 +1308,7 @@ resolution shells. `rugnux --mode refine --refine-score-only` writes this score
 In `--mode refine` the input's keys come first: `REFINE_REFLECTIONS` and `REFINE_FREE_REFLECTIONS` (those read,
 after any `--refine-d-min`), `REFINE_FRENCH_WILSON` (`TRUE` where the file had intensities only and the amplitudes
 were made by French-Wilson), `REFINE_FREE_FLAGS` (`FILE`, or `GENERATED` where the file had no free set and
-rugnux drew its own) and `REFINE_D_MIN_LIMIT` (the `--refine-d-min` given).
+rugnux drew its own), `REFINE_FREE_SET_SOURCE` and `REFINE_FREE_SET_CLASS` (where a structure-factor mmCIF's
+flags are classes: the class taken as the test set, and whether the file's own mark (`STATUS`), the model file's
+test-set count with the worst fit (`REMARK3_COUNT`) or the worst fit alone (`WORST_FIT`) chose it) and
+`REFINE_D_MIN_LIMIT` (the `--refine-d-min` given).

image_analysis/structure_refinement/MolecularReplacement.cpp

@@ -2994,9 +2994,12 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec
                     copy_g = {{}};
                     copy_rg = {{}};
                 }
-                glog.Info("Molecular replacement:   {:<12} RF Z {:5.1f} TF Z {:5.1f} clash {:4.1f}%  R-work {:.3f} R-free {:.3f}",
+                // a placement that clashes too much to be one is not refined, and has no R-factors
+                glog.Info("Molecular replacement:   {:<12} RF Z {:5.1f} TF Z {:5.1f} clash {:4.1f}%  {}",
                             sg.xhm(), rf_sd > 0 ? (rfs[i] - rf_mean) / rf_sd : 0.0, tf[i].z, 100 * f.clash,
-                            score.r_work, score.r_free);
+                            std::isfinite(score.r_free)
+                                ? fmt::format("R-work {:.3f} R-free {:.3f}", score.r_work, score.r_free)
+                                : std::string("not refined"));
             }
             // The likelihood of the placement kept, judged against the same orientation elsewhere - only where
             // its translation function's Z-score does not already stand out.
@@ -3166,9 +3169,13 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec
                         trial.refined = LikelihoodTranslationZ(lsetup(), *lknown, known_llg, trial.g, centroid(trial.copy),
                                                                model.error_A, nthreads);
                     }
-                    glog.Info("Molecular replacement:   {:<12} copy {} TF Z {:5.1f} LLG gain {:7.1f}, refined {:7.1f} Z {:5.1f} clash {:4.1f}%  R-work {:.3f} R-free {:.3f}",
-                                sg.xhm(), n, tf2[trial.i].z, trial.gain, trial.refined.gain,
-                                trial.refined.z, 100 * trial.clash, trial.score.r_work, trial.score.r_free);
+                    // the likelihood's own Z only where the translation function's does not already stand out
+                    glog.Info("Molecular replacement:   {:<12} copy {} TF Z {:5.1f} LLG gain {:7.1f}, {} clash {:4.1f}%  R-work {:.3f} R-free {:.3f}",
+                                sg.xhm(), n, tf2[trial.i].z, trial.gain,
+                                std::isfinite(trial.refined.z)
+                                    ? fmt::format("refined {:7.1f} Z {:5.1f}", trial.refined.gain, trial.refined.z)
+                                    : std::string("likelihood Z not needed"),
+                                100 * trial.clash, trial.score.r_work, trial.score.r_free);
                     // The first, in that order, that stands out, packs and lowers R-free.
                     if ((tf2[trial.i].z >= TF_Z_PLACED || trial.refined.z >= TF_Z_PLACED) && trial.clash <= CLASH_MAX
                         && trial.score.r_free <= cand.r_free - COPY_R_FREE_GAIN) {
@@ -3176,8 +3183,12 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec
                         break;
                     }
                 }
-                if (!best)
+                if (!best) {
+                    glog.Info("Molecular replacement:   {:<12} copy {}: none of {} trials stands out, packs and lowers "
+                              "R-free by {:.3f}; {} {} kept", sg.xhm(), n, trials.size(), COPY_R_FREE_GAIN,
+                              cand.copies, cand.copies == 1 ? "copy" : "copies");
                     break;
+                }
                 copies.push_back(best->copy);
                 copy_rotation.push_back(best->rotation);
                 copy_position.push_back(centroid(best->copy));

image_analysis/structure_refinement/Refinement.cpp

@@ -2209,6 +2209,27 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
                 result.on_gpu ? "GPU" : "CPU", r.r_work, r.r_free, ctx.k_sol, ctx.b_sol, result.score_start.b_iso,
                 result.score_start.mask_spacing);
 
+    // The model as given, scaled: a candidate of the final choice beside the trials (the restrained stage, 5.), with
+    // its atoms and residues as they are, which the solvent build edits in place.
+    TrialOutcome given;
+    given.trial.r_work = r.r_work;
+    given.trial.r_free = r.r_free;
+    given.trial.nll_free = result.score_start.nll_free;
+    given.pos = ctx.Positions();
+    given.b = ctx.B();
+    given.aniso = ctx.aniso;
+    given.fit = ctx.Fitted();
+    given.ran = true;
+    std::vector<gemmi::Atom> given_atoms;
+    std::vector<std::pair<std::string, gemmi::EntityType>> given_residues;
+    const SolventPool given_pool = pool;
+    if (settings.restrained && !settings.score_only)
+        for (const gemmi::Chain &ch : ctx.Model().chains)
+            for (const gemmi::Residue &res : ch.residues) {
+                given_residues.emplace_back(res.name, res.entity_type);
+                given_atoms.insert(given_atoms.end(), res.atoms.begin(), res.atoms.end());
+            }
+
     // The waters and ions the build relabelled and kept (SolventEdit::site their old name), for the metal records the
     // refined model is written with.
     std::vector<SolventEdit> relabels;
@@ -2385,6 +2406,11 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
             RestraintSet repulsions = BuildRepulsions(topo, ctx.Model(), cell, sg, 0.0, ctx.nthreads);
             result.geometry_start = MeasureGeometry(topo.geometry, repulsions, ctx.XyzB(), &topo.metal_bonds);
         }
+        // (the placement moves rigid bodies, which leaves the bonds and angles as given)
+        given.trial.geometry = result.geometry_start;
+        given.trial.within_bounds = given.trial.geometry.rmsd_bonds <= MAX_RMSD_BONDS &&
+                                    given.trial.geometry.rmsd_angles <= MAX_RMSD_ANGLES;
+        given.trial.adp = result.adp_input;
         double w0 = 0, adp_w0 = 0;
         {
             const RefinementStage stage("balanced_weight");
@@ -2871,8 +2897,8 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
         //    The trials that would not change the choice are left out: the model's own TLS groups are tried only where
         //    one body per chain or the finest cut lowers R-free, the cuts between those only where also the finest
         //    fits the free set better, the U at its balanced weight only where the U at a quarter of it lowers R-free
-        //    (run ahead where the slots take them), and the U at four times its balanced weight not at all (on 49
-        //    refinements of 33 structures, no choice taken changed).
+        //    or the -ln L of the free reflections (run ahead where the slots take them), and the U at four times its
+        //    balanced weight not at all (on 49 refinements of 33 structures, no choice taken changed).
         if (settings.adp_hypotheses && state.ran && !ctx.Cancelled()) {
             const RefinementStage adp_stage("adp_hypotheses");
             ctx.SetPositions(state.pos);
@@ -2881,6 +2907,8 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
             ctx.Restore(state.fit);
             PlaceRiders(ctx, topo);
             state.trial.nll_free = ctx.Score().nll_free;
+            logger.Info("Refinement, displacements: from where the weights ended, R-free {:.4f}, -LL free {:.1f}",
+                        state.trial.r_free, state.trial.nll_free);
             std::vector<AdpHypothesis> hyps = AdpHypotheses(ctx, topo, input_tls);
             const double w = state.trial.weight, adp = state.trial.adp_weight > 0 ? state.trial.adp_weight : adp_w0;
             // the U weight balanced on the model with every refinable atom given a U
@@ -2956,13 +2984,19 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
             run(now);
             // the model's own groups where one body per chain or the finest cut lowers R-free from where the weights
             // ended, the cuts between where also the finest fits the free set better than one body per chain, the U
-            // at its balanced weight where the U at a quarter of it lowers R-free (on 364 refinements of the battery
-            // that changed the choice on 7, none by more than 0.003 in R-free, and left out 6 in 10 of these trials)
+            // at its balanced weight where the U at a quarter of it lowers R-free or -ln L free (R-free alone left out
+            // the trial that would have been kept where the quarter lowered -ln L free but not R-free)
             auto lowers = [&](int k) {
                 return k >= 0 && outcomes[k].ran && outcomes[k].trial.within_bounds &&
                        outcomes[k].trial.r_free < state.trial.r_free;
             };
-            const bool tls = lowers(coarsest) || lowers(finest), u_full = lowers(u_quarter);
+            // the U at a quarter fitting the free reflections better than where the weights ended asks for the
+            // balanced weight as well: the choice below is made on -ln L free as much as on R-free
+            auto fits_better = [&](int k) {
+                return k >= 0 && outcomes[k].ran && outcomes[k].trial.within_bounds &&
+                       outcomes[k].trial.nll_free < state.trial.nll_free;
+            };
+            const bool tls = lowers(coarsest) || lowers(finest), u_full = lowers(u_quarter) || fits_better(u_quarter);
             const bool cuts = tls && finest > coarsest && outcomes[finest].ran && outcomes[coarsest].ran &&
                               outcomes[finest].trial.nll_free < outcomes[coarsest].trial.nll_free;
             std::vector<char> needed(trials.size(), 0);
@@ -2992,6 +3026,28 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
             }
             logger.Info("Refinement: displacement model {}", result.adp_model);
         }
+        // 5. The model as given, first in order of trust: the refinement is kept where it displaces it by the rule every
+        //    trial is held to, else the model is given back as it came (the solvent build's edits undone), so that what
+        //    comes out is never worse on the free set than what went in - where the free set is free to the model as
+        //    given (RefinementSettings::as_given_competes).
+        if (settings.as_given_competes && state.ran && !ctx.Cancelled() && !displaces(state.trial, given.trial)) {
+            logger.Info("Refinement: the model as given kept (R-work {:.4f}, R-free {:.4f}), the refinement did not "
+                        "displace it (R-work {:.4f}, R-free {:.4f})", given.trial.r_work, given.trial.r_free,
+                        state.trial.r_work, state.trial.r_free);
+            size_t a = 0, k = 0;
+            for (gemmi::Chain &ch : ctx.Model().chains)
+                for (gemmi::Residue &res : ch.residues) {
+                    std::tie(res.name, res.entity_type) = given_residues[k++];
+                    for (gemmi::Atom &at : res.atoms)
+                        at = given_atoms[a++];
+                }
+            pool = given_pool;
+            relabels.clear();
+            result.solvent_build = SolventBuildResult();
+            result.adp_model = result.adp_input == "isotropic" ? "isotropic" : "the model's U held";
+            result.as_given = true;
+            state = given;
+        }
         for (size_t s = 1; s < trial_slots.size(); s++)
             ctx.evaluations += trial_slots[s].ctx->evaluations;
         ctx.evaluations -= ahead_evaluations;

image_analysis/structure_refinement/Refinement.h

@@ -78,6 +78,10 @@ struct RefinementSettings {
     // The data's own intensity statistics read twinned - the caller's verdict, from the merge (the L-test): the law
     // is then taken from 2 standard errors over its null instead of 3 (AdoptTwinLaw).
     bool twin_data_indicated = false;
+    // The model as given a candidate of the final choice (the restrained stage, 5.): only where its R-free is a free
+    // statistic - false where the model has most likely been refined against the free reflections too (a deposited
+    // model of this crystal scored on a free set drawn here).
+    bool as_given_competes = true;
     bool gpu = true;
     size_t nthreads = 1;
     size_t gpu_resident_bytes = 0;   // held on the card by what runs beside
@@ -159,6 +163,7 @@ struct RefinementResult {
     double r_work_placement = 1, r_free_placement = 1;
     // After the restrained stage (NaN where it did not run).
     bool restrained = false;
+    bool as_given = false;           // the model as given kept: no trial displaced it (the restrained stage ran)
     bool cancelled = false;          // cut short by RefinementSettings::cancel
     double r_work = NAN, r_free = NAN;
     double weight = 0;

rugnux/ResultReport.cpp

@@ -569,6 +569,7 @@ namespace {
             if (rf.restrained) {
                 Add(s, KeyReal("REFINE_R_WORK", rf.r_work, "{:.4f}"));
                 Add(s, KeyReal("REFINE_R_FREE", rf.r_free, "{:.4f}"));
+                Add(s, KeyBool("REFINE_AS_GIVEN", rf.as_given));
                 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}"));
@@ -3438,6 +3439,10 @@ ReportDocument BuildRefineReportDocument(const std::string &output_prefix, const
     Add(s, KeyInt("REFINE_FREE_REFLECTIONS", r.free_reflections));
     Add(s, KeyBool("REFINE_FRENCH_WILSON", r.french_wilson));
     Add(s, KeyEnum("REFINE_FREE_FLAGS", r.free_flags_generated ? "GENERATED" : "FILE", {"FILE", "GENERATED"}));
+    if (!r.free_set_source.empty()) {
+        Add(s, KeyEnum("REFINE_FREE_SET_SOURCE", r.free_set_source, {"STATUS", "REMARK3_COUNT", "WORST_FIT"}));
+        Add(s, KeyInt("REFINE_FREE_SET_CLASS", r.free_class));
+    }
     if (r.d_min_limit > 0)
         Add(s, KeyReal("REFINE_D_MIN_LIMIT", r.d_min_limit, "{:.3f}"));
     AddRefinementKeys(s, rf, output_prefix);

rugnux/RugnuxRefine.cpp

@@ -3,13 +3,134 @@
 
 #include "RugnuxRefine.h"
 
+#include <algorithm>
+#include <array>
 #include <cmath>
+#include <map>
+
+#include "gemmi/symmetry.hpp"
 
 #include "../common/UnitCell.h"
 #include "../image_analysis/scale_merge/FrenchWilson.h"
 #include "../image_analysis/scale_merge/RfreeFlags.h"
 #include "../image_analysis/structure_refinement/ModelValidation.h"   // SameIntensities
 
+namespace {
+
+// The test set of a structure-factor mmCIF whose free flags are classes (0-19): the file marks one (_refln.status
+// 'f', or flag 0), but that mark is the archive's reading of the flags and is not always the class the depositors'
+// refinement held out. The model as given tells: it fits worst the class it was not refined against, and the
+// test-set count its file records (REMARK 3 FREE R VALUE TEST SET COUNT, _refine.ls_number_reflns_R_free) is that
+// class's count within the refinement's resolution range. The count alone does not single a class out (of 20 classes
+// of about the same size another often has it too, or the marked one is off by a reflection or two), so the mark is
+// kept ("STATUS") where the model fits it worst or the count is nearest it. Else the class the model fits worst is
+// taken where the count is nearest it ("REMARK3_COUNT"), or where it is above the marked one by more than twice the
+// scatter of the R of the other classes, so that the marked class is fitted like a working one ("WORST_FIT"; the
+// model file gives no count, or the count points nowhere). `class_of`: the class of each reflection of the file; `counts`: per class, its reflections within that range.
+void ChooseFreeClass(std::vector<MergedReflection> &refl, const std::map<gemmi::Miller, int> &class_of,
+                     const std::map<int, int> &counts, double d_hi, double d_lo, const gemmi::Structure &model, const gemmi::UnitCell &cell,
+                     const gemmi::SpaceGroup &sg, double wavelength_A, const RefinementSettings &settings,
+                     Logger &logger, RefineOnlyResult &out) {
+    std::vector<int> cls(refl.size(), -1);
+    for (size_t i = 0; i < refl.size(); i++)
+        if (const auto c = class_of.find({{refl[i].h, refl[i].k, refl[i].l}}); c != class_of.end())
+            cls[i] = c->second;
+    // the class the file marks: that of most reflections flagged free
+    std::map<int, int> flagged;
+    for (size_t i = 0; i < refl.size(); i++)
+        if (refl[i].rfree_flag && cls[i] >= 0)
+            flagged[cls[i]]++;
+    int marked = -1;
+    for (const auto &[c, n] : flagged)
+        if (marked < 0 || n > flagged[marked])
+            marked = c;
+
+    // the model as given, scored against these flags: R per class
+    RefinementSettings score = settings;
+    score.score_only = true;
+    score.twin = false;
+    score.solvent = false;
+    logger.Info("Refinement input: the free flags are {} classes; the model as given scored on each", counts.size());
+    const RefinementResult s = RefineModel(model, refl, cell, sg, wavelength_A, logger, score);
+    if (!s.ok)
+        return;
+    const gemmi::GroupOps gops = sg.operations();
+    const gemmi::ReciprocalAsu asu(&sg);
+    auto in_asu = [&](gemmi::Miller h) { return asu.is_in(h) ? h : asu.to_asu(h, gops).first; };
+    std::map<gemmi::Miller, float> fmodel;
+    for (size_t m = 0; m < s.score_start.hkl.size(); m++)
+        fmodel[in_asu(s.score_start.hkl[m])] =
+                s.score_start.k_total[m] * std::abs(s.score_start.fcalc[m] + s.score_start.fmask[m]);
+    std::map<int, std::array<double, 2>> sums;   // class: sum |Fo - |Fmodel||, sum Fo
+    for (size_t i = 0; i < refl.size(); i++) {
+        const auto f = fmodel.find(in_asu({{refl[i].h, refl[i].k, refl[i].l}}));
+        if (cls[i] < 0 || f == fmodel.end() || refl[i].d < d_hi || refl[i].d > d_lo)
+            continue;
+        sums[cls[i]][0] += std::fabs(refl[i].F - f->second);
+        sums[cls[i]][1] += refl[i].F;
+    }
+    int worst = -1;
+    std::map<int, double> r;
+    for (const auto &[c, v] : sums) {
+        r[c] = v[1] > 0 ? v[0] / v[1] : 0;
+        if (worst < 0 || r[c] > r[worst])
+            worst = c;
+    }
+
+    // the classes whose count is nearest the model file's
+    const int n_test = model.meta.refinement.empty() ? -1 : model.meta.refinement[0].rfree_set_count;
+    std::vector<int> nearest;
+    if (n_test > 0) {
+        int best = -1;
+        for (const auto &[c, n] : counts) {
+            const int dist = std::abs(n - n_test);
+            if (best < 0 || dist < best) {
+                best = dist;
+                nearest.clear();
+            }
+            if (dist == best)
+                nearest.push_back(c);
+        }
+    }
+    auto is_nearest = [&](int c) { return std::find(nearest.begin(), nearest.end(), c) != nearest.end(); };
+
+    // the scatter of the classes but the worst: how much a class the model was refined against varies
+    double mean = 0, var = 0;
+    for (const auto &[c, v] : r)
+        if (c != worst)
+            mean += v / static_cast<double>(r.size() - 1);
+    for (const auto &[c, v] : r)
+        if (c != worst)
+            var += (v - mean) * (v - mean) / static_cast<double>(r.size() - 2);
+    const double scatter = std::sqrt(var);
+    const bool fitted_like_work = r.count(marked) > 0 && r[worst] - r[marked] > 2 * scatter;
+
+    int chosen = marked;
+    out.free_set_source = "STATUS";
+    if (worst != marked && !is_nearest(marked) && (is_nearest(worst) || fitted_like_work)) {
+        chosen = worst;
+        out.free_set_source = is_nearest(worst) ? "REMARK3_COUNT" : "WORST_FIT";
+    }
+    out.free_class = chosen;
+
+    std::string by_class, near;
+    for (const auto &[c, v] : r)
+        by_class += fmt::format("{}{} {:.3f}", by_class.empty() ? "" : ", ", c, v);
+    for (int c : nearest)
+        near += fmt::format("{}{}", near.empty() ? "" : ", ", c);
+    logger.Info("Refinement input: R as given by free-flag class: {}; the file marks class {}, the model fits class {} "
+                "worst ({:.1f} times the scatter {:.4f} of the others above the marked one){} -> class {} is the test "
+                "set ({})", by_class, marked, worst, scatter > 0 ? (r[worst] - r[marked]) / scatter : 0.0, scatter,
+                n_test > 0 ? fmt::format(", the model file's test-set count {} is nearest class {}", n_test, near)
+                           : std::string(", the model file gives no test-set count"),
+                chosen, out.free_set_source);
+    if (chosen != marked)
+        for (size_t i = 0; i < refl.size(); i++)
+            refl[i].rfree_flag = cls[i] == chosen;
+}
+
+}  // namespace
+
 RefineOnlyResult RunRefineOnly(RugnuxReflections data, const gemmi::Structure &model, double d_min_limit,
                                const RefinementSettings &settings, Logger &logger) {
     RefineOnlyResult out;
@@ -31,6 +152,26 @@ RefineOnlyResult RunRefineOnly(RugnuxReflections data, const gemmi::Structure &m
     out.space_group = sg->xhm();
     out.wavelength_A = data.wavelength_A > 0 ? data.wavelength_A : 1.0;
     std::vector<MergedReflection> &refl = data.reflections;
+    // free-flag classes, by reflection, and the reflections of each within the resolution range the model file
+    // records its refinement to, counted before any limit here
+    std::map<gemmi::Miller, int> class_of;
+    std::map<int, int> class_counts;
+    double d_hi = 0, d_lo = INFINITY;
+    if (!data.free_column.empty() && data.free_class.size() == refl.size()) {
+        if (!model.meta.refinement.empty()) {
+            const gemmi::RefinementInfo &ri = model.meta.refinement[0];
+            if (std::isfinite(ri.resolution_high) && ri.resolution_high > 0)
+                d_hi = ri.resolution_high - 1e-3;
+            if (std::isfinite(ri.resolution_low) && ri.resolution_low > 0)
+                d_lo = ri.resolution_low + 1e-3;
+        }
+        for (size_t i = 0; i < refl.size(); i++)
+            if (data.free_class[i] >= 0) {
+                class_of[{{refl[i].h, refl[i].k, refl[i].l}}] = data.free_class[i];
+                if (refl[i].d >= d_hi && refl[i].d <= d_lo)
+                    class_counts[data.free_class[i]]++;
+            }
+    }
     if (d_min_limit > 0)
         std::erase_if(refl, [&](const MergedReflection &r) { return r.d < d_min_limit; });
     if (!data.amplitudes) {
@@ -65,7 +206,17 @@ RefineOnlyResult RunRefineOnly(RugnuxReflections data, const gemmi::Structure &m
         if (!st.models.empty() && st.models[0].chains.size() != chains)
             logger.Info("The model's strict-NCS operators generate the rest of the asymmetric unit: {} chains -> {}",
                         chains, st.models[0].chains.size());
-        out.refinement = RefineModel(st, refl, data.cell, *sg, out.wavelength_A, logger, settings);
+        if (class_counts.size() > 2) {
+            ChooseFreeClass(refl, class_of, class_counts, d_hi, d_lo, st, data.cell, *sg, out.wavelength_A, settings,
+                            logger, out);
+            out.free_reflections = 0;
+            for (const MergedReflection &r : refl)
+                out.free_reflections += r.rfree_flag;
+        }
+        // the model has most likely been refined against every reflection of a free set drawn here
+        RefinementSettings refine = settings;
+        refine.as_given_competes = !out.free_flags_generated;
+        out.refinement = RefineModel(st, refl, data.cell, *sg, out.wavelength_A, logger, refine);
     } catch (const std::exception &e) {
         out.refinement = RefinementResult();
         out.refinement.failure_reason = std::string("refinement failed: ") + e.what();

rugnux/RugnuxRefine.h

@@ -20,6 +20,11 @@ struct RefineOnlyResult {
     double wavelength_A = 0;         // the file's, or 1 A where it states none
     bool french_wilson = false;      // amplitudes made here from the file's intensities
     bool free_flags_generated = false;   // the file has no free set: one was drawn here (AssignRfreeFlags)
+    // Where the file's free flags are classes (0-19), the class taken as the test set and what chose it (ChooseFreeClass
+    // in RugnuxRefine.cpp): "STATUS" the file's own mark, "REMARK3_COUNT" the class whose count the model file
+    // records and that the model fits worst, "WORST_FIT" the class the model fits worst. -1 / empty: not classes.
+    int free_class = -1;
+    std::string free_set_source;
     int reflections = 0;             // amplitudes refined against, after the resolution limit
     int free_reflections = 0;
     double d_min_limit = 0;          // the limit asked for (0: none)

rugnux/RugnuxScaleMerge.cpp

@@ -5424,6 +5424,9 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
             // refinement takes a law at a smaller lead over its null (AdoptTwinLaw).
             settings.twin_data_indicated = TwinningVerdictKey(result) == "INDICATED"
                                            && result.twinning.twin_fraction_source == TwinFractionSource::LTest;
+            // a model placed as given on a free set drawn here has seen those reflections (MODEL_FREE_SET): its R-free
+            // as given is no free statistic to keep it by
+            settings.as_given_competes = !(result.model_placement == "RIGID_BODY" && result.free_flags_source == "DRAWN");
             const PlacedModelOutput p = *placed_model_out;   // a copy: a rival taken below replaces *placed_model_out
             // In the group RefinableIn chose: the model's own where it has the rotations of the group written -
             // the screws are then the model's, as its fit already took them - so a model whose screws the merge

tests/RefineOnlyTest.cpp

@@ -221,6 +221,106 @@ TEST_CASE("RunRefineOnly scores a model against intensities, with amplitudes and
     CHECK(text.find("REFINE_FRENCH_WILSON= TRUE") != std::string::npos);
 }
 
+TEST_CASE("ParseRugnuxReflections keeps a structure-factor mmCIF's free-flag classes beside its status", "[refine_only]") {
+    // status marks class 5 free; the classes are h % 3 (three of them)
+    std::string cif = "data_r\n_cell.length_a 31\n_cell.length_b 37\n_cell.length_c 43\n_cell.angle_alpha 90\n"
+                      "_cell.angle_beta 90\n_cell.angle_gamma 90\n_symmetry.space_group_name_H-M 'P 21 21 21'\n"
+                      "loop_\n_refln.index_h\n_refln.index_k\n_refln.index_l\n_refln.status\n"
+                      "_refln.pdbx_r_free_flag\n_refln.F_meas_au\n_refln.F_meas_sigma_au\n";
+    for (int h = 1; h <= 6; h++)
+        for (int k = 1; k <= 3; k++)
+            cif += fmt::format("{} {} 1 {} {} {} 1\n", h, k, h % 3 == 2 ? 'f' : 'o', h % 3 == 2 ? 5 : h % 3, 10 * h);
+    const RugnuxReflections r = ParseRugnuxReflections(cif, "synthetic-sf.cif");
+    REQUIRE(r.reflections.size() == 18);
+    REQUIRE(r.free_class.size() == 18);
+    CHECK(r.n_free == 6);
+    for (size_t i = 0; i < r.reflections.size(); i++) {
+        CHECK(r.reflections[i].rfree_flag == (r.reflections[i].h % 3 == 2));
+        CHECK(r.free_class[i] == (r.reflections[i].h % 3 == 2 ? 5 : r.reflections[i].h % 3));
+    }
+}
+
+namespace {
+
+// The synthetic input with its reflections in 20 free-flag classes, class 5 marked free, and the intensities of
+// class `held_out` disturbed by 30 % at random: the class a model was not refined against, which it fits worst.
+// `count`: the reflections of that class.
+void ClassedRefineInput(gemmi::Structure &st, RugnuxReflections &data, int held_out, int &count) {
+    SyntheticRefineInput(st, data);
+    std::mt19937 rng(11);
+    std::uniform_real_distribution<float> noise(0.7f, 1.3f);
+    data.free_column = "FreeR_flag";
+    count = 0;
+    for (MergedReflection &m : data.reflections) {
+        const int c = (m.h * 7 + m.k * 13 + m.l * 17 + 1000) % 20;
+        data.free_class.push_back(c);
+        m.rfree_flag = c == 5;
+        data.n_free += m.rfree_flag;
+        if (c == held_out) {
+            m.I *= noise(rng) * noise(rng);
+            count++;
+        }
+    }
+}
+
+}  // namespace
+
+TEST_CASE("RunRefineOnly takes the free-flag class the model fits worst, by the model file's test-set count",
+          "[refine_only]") {
+    gemmi::Structure st;
+    RugnuxReflections data;
+    int count = 0;
+    ClassedRefineInput(st, data, 3, count);
+    st.meta.refinement.emplace_back();
+    st.meta.refinement[0].rfree_set_count = count;
+    RefinementSettings settings;
+    settings.score_only = true;
+    settings.gpu = false;
+    settings.nthreads = 2;
+    Logger logger("RefineOnlyTest");
+    const RefineOnlyResult r = RunRefineOnly(data, st, 0, settings, logger);
+    REQUIRE(r.refinement.ok);
+    CHECK(r.free_class == 3);
+    CHECK(r.free_set_source == "REMARK3_COUNT");
+    CHECK(r.free_reflections == count);
+    const std::string text = RenderReportText(BuildRefineReportDocument("out", "x-sf.cif", "model.cif", r), true);
+    CHECK(text.find("REFINE_FREE_SET_SOURCE= REMARK3_COUNT") != std::string::npos);
+    CHECK(text.find("REFINE_FREE_SET_CLASS= 3") != std::string::npos);
+}
+
+TEST_CASE("RunRefineOnly takes the worst-fit free-flag class where the model file gives no count", "[refine_only]") {
+    gemmi::Structure st;
+    RugnuxReflections data;
+    int count = 0;
+    ClassedRefineInput(st, data, 3, count);
+    RefinementSettings settings;
+    settings.score_only = true;
+    settings.gpu = false;
+    settings.nthreads = 2;
+    Logger logger("RefineOnlyTest");
+    const RefineOnlyResult r = RunRefineOnly(data, st, 0, settings, logger);
+    REQUIRE(r.refinement.ok);
+    CHECK(r.free_class == 3);
+    CHECK(r.free_set_source == "WORST_FIT");
+}
+
+TEST_CASE("RunRefineOnly keeps the marked free-flag class where the model fits it worst", "[refine_only]") {
+    gemmi::Structure st;
+    RugnuxReflections data;
+    int count = 0;
+    ClassedRefineInput(st, data, 5, count);
+    RefinementSettings settings;
+    settings.score_only = true;
+    settings.gpu = false;
+    settings.nthreads = 2;
+    Logger logger("RefineOnlyTest");
+    const RefineOnlyResult r = RunRefineOnly(data, st, 0, settings, logger);
+    REQUIRE(r.refinement.ok);
+    CHECK(r.free_class == 5);
+    CHECK(r.free_set_source == "STATUS");
+    CHECK(r.free_reflections == count);
+}
+
 TEST_CASE("RunRefineOnly scores a model with strict-NCS operators as the expanded model", "[refine_only]") {
     gemmi::Structure st;
     RugnuxReflections data;

tools/battery/model_check.py

@@ -454,12 +454,59 @@ def french_wilson(workdir, hkl, i, sigi, flags, cell, sg):
     return out.make_miller_array(), col("F"), col("SIGF"), col("FreeR_flag").astype(int)
 
 
+def free_class(hkl, f, cls, marked, xyz, cell, sg):
+    """The test set of deposited structure factors whose free flags are classes (0-19), as rugnux
+    --mode refine takes it (RugnuxRefine.cpp, ChooseFreeClass): the class the file marks (status 'f',
+    else flag 0) is not always the one the depositors held out. The deposited model fits the class it
+    was not refined against worst, and the test-set count its file records is that class's count,
+    both within the refinement's resolution range. The mark is kept where the model fits it worst or the
+    count is nearest it (STATUS); else the worst-fit class is taken where the count is nearest it
+    (REMARK3_COUNT) or where it is above the marked one by more than twice the scatter of the other
+    classes' R (WORST_FIT). The fit here: |Fcalc| without bulk solvent, an overall scale and B, on
+    the reflections finer than 6 A (the finer half on data coarser than that). Returns (class, source)."""
+    st = read_model(xyz)
+    ref = st.meta.refinement[0] if st.meta.refinement else None
+    d = cell.calculate_d_array(hkl)
+    d_hi = ref.resolution_high if ref is not None and ref.resolution_high > 0 else 0
+    d_lo = ref.resolution_low if ref is not None and ref.resolution_low > 0 else np.inf
+    in_range = (d >= d_hi - 1e-3) & (d <= d_lo + 1e-3)
+    classes = np.unique(cls[cls >= 0])
+    counts = {c: int(np.sum(in_range & (cls == c))) for c in classes}
+    n_test = ref.rfree_set_count if ref is not None else -1
+    nearest = []
+    if n_test > 0:
+        dist = {c: abs(n - n_test) for c, n in counts.items()}
+        nearest = [c for c in classes if dist[c] == min(dist.values())]
+    # the model as given on each class
+    ghkl, gf, _, gcls = to_group(hkl, f, np.ones(len(f)), cls, cell, sg, np.arange(len(f)))
+    fc_hkl, fc = fcalc(st, max(float(np.min(d)), 1.0))
+    _, i, j = np.intersect1d(hkl_key(ghkl), hkl_key(fc_hkl), return_indices=True)
+    fo, fcm, c_of = gf[i], fc[j], gcls[i]
+    dg = cell.calculate_d_array(ghkl[i])
+    s2 = 1.0 / dg ** 2
+    sel = (s2 > min(1 / 36.0, np.median(s2))) & (fo > 0) & (fcm > 0) & (c_of >= 0) & \
+        (dg >= d_hi - 1e-3) & (dg <= d_lo + 1e-3)
+    slope, lnk = np.polyfit(s2[sel], np.log(fo[sel] / fcm[sel]), 1)
+    resid = np.abs(fo - np.exp(lnk + slope * s2) * fcm)
+    r = {c: np.sum(resid[sel & (c_of == c)]) / np.sum(fo[sel & (c_of == c)]) for c in classes}
+    worst = max(r, key=r.get)
+    scatter = np.std([x for c, x in r.items() if c != worst], ddof=1)
+    if worst == marked or marked in nearest:
+        return marked, "STATUS"
+    if worst in nearest:
+        return worst, "REMARK3_COUNT"
+    if r.get(marked, np.inf) < r[worst] - 2 * scatter:
+        return worst, "WORST_FIT"
+    return marked, "STATUS"
+
+
 def depositor_data(pdb_id, cache_dir, workdir, cell_m, sg_m):
     """The first reflection block of the entry's -sf.cif as an MTZ (F, SIGF, FreeR_flag) in the
     model's cell and group, and what the amplitudes were made from. Amplitudes: F_meas, else the
     mean of F(+)/F(-); where only intensities were deposited (I, else the mean of I(+)/I(-)) they
-    are converted by French-Wilson. Free flags from _refln.status ('f' free), else from a numeric
-    _refln.pdbx_r_free_flag (0 when it takes more than two values - the CCP4 convention - else the
+    are converted by French-Wilson. Where a numeric _refln.pdbx_r_free_flag takes more than two
+    values, the free set is the class free_class decides; else free flags from _refln.status ('f'
+    free), else from a numeric _refln.pdbx_r_free_flag (0 when it takes more than two values - the CCP4 convention - else the
     minority value). Where neither marks a free set (e.g. a SHELXL entry whose free set was not
     deposited), a synthetic one is made by synthetic_free and noted. Returns None, reason when there
     is nothing usable, else the data and a note (None, or what was synthesised)."""
@@ -483,7 +530,19 @@ def depositor_data(pdb_id, cache_dir, workdir, cell_m, sg_m):
         if "status" in labels else np.array([])
     note = None
     v = b.make_int_array("pdbx_r_free_flag", -1) if "pdbx_r_free_flag" in labels else None
-    if len(status) == len(hkl) and np.any(status == "f"):
+    classes = np.unique(v[v >= 0]) if v is not None else np.array([])
+    if len(classes) > 2 and (len(status) != len(hkl) or np.any(status == "f")):
+        # free flags in classes: the test set by evidence (free_class)
+        marked = np.bincount(v[(status == "f") & (v >= 0)]).argmax() if len(status) == len(hkl) else 0
+        use = (v >= 0) & ((status == "f") | (status == "o") if len(status) == len(hkl) else True)
+        ok = use & np.isfinite(f) & (f > 0)
+        amp = f[ok] if amplitudes else np.sqrt(f[ok])
+        free, source = free_class(hkl[ok], amp, v[ok], marked,
+                                  deposition(pdb_id, cache_dir)["xyz"], cell_m, sg_m)
+        flags = np.where(v == free, 0, 1)
+        if free != marked:
+            note = f"free set: flag class {free} ({source}), not the class {marked} the file marks"
+    elif len(status) == len(hkl) and np.any(status == "f"):
         use = (status == "f") | (status == "o")
         flags = np.where(status == "f", 0, 1)
     elif v is not None and len(np.unique(v[v >= 0])) >= 2:
@@ -493,14 +552,17 @@ def depositor_data(pdb_id, cache_dir, workdir, cell_m, sg_m):
         flags = np.where(v == free, 0, 1)
     else:
         use = np.ones(len(hkl), bool) if len(status) != len(hkl) else (status == "o")
-        flags = np.ones(len(hkl), int)
+        flags = None
         note = f"deposited structure factors carry no free set: a synthetic {SYNTHETIC_FREE:.0%} one is used"
     use &= np.isfinite(f) & np.isfinite(s) & (s > 0)
     if not np.allclose(b.cell.parameters, cell_m.parameters, rtol=0.02, atol=0.5):
         return None, "deposited structure factors in another cell"
+    synthetic = flags is None
+    if synthetic:
+        flags = np.ones(len(hkl), int)
     hkl, f, s, flags = to_group(hkl[use], f[use], s[use], flags[use], cell_m, sg_m,
                                 np.arange(use.sum()))
-    if note:
+    if synthetic:
         flags = synthetic_free(hkl)
     if not amplitudes:
         hkl, f, s, flags = french_wilson(workdir, hkl, f, s, flags, cell_m, sg_m)

writer/RugnuxReflectionFile.cpp

@@ -107,6 +107,8 @@ std::vector<std::string> SfCifSpec(const gemmi::cif::Loop &loop) {
             continue;
         spec.push_back(l);
     }
+    // the numeric classes as well, whichever marks the free set (RugnuxReflections::free_class)
+    spec.emplace_back("pdbx_r_free_flag FreeR_class I 0");
     return spec;
 }
 
@@ -143,6 +145,7 @@ RugnuxReflections RugnuxReflectionsFromMtz(const gemmi::Mtz &mtz, const std::str
     }
 
     const gemmi::Mtz::Column *free = FindFree(mtz);
+    const gemmi::Mtz::Column *free_class = mtz.column_with_label("FreeR_class", nullptr, 'I');
     std::vector<int> raw_free;
     const size_t stride = mtz.columns.size();
     for (size_t r = 0; r < static_cast<size_t>(mtz.nreflections); r++) {
@@ -184,6 +187,10 @@ RugnuxReflections RugnuxReflectionsFromMtz(const gemmi::Mtz &mtz, const std::str
             const float f = row[free->idx];
             raw_free.push_back(std::isfinite(f) ? static_cast<int>(std::lround(f)) : -1);
         }
+        if (free_class != nullptr) {
+            const float f = row[free_class->idx];
+            out.free_class.push_back(std::isfinite(f) ? static_cast<int>(std::lround(f)) : -1);
+        }
     }
     if (out.reflections.empty())
         throw std::runtime_error(name + " holds no reflection with a value in " + out.value_column);

writer/RugnuxReflectionFile.h

@@ -29,6 +29,10 @@ struct RugnuxReflections {
     std::string value_column, sigma_column;
     std::string free_column;         // empty: the file has no free-set flags
     int n_free = 0;
+    // A structure-factor mmCIF's numeric free-set classes (_refln.pdbx_r_free_flag) per reflection, -1 where it
+    // gives none; empty for an MTZ or where the file has no such column. free_column marks one class of them
+    // (_refln.status 'f', or the flag by the convention above), which need not be the depositors' test set.
+    std::vector<int> free_class;
 };
 
 // The columns of `mtz` (an MTZ, or an mmCIF converted to one), found as RugnuxReflections describes: