From d9d939ae44e5f9eacee7170c45cb0f70fa5e51ee Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 16:06:50 +0100 Subject: [PATCH 1/5] RESULT-CHANGING: --mode refine: the test set of a structure-factor mmCIF whose free flags are classes, decided by the model as given A deposited structure-factor mmCIF marks the free set with _refln.status 'f', and the reader took that mark over _refln.pdbx_r_free_flag. Where the flags are classes (0-19), the status can mark a class the depositors' refinement did not hold out: the deposited model then fits the marked class like working data, R-work comes out above R-free, and refinement seems to make the model worse. The reader now keeps the numeric classes beside the status (RugnuxReflections::free_class, an extra FreeR_class column in the conversion). Where there are more than two classes, RunRefineOnly scores the model as given on each class (score-only, untwinned, no solvent build), within the resolution range the model file records its refinement to, and counts each class's reflections in that range: - STATUS: the marked class is the one the model fits worst, or its count is the nearest to the model file's test-set count (REMARK 3 FREE R VALUE TEST SET COUNT / _refine.ls_number_reflns_R_free); - REMARK3_COUNT: else the worst-fit class, where its count is the nearest; - WORST_FIT: else the worst-fit class, where it is above the marked one by more than twice the scatter of the other classes' R (the marked class is fitted like a working one); otherwise the mark stands. The count is corroboration only: over the battery's PDB cache (106 files with 20 classes) the marked class is within a few reflections of the count on nearly all, and another class has the count exactly on about 7 of them (20 classes of about the same size), so the count alone would move correct marks. The choice is logged with R per class and reported as REFINE_FREE_SET_SOURCE / REFINE_FREE_SET_CLASS. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/RUGNUX_MODEL_AND_REFERENCE.md | 8 +- docs/RUGNUX_REPORT.md | 5 +- rugnux/ResultReport.cpp | 4 + rugnux/RugnuxRefine.cpp | 148 +++++++++++++++++++++++++++++ rugnux/RugnuxRefine.h | 5 + tests/RefineOnlyTest.cpp | 100 +++++++++++++++++++ writer/RugnuxReflectionFile.cpp | 7 ++ writer/RugnuxReflectionFile.h | 4 + 8 files changed, 279 insertions(+), 2 deletions(-) diff --git a/docs/RUGNUX_MODEL_AND_REFERENCE.md b/docs/RUGNUX_MODEL_AND_REFERENCE.md index 7aca9af63..8c3da7c74 100644 --- a/docs/RUGNUX_MODEL_AND_REFERENCE.md +++ b/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 - `_refined_model.cif`/`.pdb`, `_refined_maps.mtz`, `_blobs.json` - with the `REFINE_` keys in `_report.txt`, beside `REFINE_FRENCH_WILSON=` and `REFINE_FREE_FLAGS= FILE|GENERATED` saying what had to be made of the input. diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index e3b9ceb6f..4351c3243 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -1300,4 +1300,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). diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 1a85a59ba..73171fbb9 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -3393,6 +3393,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); diff --git a/rugnux/RugnuxRefine.cpp b/rugnux/RugnuxRefine.cpp index fd8e51d90..eb75ba126 100644 --- a/rugnux/RugnuxRefine.cpp +++ b/rugnux/RugnuxRefine.cpp @@ -3,13 +3,134 @@ #include "RugnuxRefine.h" +#include +#include #include +#include + +#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 &refl, const std::map &class_of, + const std::map &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 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 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 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> 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 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 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(r.size() - 1); + for (const auto &[c, v] : r) + if (c != worst) + var += (v - mean) * (v - mean) / static_cast(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 &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 class_of; + std::map 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,6 +206,13 @@ 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()); + 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; + } out.refinement = RefineModel(st, refl, data.cell, *sg, out.wavelength_A, logger, settings); } catch (const std::exception &e) { out.refinement = RefinementResult(); diff --git a/rugnux/RugnuxRefine.h b/rugnux/RugnuxRefine.h index 78e2bf2af..64336d304 100644 --- a/rugnux/RugnuxRefine.h +++ b/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) diff --git a/tests/RefineOnlyTest.cpp b/tests/RefineOnlyTest.cpp index 37a62d6a9..1377319aa 100644 --- a/tests/RefineOnlyTest.cpp +++ b/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 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; diff --git a/writer/RugnuxReflectionFile.cpp b/writer/RugnuxReflectionFile.cpp index 5fd297502..a6f8a1aaf 100644 --- a/writer/RugnuxReflectionFile.cpp +++ b/writer/RugnuxReflectionFile.cpp @@ -107,6 +107,8 @@ std::vector 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 raw_free; const size_t stride = mtz.columns.size(); for (size_t r = 0; r < static_cast(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(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(std::lround(f)) : -1); + } } if (out.reflections.empty()) throw std::runtime_error(name + " holds no reflection with a value in " + out.value_column); diff --git a/writer/RugnuxReflectionFile.h b/writer/RugnuxReflectionFile.h index 08a2de429..7d080db88 100644 --- a/writer/RugnuxReflectionFile.h +++ b/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 free_class; }; // The columns of `mtz` (an MTZ, or an mmCIF converted to one), found as RugnuxReflections describes: From 4dce68e0076833ea90b88c4917c468b9e94359e6 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 16:06:58 +0100 Subject: [PATCH 2/5] model_check: the depositors' test set of a structure-factor file whose free flags are classes, by the rule --mode refine uses depositor_data took _refln.status 'f' first, so where the flags are classes (0-19) and the status marks a class the depositors did not hold out, REFMAC's depositor-flag and depositor-data R-free were taken on working reflections. free_class applies the rule of RunRefineOnly (ChooseFreeClass): the mark stands where the deposited model fits it worst or its count is nearest the entry's test-set count, else the worst-fit class is taken where the count is nearest it or it stands above the mark by twice the classes' scatter. The fit here is gemmi's Fcalc without bulk solvent, an overall scale and B, on the reflections finer than 6 A within the refinement's range. Over the 119 multi-class structure-factor files in the battery's PDB cache this changes the free class of one entry (5nw5: the class whose measured reflections match the entry's test-set count exactly, where the status marks another); the depdata/depflags REFMAC numbers of the others are unchanged. Kept as its own commit: it moves battery REFMAC numbers. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- tools/battery/model_check.py | 72 +++++++++++++++++++++++++++++++++--- 1 file changed, 67 insertions(+), 5 deletions(-) diff --git a/tools/battery/model_check.py b/tools/battery/model_check.py index 686ae4eb2..58efc4444 100644 --- a/tools/battery/model_check.py +++ b/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) From 022ebf8d773c88fba1e1a89dde624a9a640c0d82 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 16:07:07 +0100 Subject: [PATCH 3/5] MolecularReplacement: "not refined" and "likelihood Z not needed" in the log for the NaN sentinels, and a line where no further copy is kept A first-copy placement that clashes too much is never refined, and its R-factors were printed as "R-work nan R-free nan"; a further copy whose translation-function Z already stands out has its likelihood Z skipped, printed as "refined nan Z nan". Neither is a failure, so the log now says so. Where none of a copy search's trials stands out, packs and lowers R-free, the search stopped without a word; it now says how many trials there were and how many copies were kept. Log text only. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- .../MolecularReplacement.cpp | 23 ++++++++++++++----- 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp index dc18038e8..59b52ed94 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.cpp +++ b/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)); From f1480ac7dc401af42f9ce2f53068176a5a377030 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 16:07:32 +0100 Subject: [PATCH 4/5] RESULT-CHANGING: Refinement: the model as given is a candidate of the final choice, where its R-free is a free statistic The choice was made among the trials alone: the R-free of the model as given was recorded but never competed, so a refinement whose every trial raised R-free still replaced the model (a related model that cannot be refined to a good R came out above its starting R-free). The model as given, scaled to the data (the scale and bulk-solvent fit of the start), is now the first candidate in order of trust, judged by the rule every trial is held to: the refinement is kept where its geometry is within bounds and its R-free lower beyond the tie, else the model is given back as it came. The solvent build edits the model in place (waters put in, sites relabelled), so its atoms and residues as given are kept and put back, with the placeholder pool and the relabels, and the solvent build is reported as not kept. REFINE_AS_GIVEN says which; the twin decision then compares the arms as they ended, as before. Only where the free set is free to the model as given (RefinementSettings::as_given_competes): a model placed as given on a free set drawn by the run (the MODEL_FREE_SET case; in --mode refine, a file without free flags) has most likely been refined against those reflections, its R-free as given is close to a work R, and keeping it by that number would hand back the deposited model. On a first subset battery without this condition the model as given was kept on four open-arm sets with R-work then above R-free. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/CPU_DATA_ANALYSIS_REFINEMENT.md | 6 +++ docs/RUGNUX_REPORT.md | 1 + .../structure_refinement/Refinement.cpp | 48 +++++++++++++++++++ .../structure_refinement/Refinement.h | 5 ++ rugnux/ResultReport.cpp | 1 + rugnux/RugnuxRefine.cpp | 5 +- rugnux/RugnuxScaleMerge.cpp | 3 ++ 7 files changed, 68 insertions(+), 1 deletion(-) diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index df3a89029..a5cc1398f 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -322,6 +322,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 diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 4351c3243..f8e97e218 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -1229,6 +1229,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 `_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 | diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index 11ea25145..c61cf77d0 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -2182,6 +2182,27 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector given_atoms; + std::vector> 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 relabels; @@ -2358,6 +2379,11 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vectorevaluations; ctx.evaluations -= ahead_evaluations; diff --git a/image_analysis/structure_refinement/Refinement.h b/image_analysis/structure_refinement/Refinement.h index cb554f74f..8f0da664f 100644 --- a/image_analysis/structure_refinement/Refinement.h +++ b/image_analysis/structure_refinement/Refinement.h @@ -75,6 +75,10 @@ struct RefinementSettings { // A refinement against the data as twinned by each twin law the lattice admits (RefinementTwinLaws), beside the // untwinned one, and the model of whichever the evidence favours (RefineModel). bool twin = true; + // 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 @@ -156,6 +160,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; diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 73171fbb9..09df06751 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -547,6 +547,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}")); diff --git a/rugnux/RugnuxRefine.cpp b/rugnux/RugnuxRefine.cpp index eb75ba126..fe217da2b 100644 --- a/rugnux/RugnuxRefine.cpp +++ b/rugnux/RugnuxRefine.cpp @@ -213,7 +213,10 @@ RefineOnlyResult RunRefineOnly(RugnuxReflections data, const gemmi::Structure &m for (const MergedReflection &r : refl) out.free_reflections += r.rfree_flag; } - out.refinement = RefineModel(st, refl, data.cell, *sg, out.wavelength_A, logger, settings); + // 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(); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 0d7d38991..3d892d08b 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -5380,6 +5380,9 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { settings.gpu_resident_bytes = rsm && rsm->OnGPU() ? rsm->DeviceBytes() : 0; // A cut (the broker's deadline or cancel, Ctrl+C) cuts a running refinement short, keeping what it has settings.cancel = cancelled_.RunFlag(); + // 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 From 27e386bfd2a469ac2f327bf8e6cc760ea785a3da Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 18:12:23 +0100 Subject: [PATCH 5/5] RESULT-CHANGING: Refinement: the U at its balanced weight also where the U at a quarter lowers -ln L free The displacement model is chosen on R-free (beyond the tie) and -ln L of the free reflections together, but the fewer-trials rule ran the U at its balanced weight only where the U at a quarter of it lowered R-free. Where the quarter fitted the free reflections better than where the weights ended without lowering R-free, the balanced weight was left out, and with it the trial the choice would have kept: on a twinned refine-only case the quarter had the best -ln L free of the displacement trials and an R-free a little above the isotropic one, so isotropic was kept, while the balanced weight - run before the rule - lowers both and is kept. The condition is now the decision's own: R-free or -ln L free lower than where the weights ended. The -ln L free there is logged. Subset battery (5 refine-arm, 9 open, 2 inhouse sets) against rc181 e56b19abb: refined R-free identical on all 16; the balanced-weight U trial ran on one more set (both twin arms, no choice changed); target evaluations 49023 -> 49339 (+0.6%). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/CPU_DATA_ANALYSIS_REFINEMENT.md | 6 ++++-- .../structure_refinement/Refinement.cpp | 18 +++++++++++++----- 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index a5cc1398f..124e2d628 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/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). diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index c61cf77d0..97f44af58 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -2870,8 +2870,8 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector 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 @@ -2955,13 +2957,19 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector= 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 needed(trials.size(), 0);