diff --git a/image_analysis/structure_refinement/ModelScaleGPU.cu b/image_analysis/structure_refinement/ModelScaleGPU.cu index c5f7d054b..d77e4a6ed 100644 --- a/image_analysis/structure_refinement/ModelScaleGPU.cu +++ b/image_analysis/structure_refinement/ModelScaleGPU.cu @@ -469,17 +469,18 @@ ModelScaleParams ModelScaleGPU::Fit(const float2 *d_fcmol, const float2 *d_fmask return out; } -// FitModelScale's grid (ModelScaling.cpp, default ModelScaleBox): a coarse pass over the box and a -// finer one around its winner, every grid point its own fit, the winner the lowest finite R in grid -// order. The fits of a pass run together, one launch per Levenberg-Marquardt step. -ModelSolventFit ModelScaleGPU::FitSolvent(const float2 *d_fcmol, const float2 *d_fmask) { +// FitModelScale's grid (ModelScaling.cpp, default ModelScaleBox with b_sol up to b_hi): a coarse pass +// over the box and a finer one around its winner, every grid point its own fit, the winner the lowest +// finite R in grid order. The fits of a pass run together, one launch per Levenberg-Marquardt step. +ModelSolventFit ModelScaleGPU::FitSolvent(const float2 *d_fcmol, const float2 *d_fmask, double b_hi) { ModelSolventFit out; if (n_ < 20) return out; if (n_strong_ <= 5) throw ModelScaleGPUTooFewReflections(); - const double k_lo = 0.10, k_hi = 0.60, b_lo = 10.0, b_hi = 80.0; // ModelScaleBox{} + const double k_lo = 0.10, k_hi = 0.60, b_lo = 10.0; // ModelScaleBox{} + const double b_step = std::max(10.0, (b_hi - b_lo) / 7.0); double best_r = -1; auto try_points = [&](std::vector &fits) { FitBatch(d_fcmol, d_fmask, fits); @@ -505,15 +506,15 @@ ModelSolventFit ModelScaleGPU::FitSolvent(const float2 *d_fcmol, const float2 *d std::vector coarse; for (double ks = k_lo; ks <= k_hi + 1e-9; ks += 0.05) - for (double bs = b_lo; bs <= b_hi + 1e-9; bs += 10.0) + for (double bs = b_lo; bs <= b_hi + 1e-9; bs += b_step) coarse.push_back(point(ks, bs)); try_points(coarse); const double k0 = out.k_sol, b0 = out.b_sol; const double k_hi2 = std::min(k_hi, k0 + 0.05); - const double b_hi2 = std::min(b_hi, b0 + 10.0); + const double b_hi2 = std::min(b_hi, b0 + b_step); std::vector fine; for (double ks = std::max(k_lo, k0 - 0.05); ks <= k_hi2 + 1e-9; ks += 0.025) - for (double bs = std::max(b_lo, b0 - 10.0); bs <= b_hi2 + 1e-9; bs += 5.0) + for (double bs = std::max(b_lo, b0 - b_step); bs <= b_hi2 + 1e-9; bs += b_step / 2.0) fine.push_back(point(ks, bs)); try_points(fine); out.r = best_r; diff --git a/image_analysis/structure_refinement/ModelScaleGPU.h b/image_analysis/structure_refinement/ModelScaleGPU.h index f0f4bd7bc..3906b5c53 100644 --- a/image_analysis/structure_refinement/ModelScaleGPU.h +++ b/image_analysis/structure_refinement/ModelScaleGPU.h @@ -56,10 +56,10 @@ public: // b_star = 0. d_fcmol, d_fmask: one float2 per point on the device. ModelScaleParams Fit(const float2 *d_fcmol, const float2 *d_fmask, double k_sol, double b_sol); - // FitModelScale (ModelScaling.cpp) with the default box: the same grid, the same winner. Throws - // ModelScaleGPUTooFewReflections where FitModelScale would chain the grid points (five or fewer reflections for the - // isotropic fit). - ModelSolventFit FitSolvent(const float2 *d_fcmol, const float2 *d_fmask); + // FitModelScale (ModelScaling.cpp) with the default box, b_sol up to b_hi: the same grid, the same + // winner. Throws ModelScaleGPUTooFewReflections where FitModelScale would chain the grid points + // (five or fewer reflections for the isotropic fit). + ModelSolventFit FitSolvent(const float2 *d_fcmol, const float2 *d_fmask, double b_hi = 80.0); private: struct LevMarRun; diff --git a/image_analysis/structure_refinement/ModelScaling.cpp b/image_analysis/structure_refinement/ModelScaling.cpp index ba9c42df6..84575bd22 100644 --- a/image_analysis/structure_refinement/ModelScaling.cpp +++ b/image_analysis/structure_refinement/ModelScaling.cpp @@ -376,18 +376,21 @@ ModelScaleReport FitModelScale(gemmi::Scaling &scaling, ModelScaleBox box fitted.insert(fitted.end(), pts.begin(), pts.end()); }; - // Coarse pass over the whole box, then one refinement pass around the winner. + // Coarse pass over the whole box, then one refinement pass around the winner. b_sol is stepped by + // 10 A^2, or in eight values where the box is wider than the default, so a raised ceiling costs + // no more grid points. + const double b_step = std::max(10.0, (box.b_hi - box.b_lo) / 7.0); std::vector coarse; for (double ks = box.k_lo; ks <= box.k_hi + 1e-9; ks += 0.05) - for (double bs = box.b_lo; bs <= box.b_hi + 1e-9; bs += 10.0) + for (double bs = box.b_lo; bs <= box.b_hi + 1e-9; bs += b_step) coarse.push_back(PointFit{ks, bs}); try_points(coarse); const double k0 = best_k_sol, b0 = best_b_sol; const double k_hi2 = std::min(box.k_hi, k0 + 0.05); - const double b_hi2 = std::min(box.b_hi, b0 + 10.0); + const double b_hi2 = std::min(box.b_hi, b0 + b_step); std::vector fine; for (double ks = std::max(box.k_lo, k0 - 0.05); ks <= k_hi2 + 1e-9; ks += 0.025) - for (double bs = std::max(box.b_lo, b0 - 10.0); bs <= b_hi2 + 1e-9; bs += 5.0) + for (double bs = std::max(box.b_lo, b0 - b_step); bs <= b_hi2 + 1e-9; bs += b_step / 2.0) fine.push_back(PointFit{ks, bs}); try_points(fine); @@ -399,3 +402,16 @@ ModelScaleReport FitModelScale(gemmi::Scaling &scaling, ModelScaleBox box report.r_work_fit = best_r; return report; } + +double SolventBMax(const gemmi::Model &model) { + double sum = 0; + size_t n = 0; + for (const gemmi::Chain &chain : model.chains) + for (const gemmi::Residue &residue : chain.residues) + for (const gemmi::Atom &atom : residue.atoms) { + sum += atom.b_iso; + ++n; + } + const double b_max = ModelScaleBox{}.b_hi; + return n > 0 ? std::max(b_max, 2.0 * sum / static_cast(n)) : b_max; +} diff --git a/image_analysis/structure_refinement/ModelScaling.h b/image_analysis/structure_refinement/ModelScaling.h index 9217d6cfc..039fcb64e 100644 --- a/image_analysis/structure_refinement/ModelScaling.h +++ b/image_analysis/structure_refinement/ModelScaling.h @@ -5,6 +5,7 @@ #include +#include "gemmi/model.hpp" #include "gemmi/scaling.hpp" // Fit the overall scale, the anisotropic B and the flat bulk solvent of a model to observed @@ -21,12 +22,21 @@ // from a grid search over the physical ranges, so an unphysical answer cannot be represented. // The box the solvent parameters are searched in. The defaults are the physically reasonable // range reported by Fokine & Urzhumtsev and adopted as phenix's default grid: k_sol in (0.1, 0.8) -// and b_sol in (10, 80) A^2, of which phenix searches b_sol 10-80 in steps of 5. +// and b_sol in (10, 80) A^2, of which phenix searches b_sol 10-80 in steps of 5. A model with high +// atomic B gets a higher b_sol ceiling (SolventBMax). struct ModelScaleBox { double k_lo = 0.10, k_hi = 0.60; double b_lo = 10.0, b_hi = 80.0; }; +// The b_sol ceiling for `model`: 80 A^2, or twice the model's mean atomic B where that is higher. +// b_sol is how far the solvent's edge is smeared, and the edge can be no sharper than the atoms it +// lies against: a model at B of 300 A^2 beside a solvent held at 80 gets a solvent term that reaches +// far too deep in resolution, which only a too-small k_sol can then offset. Fokine & Urzhumtsev's +// 10-80 A^2 is for models at the B of ordinary entries (their mean b_sol 46 A^2), where the ceiling +// stays at 80. +double SolventBMax(const gemmi::Model &model); + struct ModelScaleReport { int n_points = 0; // reflections the fit actually used int n_grid = 0; // grid points evaluated diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index 95038d9c5..ea927418c 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -760,7 +760,7 @@ ModelValidationResult Validate(const std::vector &merged, gemmi::Scaling scaling(ucell, sg); scaling.use_solvent = true; scaling.prepare_points(out.fmodel, fobs_work, &ms.fmask); - FitModelScale(scaling, {}, nthreads); + FitModelScale(scaling, {.b_hi = SolventBMax(ms.st.models[0])}, nthreads); scaling.scale_data(out.fmodel, &ms.fmask); // out.fmodel now holds the scaled, solvent-corrected Fmodel out.k_sol = scaling.k_sol; out.b_sol = scaling.b_sol; diff --git a/image_analysis/structure_refinement/RigidBodyGPU.cpp b/image_analysis/structure_refinement/RigidBodyGPU.cpp index ea6b855a7..5316f580b 100644 --- a/image_analysis/structure_refinement/RigidBodyGPU.cpp +++ b/image_analysis/structure_refinement/RigidBodyGPU.cpp @@ -306,7 +306,7 @@ void RigidBodyTargetGPU::HostScale() { std::complex(fm[p][0], fm[p][1]), zone_->point_fobs[p], zone_->point_sigma[p]}); } if (!solvent_fitted_) { - FitModelScale(scaling, {}, nthreads_); + FitModelScale(scaling, {.b_hi = solvent_b_max_}, nthreads_); k_sol = scaling.k_sol; b_sol = scaling.b_sol; solvent_fitted_ = true; @@ -340,7 +340,7 @@ bool RigidBodyTargetGPU::Residuals(const double q[6], double *residuals) { // overall scale and the anisotropic B at every evaluation. if (!solvent_fitted_ && !host_scale_) { try { - engine_.FitSolvent(k_sol, b_sol); + engine_.FitSolvent(solvent_b_max_, k_sol, b_sol); solvent_fitted_ = true; } catch (const ModelScaleGPUTooFewReflections &) { host_scale_ = true; diff --git a/image_analysis/structure_refinement/RigidBodyGPU.cu b/image_analysis/structure_refinement/RigidBodyGPU.cu index 6740a4710..b8e316280 100644 --- a/image_analysis/structure_refinement/RigidBodyGPU.cu +++ b/image_analysis/structure_refinement/RigidBodyGPU.cu @@ -636,11 +636,11 @@ void RigidBodyGPUEngine::DownloadPoints(std::vector> &fcmol cuda_err(cudaStreamSynchronize(e.stream)); } -void RigidBodyGPUEngine::FitSolvent(double &k_sol, double &b_sol) { +void RigidBodyGPUEngine::FitSolvent(double b_hi, double &k_sol, double &b_sol) { RigidBodyGPUEngineImpl &e = *impl_; cuda_err(cudaSetDevice(e.device)); e.PackPoints(); - const ModelSolventFit fit = e.scale->FitSolvent(e.point_fcmol, e.point_fmask); + const ModelSolventFit fit = e.scale->FitSolvent(e.point_fcmol, e.point_fmask, b_hi); k_sol = fit.k_sol; b_sol = fit.b_sol; } diff --git a/image_analysis/structure_refinement/RigidBodyGPUEngine.h b/image_analysis/structure_refinement/RigidBodyGPUEngine.h index 8c86b20a6..a390fe781 100644 --- a/image_analysis/structure_refinement/RigidBodyGPUEngine.h +++ b/image_analysis/structure_refinement/RigidBodyGPUEngine.h @@ -112,7 +112,7 @@ public: // The scale at this evaluation's Fcalc and Fmask (ModelScaleGPU): FitModelScale()'s solvent grid, and // the overall scale and anisotropic B at a fixed solvent. FitSolvent() throws ModelScaleGPUTooFewReflections where // the grid's fits would chain (too few reflections for the isotropic start), which the host then does. - void FitSolvent(double &k_sol, double &b_sol); + void FitSolvent(double b_hi, double &k_sol, double &b_sol); void FitScale(double k_sol, double b_sol, double &k_overall, double b_star[6]); // The scale's points: fcmol (Fcalc as complex) and fmask, downloaded for the host's fit. void DownloadPoints(std::vector> &fcmol, std::vector> &fmask); diff --git a/image_analysis/structure_refinement/RigidBodyRefine.cpp b/image_analysis/structure_refinement/RigidBodyRefine.cpp index dc9e71b2d..51e3c9a24 100644 --- a/image_analysis/structure_refinement/RigidBodyRefine.cpp +++ b/image_analysis/structure_refinement/RigidBodyRefine.cpp @@ -221,7 +221,8 @@ void SymmetryComposition::Compose(const gemmi::FPhiGrid &f1, double unblu }); } -RigidBodyTargetBase::RigidBodyTargetBase(const gemmi::Model &model) : base_(ModelPositions(model)) { +RigidBodyTargetBase::RigidBodyTargetBase(const gemmi::Model &model) + : solvent_b_max_(SolventBMax(model)), base_(ModelPositions(model)) { if (base_.empty()) return; for (const gemmi::Position &p : base_) @@ -335,7 +336,7 @@ bool RigidBodyTarget::Residuals(const double q[6], double *residuals) { if (scaling.points.empty()) return false; if (!solvent_fitted_) { - FitModelScale(scaling, {}, nthreads_); + FitModelScale(scaling, {.b_hi = solvent_b_max_}, nthreads_); k_sol = scaling.k_sol; b_sol = scaling.b_sol; solvent_fitted_ = true; diff --git a/image_analysis/structure_refinement/RigidBodyRefine.h b/image_analysis/structure_refinement/RigidBodyRefine.h index bbccc538f..d334164a4 100644 --- a/image_analysis/structure_refinement/RigidBodyRefine.h +++ b/image_analysis/structure_refinement/RigidBodyRefine.h @@ -152,6 +152,7 @@ public: double k_sol = 0, b_sol = 0; // the bulk solvent the zone's target was evaluated with protected: + double solvent_b_max_; // SolventBMax() of the model: the ceiling of the zone's b_sol grid std::vector base_; gemmi::Position centre_; double rms_radius_ = 0; diff --git a/tests/ModelScalingTest.cpp b/tests/ModelScalingTest.cpp index 300f97b0f..b5b5247f7 100644 --- a/tests/ModelScalingTest.cpp +++ b/tests/ModelScalingTest.cpp @@ -11,6 +11,7 @@ #include "../image_analysis/structure_refinement/ModelScaling.h" #include "gemmi/asumask.hpp" +#include "gemmi/model.hpp" #include "gemmi/scaling.hpp" #include "gemmi/symmetry.hpp" #include "gemmi/unitcell.hpp" @@ -114,6 +115,43 @@ TEST_CASE("ModelScaling stays physical when the solvent is unidentifiable") { CHECK(std::isfinite(unbounded.b_sol)); // it converges; it is simply not bounded } +// A model at very high atomic B, where the solvent edge is smeared as much as the atoms: SolventBMax +// raises the b_sol ceiling to twice the model's mean B, and the fit then finds a b_sol far above 80 +// that the default box cannot represent - which leaves it on the box edge with a worse R. +TEST_CASE("ModelScaling finds a high b_sol for a high-B model") { + gemmi::Model model(1); + model.chains.emplace_back("A"); + model.chains[0].residues.emplace_back(); + for (double b : {280.0, 320.0}) { + gemmi::Atom atom; + atom.b_iso = static_cast(b); + model.chains[0].residues[0].atoms.push_back(atom); + } + CHECK(SolventBMax(model) == Catch::Approx(600.0)); + for (auto &atom : model.chains[0].residues[0].atoms) + atom.b_iso = 20.0f; + CHECK(SolventBMax(model) == 80.0); + + const gemmi::UnitCell cell = Cell(); + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 21 21 21"); + REQUIRE(sg != nullptr); + gemmi::Scaling scaling(cell, sg); + scaling.use_solvent = true; + MakePoints(scaling, cell, /*k_overall=*/0.5, /*b_iso=*/300.0, + /*k_sol=*/0.35, /*b_sol=*/300.0, /*d_min=*/3.5, /*d_max=*/50.0); + REQUIRE(scaling.points.size() > 500); + + gemmi::Scaling capped = scaling; + const ModelScaleReport capped_report = FitModelScale(capped); + CHECK(capped.b_sol == 80.0); + + const ModelScaleReport report = FitModelScale(scaling, {.b_hi = 600.0}); + // The b_sol steps are 590/14 A^2 in the refinement pass, so this is one step of tolerance. + CHECK(scaling.b_sol == Catch::Approx(300.0).margin(43.0)); + CHECK(scaling.k_sol == Catch::Approx(0.35).margin(0.03)); + CHECK(report.r_work_fit < capped_report.r_work_fit); +} + // A densely packed crystal with no disordered solvent channels - a small molecule, say - has an // empty solvent mask, and then Fmask is zero for every reflection. The bulk-solvent contribution // k_sol * exp(-b_sol s^2) * Fmask is then identically zero WHATEVER k_sol and b_sol come out as,