diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index 588bb140b..734d27540 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/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> 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(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 pos0 = ctx.Positions(); std::vector q(6 * bodies.size(), 0.0); auto target = std::make_unique(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(bodies.size()); + before = after; + } else { + ctx.SetPositions(pos0); + before = ctx.Refit(); + } } - ctx.Refit(); - { + if (group_b) { + const std::vector b0 = ctx.B(); std::vector q(groups.size(), 0.0); auto target = std::make_unique(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 weight_ladder = {0.0078125, 0.015625, 0.03125, 0.0625}; + std::vector 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 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; diff --git a/tools/rugnux_refine_test.cpp b/tools/rugnux_refine_test.cpp index 4350d9b7a..9c1147cd1 100644 --- a/tools/rugnux_refine_test.cpp +++ b/tools/rugnux_refine_test.cpp @@ -28,6 +28,9 @@ void usage() { " -N 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 the B refinement's weights, times the positions' chosen weight\n" " --cycles macro-cycles per weight (default as in RefinementSettings)\n" " --xyz-iterations , --b-iterations L-BFGS iterations per macro-cycle\n" " --ladder 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]);