From bda90dbdd58bcd4f1f0cf8e389dbf62150811478 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 25 Oct 2026 23:35:31 +0100 Subject: [PATCH] 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 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit ed292f35b482b696d81f71f1d97e0016ca05b144) --- docs/ACKNOWLEDGEMENT.md | 10 +- .../structure_refinement/Refinement.cpp | 221 +++++++++++------- .../structure_refinement/Refinement.h | 16 +- .../structure_refinement/RefinementGPU.cu | 8 +- .../structure_refinement/RefinementMath.h | 14 +- .../RefinementRestraints.cpp | 98 +++++--- .../RefinementRestraints.h | 5 + 7 files changed, 242 insertions(+), 130 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 389b7bbd2..29a2474ab 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -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: diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index c53ab5ccc..153988e24 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -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(std::clamp(atoms[i]->b_iso + b_iso, B_FLOOR, b_ceiling)); + b_min = std::min(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(std::exp(b_iso * refl.inv_d2[m] / 4.0)); + rows.fmask[m] *= static_cast(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 ¶m_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 Get() const { + std::vector 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(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> 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 BDiagonal(const std::vector &xray, const std::vector &restraints, const std::vector ¶m_atoms, double weight) { std::vector h; @@ -512,44 +558,18 @@ std::vector BDiagonal(const std::vector &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 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 ¶m_atoms, double weight, - int restraint_set, int iterations, RefinementMinimiser &minimiser, std::vector 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 diagonal = {}) { + std::vector p = target.Get(); + if (p.empty()) return; - const RefinementStage stage(mode == AtomTarget::Mode::Xyz ? "minimise.xyz" : "minimise.b"); - std::vector 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(p.size(), 1.0) : std::move(diagonal)); minimiser.Minimise([&](const std::vector &x, double &f, std::vector &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(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 &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 &refinable) { + std::vector> gx; + ctx.engine->Evaluate(ctx.Model(), &gx); + ctx.FoldRiders(gx); + std::vector> xyzb = ctx.XyzB(); + std::mt19937 rng(2); + std::normal_distribution shake(0.0, 0.2); + for (std::array &a : xyzb) + a[3] *= std::exp(shake(rng)); + std::vector> 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 &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 &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 &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 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 pos0 = ctx.Positions(); const std::vector b0 = ctx.B(); @@ -1290,11 +1346,11 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector> 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 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> 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> 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"); diff --git a/image_analysis/structure_refinement/Refinement.h b/image_analysis/structure_refinement/Refinement.h index e083b6000..f9eadd055 100644 --- a/image_analysis/structure_refinement/Refinement.h +++ b/image_analysis/structure_refinement/Refinement.h @@ -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 weight_ladder = {1.0 / 1024, 1.0 / 128, 1.0 / 16}; - std::vector adp_ladder = {0.0, 4.0}; + std::vector 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; diff --git a/image_analysis/structure_refinement/RefinementGPU.cu b/image_analysis/structure_refinement/RefinementGPU.cu index 5c3ce356c..4ba649599 100644 --- a/image_analysis/structure_refinement/RefinementGPU.cu +++ b/image_analysis/structure_refinement/RefinementGPU.cu @@ -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; } } diff --git a/image_analysis/structure_refinement/RefinementMath.h b/image_analysis/structure_refinement/RefinementMath.h index 996a96ac6..af87c6bc6 100644 --- a/image_analysis/structure_refinement/RefinementMath.h +++ b/image_analysis/structure_refinement/RefinementMath.h @@ -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 a, V3 b, T d0, T sigma, V3 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 +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 diff --git a/image_analysis/structure_refinement/RefinementRestraints.cpp b/image_analysis/structure_refinement/RefinementRestraints.cpp index 04e997305..5d4c5ae10 100644 --- a/image_analysis/structure_refinement/RefinementRestraints.cpp +++ b/image_analysis/structure_refinement/RefinementRestraints.cpp @@ -5,6 +5,7 @@ #include #include +#include #include #include @@ -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 names = {"N", "CA", "C", "O", "OXT", "P", "OP1", "OP2", "O5'", "C5'", @@ -303,11 +309,11 @@ double EvaluateRestraints(const RestraintSet &set, const std::vector BRestraintCurvature(const RestraintSet &set, const std::vect std::vector h(xyzb.size(), 0.0); for (const RestraintRecord &r : set.records) if (r.kind == static_cast(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 &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(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(RestraintKind::Repulsion); r.atom[0] = static_cast(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(rot[a][b]); r.tran[a] = static_cast(tran.at(a)); } - found[block].push_back(r); + if (make(static_cast(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(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(RestraintKind::BSimilarity); + r.sigma = B_KAPPA_CONTACT; + return true; + }); +} + void PlaceRidingHydrogens(RefinementTopology &topo) { gemmi::place_hydrogens_on_all_atoms(*topo.topo); } diff --git a/image_analysis/structure_refinement/RefinementRestraints.h b/image_analysis/structure_refinement/RefinementRestraints.h index 5e4be92f9..572dbbc7c 100644 --- a/image_analysis/structure_refinement/RefinementRestraints.h +++ b/image_analysis/structure_refinement/RefinementRestraints.h @@ -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;