mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases 90 Wiki Activity

Merge speed-trials into rc181

Browse Source 
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed 2026-10-28 03:17:45 +01:00
2 parents deaf5b7a20 44aa7a30ef
commit 849eff07c3
3 files changed
+153 -71
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

docs/CPU_DATA_ANALYSIS_REFINEMENT.md
+41 -33
View File
@@ -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 `<prefix>_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
image_analysis/structure_refinement/Refinement.cpp
+111 -37
View File
@@ -1513,6 +1513,7 @@ struct AdpHypothesis {
std::vector<TlsBody> bodies; // Tls, with the part of each atom's B they hold at the start (Context::tls_b)
std::vector<double> tls_b;
std::vector<std::pair<int, int>> 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<AdpHypothesis> 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<Me
};
// Of a set of trials, the one kept: they come in order of trust, the first kept unless a later one displaces
// it. -1: none ran.
auto choose = [&](std::vector<TrialOutcome> &outcomes, const char *stage) {
auto best_of = [&](const std::vector<TrialOutcome> &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<int>(k);
}
return best;
};
auto choose = [&](std::vector<TrialOutcome> &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<TrialWeights> 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<Me
ladder.push_back({w0 * rung, adp_w0 * adp, settings.ncs_sigma});
}
const std::vector<AtomAniso> aniso0 = ctx.aniso;
std::vector<TrialOutcome> outcomes = RunTrials(trial_slots, ladder, refinable, pos0, b0, aniso0, nullptr,
settings.macro_cycles, settings, ncs_p);
std::vector<TrialOutcome> outcomes(ladder.size());
auto run_ladder = [&](const std::function<bool(const TrialWeights &)> &these) {
std::vector<TrialWeights> t;
for (const TrialWeights &l : ladder)
if (these(l))
t.push_back(l);
std::vector<TrialOutcome> 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<Me
}
const std::vector<gemmi::Position> changed_pos = ctx.Positions();
const std::vector<double> 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<Me
ncs_p);
logger.Info("Solvent build, round {}: the solvent refined alone: {}", round, local[0].line);
TrialOutcome after = std::move(local[0]);
if (after.ran && !(after.trial.r_free < state.trial.r_free) && built.r_free < state.trial.r_free) {
logger.Info("Solvent build, round {}: kept as built, R-work {:.4f}, R-free {:.4f}", round,
built.r_work, built.r_free);
after.pos = changed_pos;
after.b = changed_b;
after.aniso = state.aniso;
after.fit = built_fit;
after.trial.r_work = built.r_work;
after.trial.r_free = built.r_free;
after.trial.geometry = state.trial.geometry;
after.trial.within_bounds = state.trial.within_bounds;
}
if (!after.ran || !(after.trial.r_free < state.trial.r_free)) {
logger.Info("Solvent build, round {}: not kept (R-free {:.4f}, before {:.4f})", round,
after.trial.r_free, state.trial.r_free);
@@ -2694,6 +2731,9 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
b.edits.push_back(e);
}
b.rounds = round;
// the solvent's B refined for the judge only: the rest goes on with the B factors as the trajectory
// had them
after.trial.adp_weight = state.trial.adp_weight;
state = std::move(after);
from_best.clear();
}
@@ -2732,12 +2772,16 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
if (k >= 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<Me
ctx.SetAniso(state.aniso);
ctx.Restore(state.fit);
PlaceRiders(ctx, topo);
state.trial.nll_free = ctx.Score().nll_free;
std::vector<AdpHypothesis> 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<Me
ctx.SetAniso(keep);
}
const double ncs_s = state.trial.ncs_sigma;
std::vector<TrialWeights> trials = {{w, adp, ncs_s}};
std::vector<const AdpHypothesis *> 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<TrialWeights> trials;
std::vector<const AdpHypothesis *> 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<int>(trials.size());
finest = static_cast<int>(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<char> 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<char> &these) {
std::vector<TrialWeights> t;
std::vector<const AdpHypothesis *> 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<TrialOutcome> 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<char> 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<int>(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]);
}
image_analysis/structure_refinement/Refinement.h
+1 -1
View File
@@ -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<double> final_ladder = {0.5, 1.0, 2.0};
std::vector<double> 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