Refinement: placement steps kept only where R-free drops; a wider weight ladder

On deposited models the per-residue B shift raised R-free at every resolution tried (2.56 A: 0.320 ->
0.327) and the rigid bodies rarely helped; each is now kept only where it lowers R-free by more than the
0.001 tie, otherwise the model goes on as it was. The positions' weight ladder runs from 1/1024 to 1/32 of
the gradient-balanced weight: the R-free optimum sat near 1/500 at 2.6 A and near 1/50 at 1.1-1.7 A.

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

image_analysis/structure_refinement/Refinement.cpp

@@ -411,7 +411,7 @@ private:
 };
 
 // The placement stage: rigid bodies, then group B.
-void RefinePlacement(Context &ctx, RefinementResult &result) {
+void RefinePlacement(Context &ctx, bool group_b, RefinementResult &result) {
     std::vector<std::vector<int>> bodies, groups;
     size_t i = 0;
     for (const gemmi::Chain &ch : ctx.Model().chains) {
@@ -429,8 +429,11 @@ void RefinePlacement(Context &ctx, RefinementResult &result) {
         if (body.size() >= 50)
             bodies.push_back(body);
     }
-    result.rigid_bodies = static_cast<int>(bodies.size());
+    // Each step is kept only where it lowers R-free beyond the tie: a deposited model was refined atom by atom
+    // already, and moving its chains or shifting its B factors residue by residue then only makes it worse.
+    RFactors before = ctx.Refit();
     if (!bodies.empty()) {
+        const std::vector<gemmi::Position> pos0 = ctx.Positions();
         std::vector<double> q(6 * bodies.size(), 0.0);
         auto target = std::make_unique<RigidBodiesTarget>(ctx, bodies);
         RigidBodiesTarget *rb = target.get();
@@ -438,9 +441,17 @@ void RefinePlacement(Context &ctx, RefinementResult &result) {
         ceres::GradientProblemSolver::Summary summary;
         ceres::Solve(SolverOptions(30), problem, q.data(), &summary);
         rb->Place(q.data());
+        const RFactors after = ctx.Refit();
+        if (after.r_free < before.r_free - R_FREE_TIE) {
+            result.rigid_bodies = static_cast<int>(bodies.size());
+            before = after;
+        } else {
+            ctx.SetPositions(pos0);
+            before = ctx.Refit();
+        }
     }
-    ctx.Refit();
-    {
+    if (group_b) {
+        const std::vector<double> b0 = ctx.B();
         std::vector<double> q(groups.size(), 0.0);
         auto target = std::make_unique<GroupBTarget>(ctx, groups);
         GroupBTarget *gb = target.get();
@@ -448,6 +459,9 @@ void RefinePlacement(Context &ctx, RefinementResult &result) {
         ceres::GradientProblemSolver::Summary summary;
         ceres::Solve(SolverOptions(30), problem, q.data(), &summary);
         gb->Place(q.data());
+        result.group_b = ctx.Refit().r_free < before.r_free - R_FREE_TIE;
+        if (!result.group_b)
+            ctx.SetB(b0);
     }
 }
 
@@ -642,12 +656,12 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
                 ctx.atoms.size(), ctx.refl.Size(), result.d_min, r.n_free, result.on_gpu ? "GPU" : "CPU", r.r_work,
                 r.r_free);
 
-    RefinePlacement(ctx, result);
+    RefinePlacement(ctx, settings.group_b, result);
     r = ctx.Refit();
     result.r_work_placement = r.r_work;
     result.r_free_placement = r.r_free;
-    logger.Info("Refinement, placement ({} rigid bodies, B per residue): R-work {:.4f}, R-free {:.4f}",
-                result.rigid_bodies, r.r_work, r.r_free);
+    logger.Info("Refinement, placement ({} rigid bodies kept, B per residue {}): R-work {:.4f}, R-free {:.4f}",
+                result.rigid_bodies, result.group_b ? "kept" : "not kept", r.r_work, r.r_free);
 
     if (settings.restrained) {
         RefinementTopology topo = BuildRefinementTopology(ctx.st, settings.riding_hydrogens);

image_analysis/structure_refinement/Refinement.h

@@ -22,8 +22,8 @@ class Logger;
 // Refinement of an atomic model already placed against the data (ValidateAgainstModel's placed model, or a
 // molecular-replacement placement), in two stages:
 //
-// 1. Placement: each chain as a rigid body, then one B shift per residue, then the bulk solvent and the
-//    scale. Needs no dictionary.
+// 1. Placement: each chain as a rigid body, then one B shift per residue, each kept only where it lowers
+//    R-free. Needs no dictionary.
 // 2. Restrained: every atom's position and isotropic B, against the maximum-likelihood amplitude target with
 //    the geometry the dictionary describes (MonomerLibrary.h: bonds, angles, planes, chiral volumes,
 //    torsions, repulsions; B similarity of bonded atoms). Residues it does not describe are held fixed and
@@ -35,11 +35,12 @@ class Logger;
 // statistics and the map coefficients go out. Nothing it finds feeds back into processing decisions.
 struct RefinementSettings {
     bool restrained = true;          // stage 2; false stops after the placement stage
+    bool group_b = true;             // the placement stage's B shift per residue
     int macro_cycles = 5;
     int xyz_iterations = 25;         // L-BFGS iterations per macro-cycle
     int b_iterations = 15;
     // The positions' X-ray weight, times the gradient-balanced one (BalancedWeight), tried in this order.
-    std::vector<double> weight_ladder = {0.0078125, 0.015625, 0.03125, 0.0625};
+    std::vector<double> weight_ladder = {1.0 / 1024, 1.0 / 512, 1.0 / 256, 1.0 / 128, 1.0 / 64, 1.0 / 32};
     // The X-ray weights of the B refinement tried next, times the positions' chosen weight.
     std::vector<double> adp_ladder = {1.0, 4.0, 16.0};
     bool riding_hydrogens = true;    // the dictionary's hydrogens added, riding on their parents
@@ -68,7 +69,8 @@ struct RefinementResult {
     // The model as given, scaled to these reflections (hydrogens removed, B made isotropic).
     double r_work_start = 1, r_free_start = 1;
     // After the placement stage.
-    int rigid_bodies = 0;
+    int rigid_bodies = 0;            // the chains moved as rigid bodies, where that was kept
+    bool group_b = false;            // the B shift per residue was kept
     double r_work_placement = 1, r_free_placement = 1;
     // After the restrained stage (NaN where it did not run).
     bool restrained = false;

tools/rugnux_refine_test.cpp

@@ -28,6 +28,9 @@ void usage() {
                  "  -N <threads>        CPU threads (default 8)\n"
                  "  --no-gpu            the CPU engine\n"
                  "  --placement-only    rigid bodies and group B only\n"
+                 "  --no-group-b        no B shift per residue in the placement stage\n"
+                 "  --no-hydrogens      no riding hydrogens\n"
+                 "  --adp-ladder <a,..> the B refinement's weights, times the positions' chosen weight\n"
                  "  --cycles <n>        macro-cycles per weight (default as in RefinementSettings)\n"
                  "  --xyz-iterations <n>, --b-iterations <n>  L-BFGS iterations per macro-cycle\n"
                  "  --ladder <a,b,...>  the weight ladder, times the gradient-balanced weight\n"
@@ -57,7 +60,11 @@ int main(int argc, char **argv) {
         else if (a == "--no-gpu") {
             settings.gpu = false;
             hide_gpus();
-        } else if (a == "--placement-only")
+        } else if (a == "--no-group-b")
+            settings.group_b = false;
+        else if (a == "--no-hydrogens")
+            settings.riding_hydrogens = false;
+        else if (a == "--placement-only")
             settings.restrained = false;
         else if (a == "--cycles" && i + 1 < argc)
             settings.macro_cycles = std::stoi(argv[++i]);