diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index 0e93175e8..f4cc05e10 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -303,25 +303,30 @@ are: stays within 0.001 of the best so far, and the trajectory ends at its best step. On a deposited model, whose free set the depositor's refinement has already fitted, R-work falls while R-free rises, and the trajectory ends where it began. -3. **The final ladder**: the trajectory's weights times $\{0.5, 1, 2\}$, kept where one displaces it. -4. **The displacement model**, from where the weights ended, each trial a macro-cycle longer than a step: the B - factors refined on as they are (isotropic, or the model's U held); 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 weight balanced as $w_{B,0}$ is (on the U shaken by 20% of its B) times $\{1/4, 1, 4\}$. A trial - displaces the one before it only where its geometry is within bounds and both its R-free (beyond the tie) and the - $-\ln L$ of the free reflections are lower: the free set tests the parameters a model adds, so the choice is made - on the evidence of each run, not at a resolution (Merritt 2012). The refined U are written as ANISOU and the TLS - groups as `_pdbx_refine_tls*` and `REMARK 3` (`REFINE_ADP_MODEL`, `REFINE_ANISO_ATOMS`, `REFINE_TLS_GROUPS`). - It is decided after the solvent build (§15.10), on the model with its solvent: a water a round put in refines a U - of its own in the per-atom trial as every refinable atom does, and no TLS group holds a water. +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 + 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 + weight balanced as $w_{B,0}$ is (on the U shaken by 20% of its B) times $\{1/4, 1\}$. The cuts between one group + per chain and the finest run only where the finest fits the free set better (lower $-\ln L$) than one per chain; + on 49 refinements of 33 structures neither this nor leaving out the U weight $\times 4$ changed a choice taken. A + trial displaces the one before it (the first: where the weights ended) only where its geometry is within bounds + and both its R-free (beyond the tie) and the $-\ln L$ of the free reflections are lower: the free set tests the + parameters a model adds, so the choice is made on the evidence of each run, not at a resolution (Merritt 2012). + The refined U are written as ANISOU and the TLS groups as `_pdbx_refine_tls*` and `REMARK 3` (`REFINE_ADP_MODEL`, + `REFINE_ANISO_ATOMS`, `REFINE_TLS_GROUPS`). It is decided after the solvent build (§15.10), on the model with its + solvent: a water a round put in refines a U of its own in the per-atom trial as every refinable atom does, and no + TLS group holds a water. Where the model has copies (§15.5), the NCS restraints are one more weight: every rung of the first ladder is run -without them and then with them at $\sigma = 0.1$ Å, and every step of the trajectory also tries them changed - -switched on where they are off, and in turn tightened ($\sigma$ halved) and loosened ($\sigma$ doubled, off beyond -0.4 Å) where they are on. Listed after their twin without, the restraints are kept only where they lower R-free by -more than the tie. A model without copies runs exactly the trials it ran before. +without them, and then with them at $\sigma = 0.1$ Å at the B factors' weight the rungs without chose (NCS moves the +X-ray weight a rung up as often as not, the B factors' choice it never changed on 39 refinements), and every step of +the trajectory also tries them changed - switched on where they are off, and in turn tightened ($\sigma$ halved) and +loosened ($\sigma$ doubled, off beyond 0.4 Å) where they are on. Listed after their twin without, the restraints are +kept only where they lower R-free by more than the tie. A model without copies runs exactly the trials it ran before. $w$ never exceeds $w_0/8$. The trials of one set run side by side - on one GPU, on several, or on the CPU's threads - each from the same start and taking nothing from the one before it, so the result does not depend @@ -428,22 +433,25 @@ a site denser than any of them (a cacodylate's As, a heavier metal) is called th Per site, everything above is in `_solvent.json` (with `electrons`, `b_env`, `peak_shift_A`, `anomalous_expected`, `n7_ligands`, `n7_off_plane`); the counts are the `REFINE_SOLVENT_` keys of the report. -**The solvent build.** Within the refinement, and to 2.8 Å (phenix.refine's ordered-solvent limit), waters are -built and removed and waters and ions relabelled, as phenix.refine builds waters within its macro-cycles. Where the -weight trajectory (§15.7) ends, a round reads the maps and the audit of that model and changes its solvent; then the -solvent alone - its waters and single-atom ions, the rest held - is refined at the weights the trajectory ended at -(its B too, at its balanced weight where the trajectory held B, and 10 iterations a macro-cycle), and the round is -kept where that lowers R-free and undone otherwise. The rounds stop at one that changes nothing, at one undone, or -after two (a third gained at most 0.0004 in R-free where it was kept, on twelve sets, and was refused on the others); -the final ladder of §15.7 then refines what they kept with the rest. The atom list of a refinement is fixed - its -engines are made for it - so the waters a round may add are in the model from the start as placeholders: one water -each (as many as a fifteenth of the model's heavy atoms, at least 64) at occupancy 0, which puts no density, no mask -and no contact anywhere and is not refined until a round puts one at a peak; a water taken out goes to occupancy 0 -the same way, and a relabelled site keeps its place. A round thus costs a map, an audit and a refinement of the -solvent: on twelve refine-arm and open-arm depositions the build adds 0.2-0.9 times the refinement's time (median -0.7), where a refinement of its own per round added 1-3.7 times (median 2.3), at an R-free within 0.001 of that on -eight of them (two 0.0016-0.0028 higher, one 0.0016 lower). `--refine-score-only` builds nothing, and `--no-refine-solvent` -turns the build off (the audit is still made). +**The solvent build.** Within the refinement, and to 2.8 Å (phenix.refine's ordered-solvent limit), waters are built +and removed and waters and ions relabelled, as phenix.refine builds waters within its macro-cycles. Where the weight +trajectory (§15.7) ends, a round reads the maps and the audit of that model and changes its solvent; then the solvent +alone - its waters and single-atom ions, the rest held - is refined at the weights the trajectory ended at (its B +too, at its balanced weight where the trajectory held B, and 10 iterations a macro-cycle), and the round is kept +where that lowers R-free; where it does not but the round as built does (refining the solvent alone can lose what its +sites gained: on a deposited model 41 waters lowered R-free by 0.0016 as built and raised it by 0.0036 refined), it +is kept as built; otherwise it is undone. The rest of the model goes on with its B factors held or refined as the +trajectory had them. The rounds stop at one that changes nothing, at one undone, or after two (a third gained at most +0.0004 in R-free where it was kept, on twelve sets, and was refused on the others); the final ladder of §15.7 then +refines what they kept with the rest. The atom list of a refinement is fixed - its engines are made for it - so the +waters a round may add are in the model from the start as placeholders: one water each (as many as a fifteenth of the +model's heavy atoms, at least 64) at occupancy 0, which puts no density, no mask and no contact anywhere and is not +refined until a round puts one at a peak; a water taken out goes to occupancy 0 the same way, and a relabelled site +keeps its place. A round thus costs a map, an audit and a refinement of the solvent: on twelve refine-arm and +open-arm depositions the build adds 0.2-0.9 times the refinement's time (median 0.7), where a refinement of its own +per round added 1-3.7 times (median 2.3), at an R-free within 0.001 of that on eight of them (two 0.0016-0.0028 +higher, one 0.0016 lower). `--refine-score-only` builds nothing, and `--no-refine-solvent` turns the build off (the +audit is still made). - *Relabelled* is a water or a single-atom ion (no alternative conformer) whose best-supported species is another one, not ambiguously - a water becomes an ion, an ion water or another ion - once per site. It takes the diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index 973e05507..e8cb56cf8 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -1513,6 +1513,7 @@ struct AdpHypothesis { std::vector bodies; // Tls, with the part of each atom's B they hold at the start (Context::tls_b) std::vector tls_b; std::vector> pairs; // Anisotropic: the pairs of the U restraints + int segments = 0; // Tls: the TLSMD cut (segments per chain at most), 0 the model's own }; // One weight trial: `cycles` macro-cycles from pos0 and b0, of the positions at X-ray weight w and, where adp_w is @@ -1919,6 +1920,7 @@ std::vector AdpHypotheses(const Context &ctx, const RefinementTop std::string name = k == 1 ? fmt::format("TLS, {} groups (one per chain)", groups.size()) : fmt::format("TLS, {} groups (chains in up to {} segments)", groups.size(), k); out.push_back(TlsHypothesis(ctx, std::move(name), std::move(groups))); + out.back().segments = k; } } AdpHypothesis aniso; @@ -2377,24 +2379,27 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector &outcomes, const char *stage) { + auto best_of = [&](const std::vector &outcomes) { int best = -1; - for (size_t k = 0; k < outcomes.size(); k++) { - const TrialOutcome &o = outcomes[k]; - if (!o.ran) - continue; - logger.Info("Refinement, {}: {}", stage, o.line); - result.trials.push_back(o.trial); - if (best < 0 || displaces(o.trial, outcomes[best].trial)) + for (size_t k = 0; k < outcomes.size(); k++) + if (outcomes[k].ran && (best < 0 || displaces(outcomes[k].trial, outcomes[best].trial))) best = static_cast(k); - } return best; }; + auto choose = [&](std::vector &outcomes, const char *stage) { + for (const TrialOutcome &o : outcomes) + if (o.ran) { + logger.Info("Refinement, {}: {}", stage, o.line); + result.trials.push_back(o.trial); + } + return best_of(outcomes); + }; 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 at their own balanced weight; where there are copies, each without NCS restraints and - // then with them. + // held and refined at their own balanced weight; where there are copies, each rung without NCS restraints + // and then with them - those with them at the B factors' weight the rungs without chose, which NCS does + // not change (on 39 refinements of 21 structures with copies it never did; the X-ray weight it does). std::vector ladder; for (double rung : settings.weight_ladder) for (double adp : settings.adp_ladder) { @@ -2403,8 +2408,23 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector aniso0 = ctx.aniso; - std::vector outcomes = RunTrials(trial_slots, ladder, refinable, pos0, b0, aniso0, nullptr, - settings.macro_cycles, settings, ncs_p); + std::vector outcomes(ladder.size()); + auto run_ladder = [&](const std::function &these) { + std::vector t; + for (const TrialWeights &l : ladder) + if (these(l)) + t.push_back(l); + std::vector ran = RunTrials(trial_slots, t, refinable, pos0, b0, aniso0, nullptr, + settings.macro_cycles, settings, ncs_p); + for (size_t i = 0, m = 0; i < ladder.size(); i++) + if (these(ladder[i])) + outcomes[i] = std::move(ran[m++]); + }; + run_ladder([](const TrialWeights &l) { return l.ncs == 0; }); + if (const int b = best_of(outcomes); ncs_p != nullptr && b >= 0) { + const double adp = outcomes[b].trial.adp_weight; + run_ladder([&](const TrialWeights &l) { return l.ncs > 0 && l.adp == adp; }); + } int k = choose(outcomes, "ladder"); TrialOutcome state; if (k >= 0) @@ -2645,6 +2665,11 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector changed_pos = ctx.Positions(); const std::vector changed_b = ctx.B(); + // the round as built, nothing refined: where refining the solvent alone loses what the sites gained, + // it is judged as built + ctx.Restore(state.fit); + const RFactors built = ctx.Refit(); + const Context::Fit built_fit = ctx.Fitted(); // (B refined, at its balanced weight where the trajectory held B: a new water's starting B is a // guess; fewer iterations, the waters' few parameters converging fast) RefinementSettings solvent_settings = settings; @@ -2657,6 +2682,18 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector= 0 && displaces(outcomes[k].trial, state.trial)) state = std::move(outcomes[k]); } - // 4. The displacement model, decided late, on the free set: from where the weights ended, the B factors refined - // on as they are, beside TLS bodies (the model's own groups, one per chain, the chains in TLSMD segments) with a - // residual B, and a U per atom at three weights - each trial a macro-cycle longer than a step. One displaces - // the one before it only where its geometry is within bounds and both its R-free (beyond the tie) and the -ln L - // of the free reflections are lower: the free set is the test of the parameters it adds. + // 4. The displacement model, decided late, on the free set: from where the weights ended, TLS bodies (the + // model's own groups, one per chain, the chains in TLSMD segments) with a residual B, and a U per atom at + // its balanced weight and a quarter of it - each trial a macro-cycle longer than a step. One displaces the one before it (the first + // where the weights ended, the B factors as they are) only where its geometry is within bounds and both its + // R-free (beyond the tie) and the -ln L of the free reflections are lower: the free set is the test of the + // parameters it adds. // Following Merritt (2012) Acta Cryst. D68, 468-477 (the displacement model the data support, decided by test) + // The trials that would not change the choice are left out: the TLS cuts between one body per chain and the + // finest are tried only where the finest fits the free set better than one body per chain, and the U at four + // times its balanced weight not at all (on 49 refinements of 33 structures, no choice taken changed). if (settings.adp_hypotheses && state.ran && !ctx.Cancelled()) { const RefinementStage adp_stage("adp_hypotheses"); ctx.SetPositions(state.pos); @@ -2745,9 +2789,8 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector hyps = AdpHypotheses(ctx, topo, input_tls); - AdpHypothesis iso; - iso.name = result.adp_input == "isotropic" ? "isotropic" : "the model's U held"; const double w = state.trial.weight, adp = state.trial.adp_weight > 0 ? state.trial.adp_weight : adp_w0; // the U weight balanced on the model with every refinable atom given a U double u_w0 = adp; @@ -2767,37 +2810,68 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector trials = {{w, adp, ncs_s}}; - std::vector which = {&iso}; + // every trial in order of trust; those between one TLS body per chain and the finest cut run after the + // others where those ask for them + std::vector trials; + std::vector which; + int coarsest = -1, finest = -1; for (const AdpHypothesis &h : hyps) { if (h.kind == AdpHypothesis::Kind::Anisotropic) { - for (double f : {0.25, 1.0, 4.0}) { + for (double f : {0.25, 1.0}) { trials.push_back({w, u_w0 * f, ncs_s}); which.push_back(&h); } } else { + if (h.segments > 0) { + if (coarsest < 0) + coarsest = static_cast(trials.size()); + finest = static_cast(trials.size()); + } trials.push_back({w, adp, ncs_s}); which.push_back(&h); } } - outcomes = RunTrials(trial_slots, trials, refinable, state.pos, state.b, state.aniso, &state.fit, - settings.macro_cycles + 1, settings, ncs_p, which); - int best = -1; + std::vector later(trials.size(), 0); + for (int k = coarsest + 1; k < finest; k++) + later[k] = 1; + outcomes.assign(trials.size(), TrialOutcome()); + auto run = [&](const std::vector &these) { + std::vector t; + std::vector h; + for (size_t k = 0; k < trials.size(); k++) + if (these[k]) { + t.push_back(trials[k]); + h.push_back(which[k]); + } + std::vector ran = RunTrials(trial_slots, t, refinable, state.pos, state.b, state.aniso, + &state.fit, settings.macro_cycles + 1, settings, ncs_p, h); + for (size_t k = 0, m = 0; k < trials.size(); k++) + if (these[k]) + outcomes[k] = std::move(ran[m++]); + }; + std::vector now(trials.size()); + for (size_t k = 0; k < trials.size(); k++) + now[k] = !later[k]; + run(now); + auto fits_better = [&](int a, int b) { + return outcomes[a].ran && outcomes[b].ran && outcomes[a].trial.nll_free < outcomes[b].trial.nll_free; + }; + if (finest > coarsest + 1 && fits_better(finest, coarsest)) + run(later); + int best = -1; // where the weights ended + auto trial = [&](int k) -> const RefinementWeightTrial & { return k < 0 ? state.trial : outcomes[k].trial; }; for (size_t k = 0; k < outcomes.size(); k++) { - const TrialOutcome &o = outcomes[k]; - if (!o.ran) + const RefinementWeightTrial &t = outcomes[k].trial; + if (!outcomes[k].ran) continue; - logger.Info("Refinement, displacements: {}", o.line); - result.adp_trials.push_back(o.trial); - if (best < 0 || (o.trial.within_bounds && (!outcomes[best].trial.within_bounds || - (o.trial.r_free < outcomes[best].trial.r_free - R_FREE_TIE && - o.trial.nll_free < outcomes[best].trial.nll_free)))) + logger.Info("Refinement, displacements: {}", outcomes[k].line); + result.adp_trials.push_back(t); + if (t.within_bounds && (!trial(best).within_bounds || + (t.r_free < trial(best).r_free - R_FREE_TIE && t.nll_free < trial(best).nll_free))) best = static_cast(k); } - // Each is a step further on the same weights: the one chosen is kept where it displaces where the weights - // ended, as any step (at low resolution the B factors refined on may all be worse than B held). - result.adp_model = iso.name; - if (best >= 0 && displaces(outcomes[best].trial, state.trial)) { + result.adp_model = result.adp_input == "isotropic" ? "isotropic" : "the model's U held"; + if (best >= 0) { result.adp_model = outcomes[best].trial.adp; state = std::move(outcomes[best]); } diff --git a/image_analysis/structure_refinement/Refinement.h b/image_analysis/structure_refinement/Refinement.h index 0bd776782..cb554f74f 100644 --- a/image_analysis/structure_refinement/Refinement.h +++ b/image_analysis/structure_refinement/Refinement.h @@ -57,7 +57,7 @@ struct RefinementSettings { // 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; - std::vector final_ladder = {0.5, 1.0, 2.0}; + std::vector final_ladder = {1.0, 2.0}; double max_weight = 1.0 / 8; // the positions' X-ray weight at most, times the balanced one bool riding_hydrogens = true; // the dictionary's hydrogens added, riding on their parents // Where the model has copies of a molecule (RefinementNcs.h), the trials try them restrained to each other, at