Refinement: B restraints in proportion to B (Kullback-Leibler form) with contacts, an own ADP weight, overall B into the atoms

- B similarity: E = (B_i - B_j)^2 / (kappa^2 B_i B_j), the symmetric Kullback-Leibler divergence of two isotropic Gaussians (Murshudov et al. 2011), so sigma = kappa sqrt(B_i B_j) instead of a fixed 1.5-4.5 A^2. kappa 0.1 / 0.133 / 0.2 / 0.3 for main-chain bonded / angle-related and side-chain pairs (REFMAC's logs list sigma ~0.1 B for main-chain bonded pairs at B 11-140; the ratios are those of the old sigmas, which these equal at B = 15). New class: atoms within 4 A that no bond or angle relates, waters and symmetry mates included, kappa 0.5 (BuildBContacts, sharing the neighbour search with BuildRepulsions). CPU and GPU. - The B factors' X-ray weight is balanced on its own (BalancedADPWeight: B-restraint over X-ray B-gradient on the model with B shaken by 20%) instead of being a multiple of the positional weight; the ladder tries B held and refined at that weight, the trajectory's B steps divide and multiply it by four, positional steps and the final ladder leave it. - Before each B-refining macro-cycle the overall isotropic B of the scale moves into the atoms and the solvent's B (F_model unchanged), so the two do not drift against each other. Tried and not kept: xyz and B as one L-BFGS vector (the B block scaled by its curvature, the positions by the median X-ray curvature): no better than alternation with a memory per block on 12 starts (mean R-free 0.265 vs 0.262), worse on 7arr/7d98; a refit before each B step: no gain for +20-40% time. Refinement battery MD starts, median R-free by the REFMAC 0-cycle scorer, previous commit -> this (REFMAC to convergence): 5f6m x8 0.167 -> 0.167 (0.168), 8dqb x8 0.291 -> 0.294 (0.284), 7arr x4 0.270 -> 0.278 (0.311), 11if x4 0.245 -> 0.244 (0.251), 5reo x4 0.312 -> 0.258 (0.261: the ladder no longer keeps B held where B is the model's largest error), 7d98 x4 0.342 -> 0.335 (0.361). Deposited models: unchanged (B stays held there). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit ed292f35b482b696d81f71f1d97e0016ca05b144)
2026-10-26 01:38:24 +01:00
1 parent 6369a263b1
commit bda90dbdd5
@@ -551,10 +551,12 @@ library for structural biology" (2022), J. Open Source Softw. 7, 4200
 A. Blondel and M. Karplus, "New formulation for derivatives of torsion angles and improper torsion angles
 in molecular mechanics: elimination of singularities" (1996), J. Comput. Chem. 17, 1132-1141
 [doi:10.1002/(SICI)1096-987X(19960715)17:9<1132::AID-JCC5>3.0.CO;2-T](https://doi.org/10.1002/(SICI)1096-987X(19960715)17:9%3C1132::AID-JCC5%3E3.0.CO;2-T).
-The B factors of bonded atoms are restrained to each other, as REFMAC5 restrains them (Murshudov et al.
-1997, above). The weight of the X-ray term against the geometry starts from the ratio of the two
-gradients and is chosen among several full refinements by R-free with the geometry held within bounds,
-as phenix.refine chooses it: P. V. Afonine, R. W. Grosse-Kunstleve, N. Echols, J. J. Headd,
+The B factors of bonded atoms, and of atoms in contact, are restrained to each other in the
+Kullback-Leibler form REFMAC5 uses, whose tolerance grows with B (Murshudov et al. 2011, above), and the
+overall isotropic B of the scale is moved into the atoms as phenix.refine moves it (Afonine et al. 2012,
+below). The weights of the X-ray term against the geometry and against the B restraints start from the
+ratios of the gradients and are chosen among short refinements by R-free with the geometry held within
+bounds, as phenix.refine chooses them: P. V. Afonine, R. W. Grosse-Kunstleve, N. Echols, J. J. Headd,
 N. W. Moriarty, M. Mustyakimov, T. C. Terwilliger, A. Urzhumtsev, P. H. Zwart and P. D. Adams,
 "Towards automated crystallographic structure refinement with phenix.refine" (2012), Acta Cryst. D68,
 352-367 [doi:10.1107/S0907444912001308](https://doi.org/10.1107/S0907444912001308). R-free itself:
@@ -266,6 +266,34 @@ struct Context {
        return ComputeRFactors(refl, fmodel);
    }

+    // The overall isotropic B of the last Refit() moved into the atoms, and into the solvent's B so that F_model stays
+    // as it is (but where an atom meets a bound of B): left in the scale, it and the atoms' B drift against each
+    // other as the B factors are refined - and the B restraints, in proportion to B, depend on which holds it.
+    // As phenix.refine moves the overall B into the atoms (Afonine et al. (2012) Acta Cryst. D68, 352-367)
+    void MoveOverallB() {
+        const gemmi::Mat33 b = cell.orth.mat.multiply(b_star.as_mat33()).multiply(cell.orth.mat.transpose());
+        const double b_iso = (b[0][0] + b[1][1] + b[2][2]) / 3.0;
+        double b_min = b_ceiling;
+        for (size_t i = 0; i < atoms.size(); i++)
+            if (rider_parent.empty() || rider_parent[i] < 0) {
+                atoms[i]->b_iso = static_cast<float>(std::clamp(atoms[i]->b_iso + b_iso, B_FLOOR, b_ceiling));
+                b_min = std::min<double>(b_min, atoms[i]->b_iso);
+            }
+        MoveRiders();
+        b_floor = std::max(B_FLOOR, b_min - B_STEP_DOWN);
+        engine->SetBlur(RefinementBlur(cell, *sg, d_min, b_floor));
+        // b_iso I in Cartesian is b_iso F F^T in the reciprocal basis
+        const gemmi::Mat33 f = cell.frac.mat.multiply(cell.frac.mat.transpose());
+        b_star = {b_star.u11 - b_iso * f[0][0], b_star.u22 - b_iso * f[1][1], b_star.u33 - b_iso * f[2][2],
+                  b_star.u12 - b_iso * f[0][1], b_star.u13 - b_iso * f[0][2], b_star.u23 - b_iso * f[1][2]};
+        b_sol += b_iso;
+        for (size_t m = 0; m < refl.Size(); m++) {
+            rows.k_total[m] *= static_cast<float>(std::exp(b_iso * refl.inv_d2[m] / 4.0));
+            rows.fmask[m] *= static_cast<float>(std::exp(-b_iso * refl.inv_d2[m] / 4.0));
+        }
+        engine->SetRowModel(rows);
+    }
+
    // Every term of the score of the atoms as they are, at the scale, the solvent and the likelihood parameters
    // of the last Refit().
    RefinementScore Score() {
@@ -448,8 +476,22 @@ public:
    AtomTarget(Context &ctx, Mode mode, const std::vector<int> &param_atoms, double weight, int restraint_set)
        : ctx_(ctx), mode_(mode), param_atoms_(param_atoms), weight_(weight), restraint_set_(restraint_set) {}

-    bool Evaluate(const double *p, double *cost, double *gradient) const {
-        const RefinementStage stage(mode_ == Mode::Xyz ? "eval.xyz" : "eval.b");
+    // The parameters of the atoms as they are: x y z per atom, or B per atom.
+    std::vector<double> Get() const {
+        std::vector<double> p;
+        for (int i : param_atoms_) {
+            const gemmi::Atom &a = *ctx_.atoms[i];
+            if (mode_ == Mode::Xyz) {
+                p.push_back(a.pos.x);
+                p.push_back(a.pos.y);
+                p.push_back(a.pos.z);
+            } else
+                p.push_back(a.b_iso);
+        }
+        return p;
+    }
+    // The atoms set to the parameters, their riders with them.
+    void Set(const double *p) const {
        for (size_t k = 0; k < param_atoms_.size(); k++) {
            gemmi::Atom &a = *ctx_.atoms[param_atoms_[k]];
            if (mode_ == Mode::Xyz)
@@ -458,32 +500,36 @@ public:
                a.b_iso = static_cast<float>(std::clamp(p[k], ctx_.b_floor, ctx_.b_ceiling));
        }
        ctx_.MoveRiders();
+    }
+    // the most a parameter moves in the first trial of a line search
+    double MaxStep() const { return mode_ == Mode::Xyz ? MAX_STEP_XYZ : MAX_STEP_B; }
+
+    bool Evaluate(const double *p, double *cost, double *gradient) const {
+        const RefinementStage stage(mode_ == Mode::Xyz ? "eval.xyz" : "eval.b");
+        Set(p);
        ctx_.evaluations++;
        std::vector<std::array<double, 4>> gx, gr;
        double xray = 0, restraints = 0;
        {
            const RefinementStage s1("eval.engine");
-            if (restraint_set_ >= 0)
-                xray = ctx_.engine->EvaluateWithRestraints(ctx_.Model(), restraint_set_, gradient != nullptr ? &gx : nullptr,
-                                                           gradient != nullptr ? &gr : nullptr, restraints);
-            else
-                xray = ctx_.engine->Evaluate(ctx_.Model(), gradient != nullptr ? &gx : nullptr);
+            xray = ctx_.engine->EvaluateWithRestraints(ctx_.Model(), restraint_set_, gradient != nullptr ? &gx : nullptr,
+                                                       gradient != nullptr ? &gr : nullptr, restraints);
        }
-        if (gradient != nullptr)
-            ctx_.FoldRiders(gx);
        *cost = weight_ * xray + restraints;
-        if (gradient != nullptr)
-            for (size_t k = 0; k < param_atoms_.size(); k++) {
-                const int i = param_atoms_[k];
-                if (mode_ == Mode::Xyz)
-                    for (int c = 0; c < 3; c++)
-                        gradient[3 * k + c] = weight_ * gx[i][c] + (restraint_set_ >= 0 ? gr[i][c] : 0.0);
-                else {
-                    // held at a bound, B does not move further out
-                    const bool held = (p[k] <= ctx_.b_floor && gx[i][3] > 0) || (p[k] >= ctx_.b_ceiling && gx[i][3] < 0);
-                    gradient[k] = held ? 0.0 : weight_ * gx[i][3] + (restraint_set_ >= 0 ? gr[i][3] : 0.0);
-                }
+        if (gradient == nullptr)
+            return std::isfinite(*cost);
+        ctx_.FoldRiders(gx);
+        for (size_t k = 0; k < param_atoms_.size(); k++) {
+            const int i = param_atoms_[k];
+            if (mode_ == Mode::Xyz)
+                for (int c = 0; c < 3; c++)
+                    gradient[3 * k + c] = weight_ * gx[i][c] + gr[i][c];
+            else {
+                // held at a bound, B does not move further out
+                const bool held = (p[k] <= ctx_.b_floor && gx[i][3] > 0) || (p[k] >= ctx_.b_ceiling && gx[i][3] < 0);
+                gradient[k] = held ? 0.0 : weight_ * gx[i][3] + gr[i][3];
            }
+        }
        return std::isfinite(*cost);
    }

@@ -495,10 +541,10 @@ private:
    int restraint_set_;
 };

-// The minimiser's diagonal for the B factors of `param_atoms` (RefinementMinimiser::SetDiagonal): weight times
-// the X-ray target's curvature in B plus the B restraints', each element at least a thousandth of the median. The
-// positions take none: their restraints' curvature is that of stiff bonds and angles, which move atoms together,
-// so as a diagonal it slows rather than scales the minimiser.
+// The minimiser's diagonal for the B factors of `param_atoms` (RefinementMinimiser::SetDiagonal): weight times the
+// X-ray target's curvature in B plus the B restraints', each element at least a thousandth of the median. The
+// positions take none: their restraints' curvature is that of stiff bonds and angles, which move atoms together, so
+// as a diagonal it slows rather than scales the minimiser.
 std::vector<double> BDiagonal(const std::vector<double> &xray, const std::vector<double> &restraints,
                              const std::vector<int> &param_atoms, double weight) {
    std::vector<double> h;
@@ -512,44 +558,18 @@ std::vector<double> BDiagonal(const std::vector<double> &xray, const std::vector
    return h;
 }

-// The curvature of the B factors the minimiser scales its steps by: the X-ray target's (XrayBCurvatureShape() times
-// CalibrateBCurvature()) and the B restraints' (BRestraintCurvature()), per atom.
-struct BCurvature {
-    std::vector<double> xray, restraints;
-};
-
-// `minimiser` carries its memory from the call before; `diagonal` (B only) the curvature it scales the steps by.
-void Minimise(Context &ctx, AtomTarget::Mode mode, const std::vector<int> &param_atoms, double weight,
-              int restraint_set, int iterations, RefinementMinimiser &minimiser, std::vector<double> diagonal = {}) {
-    if (param_atoms.empty())
+// `minimiser` carries its memory from the call before; `diagonal` is the curvature it scales the steps by (none: 1).
+void Minimise(Context &ctx, const AtomTarget &target, int iterations, RefinementMinimiser &minimiser,
+              std::vector<double> diagonal = {}) {
+    std::vector<double> p = target.Get();
+    if (p.empty())
        return;
-    const RefinementStage stage(mode == AtomTarget::Mode::Xyz ? "minimise.xyz" : "minimise.b");
-    std::vector<double> p;
-    for (int i : param_atoms) {
-        const gemmi::Atom &a = *ctx.atoms[i];
-        if (mode == AtomTarget::Mode::Xyz) {
-            p.push_back(a.pos.x);
-            p.push_back(a.pos.y);
-            p.push_back(a.pos.z);
-        } else
-            p.push_back(a.b_iso);
-    }
-    const AtomTarget target(ctx, mode, param_atoms, weight, restraint_set);
    minimiser.SetDiagonal(diagonal.empty() ? std::vector<double>(p.size(), 1.0) : std::move(diagonal));
    minimiser.Minimise([&](const std::vector<double> &x, double &f, std::vector<double> &g) {
                           return target.Evaluate(x.data(), &f, g.data());
                       },
-                       p, iterations, MINIMISER_F_TOLERANCE, mode == AtomTarget::Mode::Xyz ? MAX_STEP_XYZ : MAX_STEP_B,
-                       [&]() { return ctx.Cancelled(); });
-    // the atoms where the solver ended, not where its last line search probed
-    for (size_t k = 0; k < param_atoms.size(); k++) {
-        gemmi::Atom &a = *ctx.atoms[param_atoms[k]];
-        if (mode == AtomTarget::Mode::Xyz)
-            a.pos = gemmi::Position(p[3 * k], p[3 * k + 1], p[3 * k + 2]);
-        else
-            a.b_iso = static_cast<float>(std::clamp(p[k], ctx.b_floor, ctx.b_ceiling));
-    }
-    ctx.MoveRiders();
+                       p, iterations, MINIMISER_F_TOLERANCE, target.MaxStep(), [&]() { return ctx.Cancelled(); });
+    target.Set(p.data());   // where the minimiser ended, not where its last line search probed
 }

 // Groups of atoms that move as rigid bodies (one per chain) against the X-ray target: per body an angle-axis
@@ -742,6 +762,28 @@ double BalancedWeight(Context &ctx, const std::vector<int> &refinable) {
    return nx > 0 ? std::sqrt(ng / nx) : 1.0;
 }

+// The B factors' X-ray weight that balances the two gradients, as BalancedWeight() does for the positions: the B
+// similarity's gradient on the model with its B shaken by 20% (a flat B, as a model may come, gives none) over the
+// X-ray target's.
+double BalancedADPWeight(Context &ctx, const RestraintSet &b_restraints, const std::vector<int> &refinable) {
+    std::vector<std::array<double, 4>> gx;
+    ctx.engine->Evaluate(ctx.Model(), &gx);
+    ctx.FoldRiders(gx);
+    std::vector<std::array<double, 4>> xyzb = ctx.XyzB();
+    std::mt19937 rng(2);
+    std::normal_distribution<double> shake(0.0, 0.2);
+    for (std::array<double, 4> &a : xyzb)
+        a[3] *= std::exp(shake(rng));
+    std::vector<std::array<double, 4>> gr(xyzb.size(), {0, 0, 0, 0});
+    EvaluateRestraints(b_restraints, xyzb, &gr);
+    double nx = 0, ng = 0;
+    for (int i : refinable) {
+        nx += gx[i][3] * gx[i][3];
+        ng += gr[i][3] * gr[i][3];
+    }
+    return nx > 0 ? std::sqrt(ng / nx) : 1.0;
+}
+
 // The model's map coefficients and its difference map, weighted by the refinement's own likelihood: m its figure
 // of merit and D its D (2mFo - DFc, mFo - DFc).
 void MakeMaps(Context &ctx, const std::vector<MergedReflection> &merged, double wavelength_A, RefinementResult &result) {
@@ -892,6 +934,15 @@ RestraintSet RestraintsNow(Context &ctx, const RefinementTopology &topo) {
    return set;
 }

+// The B similarity of the model as it is: bonded and angle-related pairs, and atoms in contact.
+RestraintSet BRestraintsNow(Context &ctx, const RefinementTopology &topo) {
+    const RefinementStage stage("b_restraints");
+    RestraintSet set = topo.adp;
+    set.Append(BuildBContacts(topo, ctx.Model(), ctx.cell, *ctx.sg, ctx.nthreads));
+    set.Finalize(ctx.atoms.size());
+    return set;
+}
+
 // Where a weight trial runs: a model with its topology and its engine.
 struct TrialSlot {
    Context *ctx;
@@ -989,9 +1040,9 @@ TrialOutcome RunTrial(TrialSlot slot, const std::vector<int> &refinable, const s
    const RefinementStage stage(adp_w > 0 ? "rung.b" : "rung.xyz");
    ctx.SetPositions(pos0);
    ctx.SetB(b0);
-    // one memory for the positions and one for the B factors, kept from macro-cycle to macro-cycle
+    // the positions and the B factors in turn, each with a memory of its own kept from macro-cycle to macro-cycle
+    const AtomTarget xyz(ctx, AtomTarget::Mode::Xyz, refinable, w, 0), b(ctx, AtomTarget::Mode::B, refinable, adp_w, 1);
    RefinementMinimiser xyz_minimiser, b_minimiser;
-    BCurvature curvature;
    double xray_scale = 0;
    for (int cycle = 0; cycle < cycles && !ctx.Cancelled(); cycle++) {
        PlaceRiders(ctx, topo);
@@ -999,24 +1050,28 @@ TrialOutcome RunTrial(TrialSlot slot, const std::vector<int> &refinable, const s
            ctx.Restore(*fit);
        else
            ctx.Refit();
+        if (adp_w > 0)
+            ctx.MoveOverallB();
        {
            const RestraintSet set = RestraintsNow(ctx, topo);
            const RefinementStage s1("set_restraints");
            ctx.engine->SetRestraints(0, set);
        }
-        Minimise(ctx, AtomTarget::Mode::Xyz, refinable, w, 0, settings.xyz_iterations, xyz_minimiser);
+        {
+            const RefinementStage s2("minimise.xyz");
+            Minimise(ctx, xyz, settings.xyz_iterations, xyz_minimiser);
+        }
        if (adp_w > 0) {
-            {
-                const RefinementStage s2("curvature");
-                curvature.xray = ctx.XrayBCurvatureShape();
-                if (cycle == 0)
-                    xray_scale = ctx.CalibrateBCurvature(refinable, curvature.xray);
-                for (double &x : curvature.xray)
-                    x *= xray_scale;
-                curvature.restraints = BRestraintCurvature(topo.adp, ctx.XyzB());
-            }
-            Minimise(ctx, AtomTarget::Mode::B, refinable, adp_w, 1, settings.b_iterations, b_minimiser,
-                     BDiagonal(curvature.xray, curvature.restraints, refinable, adp_w));
+            const RestraintSet b_set = BRestraintsNow(ctx, topo);
+            ctx.engine->SetRestraints(1, b_set);
+            std::vector<double> xray = ctx.XrayBCurvatureShape();
+            if (cycle == 0)
+                xray_scale = ctx.CalibrateBCurvature(refinable, xray);
+            for (double &x : xray)
+                x *= xray_scale;
+            const RefinementStage s3("minimise.b");
+            Minimise(ctx, b, settings.b_iterations, b_minimiser,
+                     BDiagonal(xray, BRestraintCurvature(b_set, ctx.XyzB()), refinable, adp_w));
        }
    }
    PlaceRiders(ctx, topo);
@@ -1228,14 +1283,15 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
            RestraintSet repulsions = BuildRepulsions(topo, ctx.Model(), cell, sg, 0.0, ctx.nthreads);
            result.geometry_start = MeasureGeometry(topo.geometry, repulsions, ctx.XyzB(), &topo.metal_bonds);
        }
-        double w0 = 0;
+        double w0 = 0, adp_w0 = 0;
        {
            const RefinementStage stage("balanced_weight");
            w0 = BalancedWeight(ctx, refinable);
+            adp_w0 = BalancedADPWeight(ctx, BRestraintsNow(ctx, topo), refinable);
        }
        logger.Info("Refinement: geometry as given rmsd bonds {:.4f} A, angles {:.2f} deg, {} clashes; balanced X-ray "
-                    "weight {:.3g}", result.geometry_start.rmsd_bonds, result.geometry_start.rmsd_angles,
-                    result.geometry_start.clashes, w0);
+                    "weights {:.3g} (positions), {:.3g} (B)", result.geometry_start.rmsd_bonds,
+                    result.geometry_start.rmsd_angles, result.geometry_start.clashes, w0, adp_w0);

        const std::vector<gemmi::Position> pos0 = ctx.Positions();
        const std::vector<double> b0 = ctx.B();
@@ -1290,11 +1346,11 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
        const double w_max = w0 * settings.max_weight;

        // 1. The weights from the placed model: the ladder's rungs, tighter geometry first, each with the B factors
-        //    held and refined.
+        //    held and refined at their own balanced weight.
        std::vector<std::pair<double, double>> ladder;
        for (double rung : settings.weight_ladder)
            for (double adp : settings.adp_ladder)
-                ladder.emplace_back(w0 * rung, w0 * rung * adp);
+                ladder.emplace_back(w0 * rung, adp_w0 * adp);
        std::vector<TrialOutcome> outcomes = RunTrials(trial_slots, ladder, refinable, pos0, b0, nullptr,
                                                       settings.macro_cycles, settings);
        int k = choose(outcomes, "ladder");
@@ -1303,15 +1359,16 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
            state = std::move(outcomes[k]);

        // 2. The trajectory: steps of a few macro-cycles from where the last one ended, each trying the weights as they
-        //    are and the X-ray term's weight raised - in turn the positions' (doubled) and the B factors' (times
-        //    four) - until a step lowers neither R-work nor R-free by more than the tie.
+        //    are and, in turn, the positions' X-ray weight doubled and the B factors' divided and multiplied by four,
+        //    until a step lowers neither R-work nor R-free by more than the tie.
        for (int step = 0; step < settings.max_steps && state.ran && !ctx.Cancelled(); step++) {
            const double w = state.trial.weight, adp = state.trial.adp_weight;
            std::vector<std::pair<double, double>> trials = {{w, adp}};
-            if (step % 2 == 1 && adp > 0)
+            if (step % 2 == 1 && adp > 0) {
+                trials.emplace_back(w, adp / 4);
                trials.emplace_back(w, 4 * adp);
-            else if (2 * w <= w_max)
-                trials.emplace_back(2 * w, 2 * adp);
+            } else if (2 * w <= w_max)
+                trials.emplace_back(2 * w, adp);
            outcomes = RunTrials(trial_slots, trials, refinable, state.pos, state.b, &state.fit, settings.macro_cycles,
                                 settings);
            k = choose(outcomes, fmt::format("step {}", step + 1).c_str());
@@ -1330,7 +1387,7 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
            std::vector<std::pair<double, double>> trials;
            for (double f : settings.final_ladder)
                if (w * f <= w_max)
-                    trials.emplace_back(w * f, adp * f);
+                    trials.emplace_back(w * f, adp);
            outcomes = RunTrials(trial_slots, trials, refinable, state.pos, state.b, &state.fit, settings.macro_cycles,
                                 settings);
            k = choose(outcomes, "final");
@@ -30,11 +30,12 @@ class Logger;
 //    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
-//    stay in the structure factors and the bulk-solvent mask. The weights of the X-ray term against the
-//    geometry, and whether to refine the B factors, are chosen by R-free among trials whose geometry stays
-//    within bounds: a ladder from the placed model, then a trajectory of steps that each try the weights as
-//    they are, halved and doubled from where the last step ended until no step gains, then a short final ladder.
+//    torsions, repulsions; B similarity of bonded atoms and of atoms in contact, in proportion to B). Residues
+//    it does not describe are held fixed and stay in the structure factors and the bulk-solvent mask. The
+//    weights of the X-ray term against the geometry and against the B similarity, and whether to refine the B
+//    factors, are chosen by R-free among trials whose geometry stays within bounds: a ladder from the placed
+//    model, then a trajectory of steps that each try the weights as they are and changed, from where the last
+//    step ended, until no step gains; then a short final ladder.
 //
 // The library reads and writes nothing: the model and the reflections come in, the refined model, the
 // statistics and the map coefficients go out. Nothing it finds feeds back into processing decisions.
@@ -46,9 +47,10 @@ struct RefinementSettings {
    int xyz_iterations = 30;         // L-BFGS iterations per macro-cycle
    int b_iterations = 20;
    // The first trials, from the placed model: the positions' X-ray weight, times the gradient-balanced one
-    // (BalancedWeight), with the B factors' X-ray weight that weight times each of adp_ladder (0: B held).
+    // (BalancedWeight), each with the B factors' X-ray weight their own balanced one times each of adp_ladder (0: B
+    // held).
    std::vector<double> weight_ladder = {1.0 / 1024, 1.0 / 128, 1.0 / 16};
-    std::vector<double> adp_ladder = {0.0, 4.0};
+    std::vector<double> adp_ladder = {0.0, 1.0};
    // Then at most this many steps of the trajectory, each a set of trials from where the last one ended,
    int max_steps = 8;
    // and the final trials from where it ended, at its weights times each of these;
@@ -259,10 +259,10 @@ __global__ void restraint_kernel(const RestraintRecord *records, int n, const in
            break;
        }
        case RestraintKind::BSimilarity: {
-            const float rr = (xyzb[r.atom[0]].w - xyzb[r.atom[1]].w) / r.sigma;
-            e = rr * rr;
-            contrib[r.slot] = make_float4(0, 0, 0, 2 * rr / r.sigma);
-            contrib[r.slot + 1] = make_float4(0, 0, 0, -2 * rr / r.sigma);
+            float gi = 0, gj = 0;
+            e = refine_math::BSimilarityTerm(xyzb[r.atom[0]].w, xyzb[r.atom[1]].w, r.sigma, gi, gj);
+            contrib[r.slot] = make_float4(0, 0, 0, gi);
+            contrib[r.slot + 1] = make_float4(0, 0, 0, gj);
            break;
        }
    }
@@ -19,7 +19,7 @@
 // What a restraint is, for the evaluators on either side (RefinementRestraints.h builds them).
 enum class RestraintKind : int { Bond = 0, Angle, Torsion, Chiral, Plane, Repulsion, BSimilarity };

-// One restraint. E = ((value - ideal) / sigma)^2, summed over a set.
+// One restraint. E = ((value - ideal) / sigma)^2, summed over a set (BSimilarity: sigma is kappa, BSimilarityTerm).
 struct RestraintRecord {
    int kind = 0;
    int atom[4] = {0, 0, 0, 0};   // Plane: atom[0] the first of its atoms in plane_atoms, atom[1] their number
@@ -199,6 +199,18 @@ REFINE_HD inline T RepulsionTerm(V3<T> a, V3<T> b, T d0, T sigma, V3<T> g[2]) {
    return r * r;
 }

+// Two atoms' B alike in proportion to their size: the symmetric Kullback-Leibler divergence of two isotropic
+// Gaussians of variances B_i and B_j over kappa^2 (up to a factor), E = (B_i - B_j)^2 / (kappa^2 B_i B_j) - a
+// difference restraint whose sigma is kappa sqrt(B_i B_j). Sets dE/dB_i and dE/dB_j.
+// Following Murshudov et al. (2011) Acta Cryst. D67, 355-367 (the Kullback-Leibler form of the ADP restraints)
+template <typename T>
+REFINE_HD inline T BSimilarityTerm(T bi, T bj, T kappa, T &gi, T &gj) {
+    const T k2 = kappa * kappa, d = bi - bj;
+    gi = d * (bi + bj) / (k2 * bi * bi * bj);
+    gj = -d * (bi + bj) / (k2 * bi * bj * bj);
+    return d * d / (k2 * bi * bj);
+}
+
 // The unit eigenvector of the smallest eigenvalue of a symmetric 3x3 matrix (a00 a11 a22 a01 a02 a12), by
 // cyclic Jacobi rotations.
 template <typename T>
@@ -5,6 +5,7 @@

 #include <algorithm>
 #include <cmath>
+#include <functional>
 #include <set>
 #include <unordered_map>

@@ -25,9 +26,14 @@ constexpr float CHIRAL_SIGMA = 0.2f;                  // A^3
 constexpr float REPULSION_SIGMA = 0.2f;               // A
 constexpr float HBOND_CONTACT = 2.6f;                 // A
 constexpr float METAL_CONTACT = 1.9f;                 // A, to N, O or S
-// B similarity: main chain and the rest, for bonded pairs and pairs bonded to a common atom.
-constexpr float B_SIGMA_MAIN_BOND = 1.5f, B_SIGMA_MAIN_ANGLE = 2.0f;
-constexpr float B_SIGMA_SIDE_BOND = 3.0f, B_SIGMA_SIDE_ANGLE = 4.5f;
+// B similarity (BSimilarityTerm's kappa: sigma = kappa sqrt(B_i B_j)): main chain and the rest, for bonded pairs and
+// pairs bonded to a common atom - REFMAC's logs list sigma near 0.1 B for main-chain bonded pairs from B 11 to 140
+// A^2, the ratios between the classes are those of the fixed sigmas this replaced (1.5, 2.0, 3.0, 4.5 A^2, i.e. these
+// at B = 15) - and atoms in contact, waters among them, that no bond or angle relates.
+constexpr float B_KAPPA_MAIN_BOND = 0.1f, B_KAPPA_MAIN_ANGLE = 0.133f;
+constexpr float B_KAPPA_SIDE_BOND = 0.2f, B_KAPPA_SIDE_ANGLE = 0.3f;
+constexpr float B_KAPPA_CONTACT = 0.5f;
+constexpr double B_CONTACT_DISTANCE = 4.0;           // A

 bool IsMainChain(const gemmi::Atom &a) {
    static const std::set<std::string> names = {"N", "CA", "C", "O", "OXT", "P", "OP1", "OP2", "O5'", "C5'",
@@ -303,11 +309,11 @@ double EvaluateRestraints(const RestraintSet &set, const std::vector<std::array<
                break;
            }
            case RestraintKind::BSimilarity: {
-                const double rr = (xyzb[r.atom[0]][3] - xyzb[r.atom[1]][3]) / sigma;
-                e = rr * rr;
+                double gi = 0, gj = 0;
+                e = refine_math::BSimilarityTerm(xyzb[r.atom[0]][3], xyzb[r.atom[1]][3], sigma, gi, gj);
                if (grad != nullptr) {
-                    (*grad)[r.atom[0]][3] += 2 * rr / sigma;
-                    (*grad)[r.atom[1]][3] -= 2 * rr / sigma;
+                    (*grad)[r.atom[0]][3] += gi;
+                    (*grad)[r.atom[1]][3] += gj;
                }
                break;
            }
@@ -328,8 +334,10 @@ std::vector<double> BRestraintCurvature(const RestraintSet &set, const std::vect
    std::vector<double> h(xyzb.size(), 0.0);
    for (const RestraintRecord &r : set.records)
        if (r.kind == static_cast<int>(RestraintKind::BSimilarity)) {
-            h[r.atom[0]] += 2.0 / (r.sigma * r.sigma);
-            h[r.atom[1]] += 2.0 / (r.sigma * r.sigma);
+            // d2E/dB_i^2 of BSimilarityTerm = 2 B_j / (kappa^2 B_i^3)
+            const double bi = xyzb[r.atom[0]][3], bj = xyzb[r.atom[1]][3], k2 = double(r.sigma) * r.sigma;
+            h[r.atom[0]] += 2.0 * bj / (k2 * bi * bi * bi);
+            h[r.atom[1]] += 2.0 * bi / (k2 * bj * bj * bj);
        }
    return h;
 }
@@ -454,7 +462,7 @@ RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens)
        rb.atom[0] = i;
        rb.atom[1] = j;
        const bool main = IsMainChain(*b.atoms[0]) && IsMainChain(*b.atoms[1]);
-        rb.sigma = main ? B_SIGMA_MAIN_BOND : B_SIGMA_SIDE_BOND;
+        rb.sigma = main ? B_KAPPA_MAIN_BOND : B_KAPPA_SIDE_BOND;
        out.adp.Add(rb, 2);
    }
    for (const gemmi::Topo::Angle &an : topo->angles) {
@@ -479,7 +487,7 @@ RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens)
        rb.atom[0] = i;
        rb.atom[1] = k;
        const bool main = IsMainChain(*an.atoms[0]) && IsMainChain(*an.atoms[2]);
-        rb.sigma = main ? B_SIGMA_MAIN_ANGLE : B_SIGMA_SIDE_ANGLE;
+        rb.sigma = main ? B_KAPPA_MAIN_ANGLE : B_KAPPA_SIDE_ANGLE;
        out.adp.Add(rb, 2);
    }
    for (const gemmi::Topo::Torsion &t : topo->torsions) {
@@ -557,14 +565,17 @@ RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens)
    return out;
 }

-RestraintSet BuildRepulsions(const RefinementTopology &topo, gemmi::Model &model, const gemmi::UnitCell &cell_in,
-                             const gemmi::SpaceGroup &sg, double margin, size_t nthreads) {
+namespace {
+
+// The pairs of atoms of `model`, symmetry mates included, closer than `reach` and not excluded, where at least one of
+// the two is refinable, neither is a hydrogen and the two are not alternative conformations of each other: `make`
+// turns each into a record (the image of j as x -> rot x + tran, Cartesian, and the distance), or says to leave it out.
+// The search runs on `nthreads` threads; the list is in atom order whatever their number.
+RestraintSet Contacts(const RefinementTopology &topo, gemmi::Model &model, const gemmi::UnitCell &cell_in,
+                      const gemmi::SpaceGroup &sg, double reach, size_t nthreads,
+                      const std::function<bool(int i, int j, double d, RestraintRecord &r)> &make) {
    gemmi::UnitCell cell = cell_in;
    cell.set_cell_images_from_spacegroup(&sg);
-    float max_radius = 0;
-    for (float r : topo.vdw_radius)
-        max_radius = std::max(max_radius, r);
-    const double reach = 2 * max_radius + margin;
    gemmi::NeighborSearch ns(model, cell, std::max(reach, 5.0));
    ns.populate(false);

@@ -609,29 +620,16 @@ RestraintSet BuildRepulsions(const RefinementTopology &topo, gemmi::Model &model
                const gemmi::Mat33 rot = cell.orth.mat.multiply(ft.mat).multiply(cell.frac.mat);
                const gemmi::Vec3 tran = cell.orth.mat.multiply(ft.vec + shift);
                const gemmi::Position image(rot.multiply(pos[j]) + tran);
-                const double d = pos[i].dist(image);
-
-                float d0 = topo.vdw_radius[i] + topo.vdw_radius[j];
-                if (topo.metal[i] || topo.metal[j]) {
-                    const gemmi::El other = topo.element[topo.metal[i] ? j : static_cast<int>(i)];
-                    if (other == gemmi::El::N || other == gemmi::El::O || other == gemmi::El::S)
-                        d0 = METAL_CONTACT;
-                } else if (topo.hbond[i] && topo.hbond[j])
-                    d0 = HBOND_CONTACT;
-                if (!(d < d0 + margin))
-                    return;
                RestraintRecord r;
-                r.kind = static_cast<int>(RestraintKind::Repulsion);
                r.atom[0] = static_cast<int>(i);
                r.atom[1] = j;
-                r.ideal = d0;
-                r.sigma = REPULSION_SIGMA;
                for (int a = 0; a < 3; a++) {
                    for (int b = 0; b < 3; b++)
                        r.rot[3 * a + b] = static_cast<float>(rot[a][b]);
                    r.tran[a] = static_cast<float>(tran.at(a));
                }
-                found[block].push_back(r);
+                if (make(static_cast<int>(i), j, pos[i].dist(image), r))
+                    found[block].push_back(r);
            });
        }
    });
@@ -643,6 +641,42 @@ RestraintSet BuildRepulsions(const RefinementTopology &topo, gemmi::Model &model
    return out;
 }

+}  // namespace
+
+RestraintSet BuildRepulsions(const RefinementTopology &topo, gemmi::Model &model, const gemmi::UnitCell &cell,
+                             const gemmi::SpaceGroup &sg, double margin, size_t nthreads) {
+    float max_radius = 0;
+    for (float r : topo.vdw_radius)
+        max_radius = std::max(max_radius, r);
+    return Contacts(topo, model, cell, sg, 2 * max_radius + margin, nthreads,
+                    [&](int i, int j, double d, RestraintRecord &r) {
+                        float d0 = topo.vdw_radius[i] + topo.vdw_radius[j];
+                        if (topo.metal[i] || topo.metal[j]) {
+                            const gemmi::El other = topo.element[topo.metal[i] ? j : i];
+                            if (other == gemmi::El::N || other == gemmi::El::O || other == gemmi::El::S)
+                                d0 = METAL_CONTACT;
+                        } else if (topo.hbond[i] && topo.hbond[j])
+                            d0 = HBOND_CONTACT;
+                        if (!(d < d0 + margin))
+                            return false;
+                        r.kind = static_cast<int>(RestraintKind::Repulsion);
+                        r.ideal = d0;
+                        r.sigma = REPULSION_SIGMA;
+                        return true;
+                    });
+}
+
+RestraintSet BuildBContacts(const RefinementTopology &topo, gemmi::Model &model, const gemmi::UnitCell &cell,
+                            const gemmi::SpaceGroup &sg, size_t nthreads) {
+    return Contacts(topo, model, cell, sg, B_CONTACT_DISTANCE, nthreads, [&](int i, int j, double d, RestraintRecord &r) {
+        if (i == j || !(d < B_CONTACT_DISTANCE))
+            return false;
+        r.kind = static_cast<int>(RestraintKind::BSimilarity);
+        r.sigma = B_KAPPA_CONTACT;
+        return true;
+    });
+}
+
 void PlaceRidingHydrogens(RefinementTopology &topo) {
    gemmi::place_hydrogens_on_all_atoms(*topo.topo);
 }
@@ -96,6 +96,11 @@ void PlaceRidingHydrogens(RefinementTopology &topo);
 RestraintSet BuildRepulsions(const RefinementTopology &topo, gemmi::Model &model, const gemmi::UnitCell &cell,
                             const gemmi::SpaceGroup &sg, double margin, size_t nthreads = 1);

+// The B similarity of the atoms in contact (closer than 4 A, symmetry mates included) that no bond or angle relates,
+// waters among them, under the same conditions as BuildRepulsions().
+RestraintSet BuildBContacts(const RefinementTopology &topo, gemmi::Model &model, const gemmi::UnitCell &cell,
+                            const gemmi::SpaceGroup &sg, size_t nthreads = 1);
+
 // The usual summary of a model's geometry.
 struct GeometryStatistics {
    double rmsd_bonds = 0, rmsd_angles = 0, rmsd_chirals = 0;