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

rugnux --model: molecular replacement where the rigid body does not place the model

Browse Source 
One --model now does the right thing by itself:
- A prediction (AlphaFold/ModelCIF, read from the file) goes straight to molecular replacement.
- An experimental structure is placed by the rigid body as before. The null now runs on every
  --model run (no NOT_TESTED): "placed" = the fit beats its null (4.51 sd) in the answer's merge
  or in the twin hypothesis' merge. Placed: behaviour unchanged (hand and reindex may go into
  p.mtz); a subgroup the model fits and prefers is written as <prefix>_mr.
- Not placed (or no crystal cell): molecular replacement with the experimental-model path. It may
  relabel the hand into p.mtz past the anomalous veto, never reindexes; another group or the twin
  subgroup goes only to <prefix>_mr. Maps, R-free and CC come from validating the MR placement with
  every decision switched off; MODEL_FIT stays the rigid body's verdict.
- The twin hypothesis' merge is made in memory with a model where a metric re-ask promoted the
  group or the model's group has fewer rotations than the adopted one; written and reported only
  where the model chooses it.
- --mr is kept as "--model, searched for by MR straight away".
- Report: MODEL_PLACEMENT= RIGID_BODY|MR|NOT_PLACED, MODEL_DECISIONS_TAKEN gains MR_ENANTIOMORPH,
  MODEL_CHOICE_FILE_PREFIX/SPACE_GROUP; the MR keys move into 5. MODEL VALIDATION; MODEL_FIT absent
  where no rigid body ran. MR NOT_PLACED is never a warning.
- Battery wrong-model step: check (c) also fails a wrong model placed by MR.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3
This commit is contained in:
leonarski_f
2026-10-12 22:20:36 +02:00
co-authored by Claude Opus 5.5
parent f567771d56
commit d0d225cb99
13 changed files with 462 additions and 280 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/structure_refinement/ModelValidation.cpp
+18 -23
View File
@@ -537,7 +537,8 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
const std::vector<float> &report_shell_d_min,
const ModelValidationSchedule &schedule,
const gemmi::Op *twin_law,
bool rigid_body_gpu) {
bool rigid_body_gpu,
const gemmi::Structure *model) {
ModelValidationResult result;
result.model_path = model_path;
@@ -548,7 +549,7 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
// Detect, not the default: without it GEMMI picks the format from the extension and only
// falls back to the content when it does not recognise one. A model arrives named however
// whoever produced it named it, so the file itself is the better authority.
st = gemmi::read_structure_gz(model_path, gemmi::CoorFormat::Detect);
st = model != nullptr ? *model : gemmi::read_structure_gz(model_path, gemmi::CoorFormat::Detect);
expand_strict_ncs(st);
} catch (const std::exception &e) {
result.failure_reason = fmt::format("cannot read model {}: {}", model_path, e.what());
@@ -647,7 +648,8 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
// down instead - so that every dataset of one crystal form comes out in the model's convention,
// and their maps and reflection files can be compared directly. Moving the model would leave
// each dataset in whichever indexing its own run happened to pick.
if (data_cell.is_crystal()) {
// A model handed over already placed in the data's cell (molecular replacement) is on the data's axes.
if (data_cell.is_crystal() && model == nullptr) {
const gemmi::GroupOps model_gops = sg->operations();
const bool same_point_group =
data_space_group != nullptr && same_rotations(data_space_group->operations(), model_gops);
@@ -1018,19 +1020,13 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
identity_fit = std::move(indexing.identity_fit);
indexing.identity_fit.reset();
// --- is there anything for the model to decide? ---
// There are three: the space-group label, where the model asserts the other enantiomorph; the
// indexing, where a relabelling of the data fits it better; and the setting, where the model is
// written on other axes, which the files then take (Rugnux.cpp). A model that asserts none -
// already in the group and on the axes the data were merged in, and preferring their indexing -
// has made no claim
// that needs arbitrating, and the null below would be several seconds spent gating a decision
// nobody is making. That is the isomorphous case a screening campaign is made of, and it is the
// one that has to be fast. Nothing else changes: the R-factors, the maps and the placement are
// computed and reported exactly as in every other case.
const bool decision_pending = result.model_enantiomorph_candidate
|| !(indexing.op == gemmi::Op::identity())
|| !(result.change_of_basis_op == gemmi::Op::identity());
// --- does the model fit at all? ---
// Asked of every model, whether or not it claims anything about the written reflections: the
// answer is also what says whether the model was placed at all, and a model the rigid body could
// not place is searched for by molecular replacement instead (Rugnux). Only a validation whose
// verdict is taken from another one (KeepModelVerdict), or that is not wanted (a placement already
// judged elsewhere, scored for its maps), skips it.
const bool decision_pending = schedule.fit_test;
if (schedule.on_forecast) {
// Two validations at once take twice the engines, and are allowed where twice what this one plans -
// its rigid-body engines, its structure-factor engine and, where a null is coming, the null's,
@@ -1310,10 +1306,8 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
real_r_work, null_d_min, null_mean, null_sd, result.r_work_sigma,
result.model_fits ? "FITS" : "DOES NOT FIT");
} else {
logger.Info("Model validation: the model asserts no other enantiomorph for these data and "
"prefers the indexing they were merged in, so it claims nothing about the written "
"reflections and no null was run; R-work {:.4f}, R-free {:.4f}",
real.fit.r_work, real.fit.r_free);
logger.Info("Model validation: no null here - the verdict is another validation's; R-work {:.4f}, "
"R-free {:.4f}", real.fit.r_work, real.fit.r_free);
}
#ifdef JFJOCH_USE_CUDA
} catch (const RigidBodyGPUFailure &) {
@@ -1817,7 +1811,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
double wavelength_A,
const std::vector<float> &report_shell_d_min,
const ModelValidationSchedule &schedule,
const gemmi::Op *twin_law) {
const gemmi::Op *twin_law,
const gemmi::Structure *model) {
#ifdef JFJOCH_USE_CUDA
// A CUDA failure ends the validation, not the run, and nothing is moved to the CPU part way through:
// the validation reports why it did not finish, and a validation that did not finish decides nothing,
@@ -1825,7 +1820,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
try {
return Validate(merged, cell, model_path, output_prefix, logger, data_space_group,
probe_indexing_ambiguity, nthreads, wavelength_A, report_shell_d_min, schedule, twin_law,
true);
true, model);
} catch (const RigidBodyGPUFailure &e) {
ModelValidationResult failed;
failed.model_path = model_path;
@@ -1837,7 +1832,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
}
#else
return Validate(merged, cell, model_path, output_prefix, logger, data_space_group, probe_indexing_ambiguity,
nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, false);
nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, false, model);
#endif
}