diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index c9568d994..02ddc62b9 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -308,7 +308,9 @@ are: 3. **The final ladder**: the trajectory's weights times $\{1, 2\}$, kept where one displaces it (the weight halved was never kept on 55 refinements, and taking it away changed no choice). 4. **The displacement model**, from where the weights ended (the B factors as they are: isotropic, or the model's U - held), each trial a macro-cycle longer than a step: TLS groups with a residual B - the model's own, one per chain + held), each trial as long as a step on the positions where the weights ended, held - only the displacement + parameters refine, so the trials differ in the displacement model alone; the model chosen then takes a step with + the positions refined, kept where it displaces its trial: TLS groups with a residual B - the model's own, one per chain of at least 50 refinable atoms, and the chains cut into 2, 4, 8 and 16 contiguous segments of at least 10 residues as TLSMD cuts them (by dynamic programming over the least-squares fit of an isotropic TLS model, a quadratic B(r), to the refined B factors), each cut where it gives more groups than the one before; and a U per atom, at the U diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index 2a76ce97d..475426df9 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -1562,7 +1562,7 @@ struct TrialOutcome { TrialOutcome RunTrial(TrialSlot slot, const std::vector &refinable, const std::vector &pos0, const std::vector &b0, const std::vector &aniso0, const Context::Fit *fit, double w, double adp_w, int cycles, const RefinementSettings &settings, const NcsCopies *ncs, - double ncs_sigma, const AdpHypothesis *hyp = nullptr) { + double ncs_sigma, const AdpHypothesis *hyp = nullptr, bool hold_positions = false) { Context &ctx = *slot.ctx; RefinementTopology &topo = *slot.topo; const RefinementStage stage(hyp != nullptr ? "rung.adp" : adp_w > 0 ? "rung.b" : "rung.xyz"); @@ -1614,17 +1614,17 @@ TrialOutcome RunTrial(TrialSlot slot, const std::vector &refinable, const s } if (adp_w > 0) ctx.MoveOverallB(); - { - const RestraintSet set = RestraintsNow(ctx, topo, ncs, ncs_sigma); - const RefinementStage s1("set_restraints"); - ctx.engine->SetRestraints(0, set); - } - { + if (!hold_positions) { + { + const RestraintSet set = RestraintsNow(ctx, topo, ncs, ncs_sigma); + const RefinementStage s1("set_restraints"); + ctx.engine->SetRestraints(0, set); + } const RefinementStage s2("minimise.xyz"); Minimise(ctx, xyz, settings.xyz_iterations, xyz_minimiser); } if (kind == AdpHypothesis::Kind::Tls) { - ctx.ApplyTls(); // the bodies' U at where the atoms moved + ctx.ApplyTls(); // the bodies' U at the atoms as they are const RefinementStage s4("minimise.tls"); Minimise(ctx, tls_target, settings.b_iterations, tls_minimiser); } @@ -1721,7 +1721,8 @@ std::vector RunTrials(const std::vector &trial_slots, const std::vector &refinable, const std::vector &pos0, const std::vector &b0, const std::vector &aniso0, const Context::Fit *fit, int cycles, const RefinementSettings &settings, - const NcsCopies *ncs, const std::vector &hyps = {}) { + const NcsCopies *ncs, const std::vector &hyps = {}, + bool hold_positions = false) { std::vector out(trials.size()); std::vector errors(trial_slots.size()); // the threads shared among the slots that run (the outcome does not depend on how many each has) @@ -1747,7 +1748,8 @@ std::vector RunTrials(const std::vector &trial_slots, for (size_t j = s; j < trials.size(); j = next++) if (const size_t i = order[j]; !trial_slots[s].ctx->Cancelled()) out[i] = RunTrial(trial_slots[s], refinable, pos0, b0, aniso0, fit, trials[i].w, trials[i].adp, - cycles, settings, ncs, trials[i].ncs, hyps.empty() ? nullptr : hyps[i]); + cycles, settings, ncs, trials[i].ncs, hyps.empty() ? nullptr : hyps[i], + hold_positions); } catch (...) { errors[s] = std::current_exception(); } @@ -2889,7 +2891,11 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector ran = RunTrials(trial_slots, t, refinable, state.pos, state.b, state.aniso, - &state.fit, settings.macro_cycles + 1, settings, ncs_p, h); + &state.fit, settings.macro_cycles, settings, ncs_p, h, true); for (size_t k = 0, m = 0; k < trials.size(); k++) if (these[k]) outcomes[k] = std::move(ran[m++]); @@ -3023,6 +3029,19 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector= 0) { result.adp_model = outcomes[best].trial.adp; state = std::move(outcomes[best]); + // the model chosen refined on with the positions, a step from where its trial ended, kept where it + // displaces that + AdpHypothesis chosen = *which[best]; + chosen.bodies = state.tls; + chosen.tls_b = state.tls_b; + std::vector on = RunTrials(trial_slots, {trials[best]}, refinable, state.pos, state.b, + state.aniso, &state.fit, settings.macro_cycles, settings, + ncs_p, {&chosen}); + if (on[0].ran) { + logger.Info("Refinement, displacements, positions refined: {}", on[0].line); + if (displaces(on[0].trial, state.trial)) + state = std::move(on[0]); + } } logger.Info("Refinement: displacement model {}", result.adp_model); }