diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index cc7bc70d5..0f176560d 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -203,14 +203,14 @@ placement, which are statements about the model and cannot corrupt a reflection. therefore leaves the written files byte for byte what a run with no model would have produced, and the report says so beside the R it was rejected on: a negative result, not a failed run. -**The null is only built where a decision is actually pending.** A model can change exactly two things, -and it does not always claim either: one already in the space group the data were merged in, on a -crystal with no merohedral ambiguity, asserts no enantiomorph and prefers the data's own indexing, so -there is nothing to arbitrate and nothing for a null to gate. That case reports `MODEL_FIT= NOT_TESTED` -and skips the replicates. It is not an edge case — it is the isomorphous fragment-screening run, the -one that has to deliver a map seconds after the last image — and it is why the ordering matters: the -indexing probe is a handful of scaling fits and runs first anyway, so whether the expensive part is -worth paying for is known before it starts. +**The null is built for every model the rigid body is tried on.** A model already in the data's group, +on a crystal with no merohedral ambiguity, claims nothing about the written reflections, and the null +once skipped it. It no longer does, because its verdict has a second use: it is what says whether the +rigid body placed the model at all, and a model it did not place - a deposited structure of another +crystal form, a model in a subgroup the data were merged above - is then searched for by molecular +replacement, which may only relabel the enantiomorph of the answer (past the anomalous veto) and writes +any other group it chooses beside it as `_mr`. The cost is the null's one to two seconds on every +run with a model. The count is nine, and it is the *spread* that sets it rather than the mean: the verdict is (mean − real)/sd of the sample, the relative error on an sd from $n$ draws is $1/\sqrt{2(n-1)}$, and a diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 3bcb02528..20e38ed2d 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -109,8 +109,8 @@ The header block above section 1 records **how the result was produced**: `RUGNU Rates, per-image costs and progress remain on stdout only. Sections, in order: the `SUMMARY`, then `1. DATA SET AND GEOMETRY`, `2. CRYSTAL`, -`3. MERGED DATA`, `4. DIAGNOSTICS`, `5. MODEL VALIDATION` only with `--model`, and the unnumbered -`MOLECULAR REPLACEMENT (EXPERIMENTAL)` only with `--mr` ([below](#molecular-replacement---mr-experimental)) +`3. MERGED DATA`, `4. DIAGNOSTICS` and `5. MODEL VALIDATION` only with `--model` or `--mr` - which +also holds molecular replacement's keys where it ran ([below](#molecular-replacement)) (see [The summary](#the-summary-and-what-the-run-decided) below). The numbering is contiguous, and the fifth section appearing renumbers nothing; a stage that did not run states that inside its section — `MERGE= NOT_PERFORMED` — rather than the section disappearing. @@ -857,24 +857,22 @@ set (the placement's six parameters are), but because every null replicate is pl what they buy is bought on both sides and cancels; and it is decided on an order of magnitude more reflections than R-free. -That null is **only built where the model claims one of the two things it could change** — the -enantiomorph, or an indexing other than the one the data were merged in. A model already in the data's -space group on a crystal with no merohedral ambiguity, which is the isomorphous case a screening -campaign is made of, claims neither: `MODEL_FIT= NOT_TESTED`, `MODEL_DECISIONS_TAKEN= NONE`, and the -run does not pay for a null that would gate nothing. **`NOT_TESTED` is not `REJECTED`** — it says the -question was never put, not that the data answered it badly — and the three values are distinguishable -by grepping the one key. The `MODEL_FIT_NULL_*` and `MODEL_FIT_SIGMA` keys are absent in that case, -since there is no null to report; `R_WORK`, `R_FREE`, the maps and the rigid-body shift are all there -as usual. +That null is built for **every** model the rigid body is tried on, whether or not it claims anything +about the written reflections, because its answer is also whether the model was placed at all +(`MODEL_PLACEMENT=`): a model the rigid body does not place against its null - in the answer's merge or +in the twin hypothesis' one - is searched for by [molecular replacement](#molecular-replacement) +instead. `MODEL_FIT` is the rigid body's verdict and only ever that; it is absent where no rigid body was +tried (a prediction, `--mr`) or its null could not be built. | key | meaning | |---|---| | `MODEL_VALIDATION` | `PERFORMED` \| `NOT_PERFORMED` (with `MODEL_VALIDATION_REASON`, and no R-factors) | -| `MODEL_FIT` | `ACCEPTED` \| `REJECTED` \| `NOT_TESTED` — whether the model may decide anything, or had nothing to decide | +| `MODEL_PLACEMENT` | `RIGID_BODY` (the model fits against its null in the answer's merge or the twin hypothesis' one) \| `MR` (placed by molecular replacement) \| `NOT_PLACED` | +| `MODEL_FIT` | `ACCEPTED` \| `REJECTED` — the rigid body's verdict: whether the model may decide anything. Absent where no rigid body was tried | | `MODEL_FIT_STATISTIC` | What the verdict was taken on; `R_WORK` | -| `MODEL_FIT_VALUE`, `MODEL_FIT_NULL_MEAN`, `MODEL_FIT_NULL_SD`, `MODEL_FIT_NULL_REPLICATES` | The real fit, and the null of the same model in random orientations. Absent when `NOT_TESTED` | -| `MODEL_FIT_SIGMA` | The real fit above that null, in its sample standard deviations. Signed. `ACCEPTED` from +4.51: 3 σ of the null, allowing for its mean and spread being estimated from nine draws (Student's t). Absent when `NOT_TESTED` | -| `MODEL_DECISIONS_TAKEN` | `NONE` \| `ENANTIOMORPH` \| `INDEXING` \| `ENANTIOMORPH+INDEXING` | +| `MODEL_FIT_VALUE`, `MODEL_FIT_NULL_MEAN`, `MODEL_FIT_NULL_SD`, `MODEL_FIT_NULL_REPLICATES` | The real fit, and the null of the same model in random orientations. Absent with `MODEL_FIT` | +| `MODEL_FIT_SIGMA` | The real fit above that null, in its sample standard deviations. Signed. `ACCEPTED` from +4.51: 3 σ of the null, allowing for its mean and spread being estimated from nine draws (Student's t). Absent with `MODEL_FIT` | +| `MODEL_DECISIONS_TAKEN` | `NONE` \| `ENANTIOMORPH` \| `INDEXING` \| `ENANTIOMORPH+INDEXING` (the rigid body) \| `MR_ENANTIOMORPH` (molecular replacement, which never reindexes) | | `MODEL_ENANTIOMORPH_ADOPTED`, `MODEL_INDEXING_OPERATOR` | The two decisions individually; `x,y,z` is no reindexing | | `MODEL_CHANGE_OF_BASIS`, `MODEL_SETTING_AS_READ` | Present only where the model was written in another description of the lattice (other axes, other centring, or another point group) and was put into the data's to be scored. Where such a model fits, the data are then written in the model's setting instead (`SETTING_SOURCE=MODEL`) and these keys are absent, since the model no longer moves. An alternative indexing of the same point group is settled here only where a reference MTZ (`-z`) has already fixed the data's indexing; otherwise the data are reindexed into the model's, and that is `MODEL_INDEXING_OPERATOR` | | `MODEL_INDEXING_MARGIN`, `MODEL_INDEXING_MARGIN_NULL`, `MODEL_INDEXING_MARGIN_SIGMA` | Present only where a merohedral ambiguity was probed. The winner's lead over the runner-up in R-free, against the lead a random placement of the same model produces | @@ -885,7 +883,12 @@ come from whichever merge the model fits better on R-free**. `MAPS_FROM=` says w `TWIN_SUBGROUP_MERGE`), beside `TWIN_SUBGROUP_R_WORK`, `TWIN_SUBGROUP_R_FREE`, and, as twinned, `TWIN_SUBGROUP_R_WORK_TWINNED`, `TWIN_SUBGROUP_R_FREE_TWINNED` and `TWIN_SUBGROUP_MODEL_TWIN_FRACTION`. `R_WORK`/`R_FREE` above stay the answer's merge's. The model chooses the maps only: it never changes -the space group or `.mtz`. +the space group or `.mtz`. Where it also fits there against its null, it has chosen the subgroup, +and that merge is written once more as `_mr.mtz`/`.cif` with the model as placed there +(`MODEL_CHOICE_FILE_PREFIX=` and `MODEL_CHOICE_SPACE_GROUP=`, by either route). With a model, the +twin hypothesis' merge is also made where a metric re-ask promoted the indexed lattice's group, or where +the model's group has fewer rotations than the adopted one; made for that reason alone, it is written and +reported only where the model chooses it. The null's own numbers — `MODEL_FIT_STATISTIC`, `MODEL_FIT_VALUE`, the `MODEL_FIT_NULL_*` keys and the `MODEL_INDEXING_MARGIN*` keys — are written with `--developer`; the default report carries the @@ -965,14 +968,21 @@ A *significantly negative* correlation in a shell is worth chasing rather than i signal, so it points at a systematic error — an indexing the model disagrees with, or an outer shell the scaling has mistreated. -`MODEL_DECISIONS_TAKEN= NONE` — whether the model was rejected or never tested — means the reflection +`MODEL_DECISIONS_TAKEN= NONE` means the reflection files are **byte for byte** what a run with no model would have written — same space group, same indexing, same `.mtz`, `.cif`, `.hkl` and `_unmerged.mtz`. A rejected model is therefore safe to try: it costs the null's compute and changes nothing else. -## Molecular replacement (`--mr`, experimental) +## Molecular replacement -`--mr ` takes a model and places it by molecular replacement. What kind of model it is is read +A `--model` is placed by molecular replacement where it has no placement of its own - a prediction - or +where the rigid body did not place it (`MODEL_PLACEMENT=` above), and `--mr ` takes a model there +straight away, without trying the rigid body first. The search's keys (`MR_*`) and its table are in +`5. MODEL VALIDATION`. Where it places the model, the R-factors, the maps and CC(model, data) of that +section are of its placement, validated against the merge it was chosen on with nothing left to decide +(the maps `_mr_*` where that merge is `_mr`); `MODEL_FIT` stays the rigid body's. Where +it does not, they stay the rigid body's. Not placing a model is a result, never a warning. What kind of +model it is is read from the file (`MR_MODEL_KIND=`): a prediction - ModelCIF's `_ma_qa_metric`, or no cell and one B value per residue from 0 to 100 (an AlphaFold PDB or mmCIF) - has its residues below pLDDT 70 removed and its B-factors made from the pLDDT; anything else is an experimental structure, taken as its polymer only (no @@ -1052,9 +1062,9 @@ GPU where there is one (`MR_SEARCH_DEVICE=`). `MR_DECISION=` is then: is scored against the model searched for afresh in the subgroup; the subgroup is taken where its own placement leads by 0.03 in R-free and passes the subgroup's own choice (its null, its lead over the subgroup's other groups). Written as `_mr.mtz`/`.cif` from the twin hypothesis' merge, with - its law and fraction. With `--mr`, the twin hypothesis' merge is made wherever a metric re-ask + its law and fraction. With a model, the twin hypothesis' merge is made wherever a metric re-ask promoted the indexed lattice's group, so that the comparison can run - and, where that is the only - reason it was made, written and reported only where the decision is `TWIN_SUBGROUP`; + reason it was made, written and reported only where the model chooses it; - `UNDECIDED` or `NOT_PLACED` - nothing is written. With `SPACE_GROUP` and `TWIN_SUBGROUP` the placed model is written as `_mr_model.pdb`/`.cif`. diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index 0b971faa2..8af7c512e 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -537,7 +537,8 @@ ModelValidationResult Validate(const std::vector &merged, const std::vector &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 &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 &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 &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 &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 & double wavelength_A, const std::vector &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 & 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 & } #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 } diff --git a/image_analysis/structure_refinement/ModelValidation.h b/image_analysis/structure_refinement/ModelValidation.h index 42f97dc37..5bf87306f 100644 --- a/image_analysis/structure_refinement/ModelValidation.h +++ b/image_analysis/structure_refinement/ModelValidation.h @@ -149,9 +149,9 @@ struct ModelValidationResult { // sits. R-work rather than R-free: nothing is refined here, so R-work carries no optimism, and it // is decided on an order of magnitude more reflections. // - // The null is only run where the model actually claims something - the enantiomorph, or an - // indexing other than the data's own. A model that claims neither needs no arbitration, so - // fit_tested stays false and model_fits says nothing: not tested is not the same answer as failed. + // The null is run on every validation that judges a model (ModelValidationSchedule::fit_test): + // model_fits also says whether the rigid body placed the model at all. fit_tested is false only + // where the null could not be built, or where the verdict is another validation's. bool fit_tested = false; bool model_fits = false; double null_r_work_mean = 0.0; @@ -227,9 +227,7 @@ struct ModelValidationResult { // intensities), the model's group is adopted as a label, without touching the reflections - but only // where the model was shown to fit these data in the first place. A model is a hypothesis: it is // always scored, and its R-factors and maps are always reported, but it decides nothing about the -// written reflections until it has beaten a null made of itself in random orientations. That null is -// only built where the model claims one of those two things in the first place; where it claims -// neither, there is nothing to arbitrate and the run does not pay for it. +// written reflections until it has beaten a null made of itself in random orientations. // // probe_indexing_ambiguity controls how a merohedral (alternative-indexing) ambiguity is resolved. // Pass false when a reference (MTZ) already fixed the indexing - the data are then kept in that @@ -244,6 +242,9 @@ struct ModelValidationResult { // twin_law, where given, is a twin law of the data's group (a Miller-index operator in the data's // basis): the model is then also scored against the data as twinned by it (twin_target and below). // +// model, where given, is the model itself - one already placed in the data's cell, by molecular +// replacement - and model_path only names it; no other setting is then probed for it. +// // report_shell_d_min is the merge statistics' own shell bounds, coarse to fine, and it is what // cc_model_shells is binned on. Sharing the grid is the point: a reader has to be able to put a // CC(model, data) row beside that shell's CC1/2 and know the two describe the same reflections. @@ -269,6 +270,9 @@ struct ModelValidationSchedule { // Asked once the verdict, the setting and the indexing are decided, before any map is made: true // where another validation will write this one's maps and map coefficients, which are then not made. std::function maps_superseded; + // The null that says whether the model fits. False where the verdict is another validation's + // (KeepModelVerdict) or was reached elsewhere (a molecular-replacement placement, scored for its maps). + bool fit_test = true; }; ModelValidationResult ValidateAgainstModel(const std::vector &merged, @@ -282,7 +286,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector & double wavelength_A = 0.0, const std::vector &report_shell_d_min = {}, const ModelValidationSchedule &schedule = {}, - const gemmi::Op *twin_law = nullptr); + const gemmi::Op *twin_law = nullptr, + const gemmi::Structure *model = nullptr); // Reindex `merged` into the frame ValidateAgainstModel reported, so the reflection files that are // written describe the same indexing as the R-factors and the maps. Returns the space group they are diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index f90a6bf6e..d2bb0ef7a 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -911,9 +911,15 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, // ASSUMED, not determined: nothing here measured the hand. A model that fits carries prior // chemical knowledge these intensities do not - and cannot, since |Fcalc| is invariant // under the change of hand - so taking its group is an assertion, not a measurement. - if (mv.has_value() && mv->ok && mv->model_fits - && (mv->adopted_model_enantiomorph - || mv->model_space_group_number == result.space_group->number)) + // Molecular replacement tells the two groups apart by where it places the copies, as surely as + // the rigid body does by its fit: where it placed the model in the group written, the same. + const bool mr_placed_here = result.model_placement == "MR" + && result.model_prior->candidates[result.model_prior->chosen].space_group.number + == result.space_group->number; + if ((mv.has_value() && mv->ok && mv->model_fits + && (mv->adopted_model_enantiomorph + || mv->model_space_group_number == result.space_group->number)) + || mr_placed_here) enantiomorph = "ASSUMED_FROM_MODEL"; else if (!result.space_group_search.has_value()) enantiomorph = "GIVEN"; // -S, or a reference MTZ: the user's assertion @@ -2318,17 +2324,29 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, } // -------------------------------------------------- 5. MODEL VALIDATION - if (result.model_validation.has_value()) { - const auto &mv = *result.model_validation; + if (result.model_validation.has_value() || result.model_prior.has_value()) { + // A prediction molecular replacement did not place has no validation at all: read as one not performed. + ModelValidationResult unplaced; + if (!result.model_validation) { + const auto &mp = *result.model_prior; + unplaced.model_path = mp.model_path; + unplaced.failure_reason = mp.ok ? "molecular replacement did not place the model" + : "the model could not be searched for: " + mp.failure_reason; + } + const auto &mv = result.model_validation ? *result.model_validation : unplaced; + const bool by_mr = result.model_placement == "MR"; + const bool mr_hand = by_mr && result.model_prior->hand_check.adopted; ReportSection s; s.title = "5. MODEL VALIDATION"; Add(s, KeyText("MODEL_FILE", mv.model_path)); + if (!result.model_placement.empty()) + Add(s, KeyEnum("MODEL_PLACEMENT", result.model_placement, {"RIGID_BODY", "MR", "NOT_PLACED"})); if (!mv.ok) { Add(s, KeyEnum("MODEL_VALIDATION", "NOT_PERFORMED", {"PERFORMED", "NOT_PERFORMED"})); Add(s, KeyText("MODEL_VALIDATION_REASON", mv.failure_reason)); Add(s, Blank()); - Add(s, Prose(" A model was given but could not be used, so there are no R-factors and no maps.\n" - " Everything else in this report is unaffected: the merge does not depend on the model.")); + Add(s, Prose(" A model was given but could not be used or placed, so there are no R-factors and no\n" + " maps. Everything else in this report is unaffected: the merge does not depend on the model.")); } else { Add(s, KeyEnum("MODEL_VALIDATION", "PERFORMED", {"PERFORMED", "NOT_PERFORMED"})); Add(s, KeyInt("MODEL_SPACE_GROUP_NUMBER", mv.model_space_group_number)); @@ -2348,12 +2366,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyReal("BULK_SOLVENT_K_SOL", mv.k_sol, "{:.3f}", true)); Add(s, KeyReal("BULK_SOLVENT_B_SOL", mv.b_sol, "{:.1f}", true)); Add(s, KeyReal("SCALE_OVERALL", mv.k_overall, "{:.4f}", true)); - // NOT_TESTED is a third answer and not a missing one: the model claimed nothing that - // needed arbitrating, so the null was not built. - Add(s, KeyEnum("MODEL_FIT", - !mv.fit_tested ? "NOT_TESTED" : (mv.model_fits ? "ACCEPTED" : "REJECTED"), - {"ACCEPTED", "REJECTED", "NOT_TESTED"})); + // The rigid body's verdict, against its null; absent where no rigid body was tried (a prediction, + // --mr) or the null could not be built. A molecular-replacement placement never makes it. if (mv.fit_tested) { + Add(s, KeyEnum("MODEL_FIT", mv.model_fits ? "ACCEPTED" : "REJECTED", {"ACCEPTED", "REJECTED"})); Add(s, KeyText("MODEL_FIT_STATISTIC", "R_WORK", true)); Add(s, KeyReal("MODEL_FIT_VALUE", mv.r_work, "{:.4f}", true)); Add(s, KeyReal("MODEL_FIT_NULL_MEAN", mv.null_r_work_mean, "{:.4f}", true)); @@ -2404,12 +2420,15 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, mv.anomalous_sites[i].sigma))); } const bool took_indexing = !(mv.indexing_op == gemmi::Op::identity()); + // What the model wrote into the answer's files, and by which route: the rigid body may take the + // enantiomorph and the indexing, molecular replacement only the enantiomorph. Add(s, KeyEnum("MODEL_DECISIONS_TAKEN", - mv.adopted_model_enantiomorph + mr_hand ? "MR_ENANTIOMORPH" + : mv.adopted_model_enantiomorph ? (took_indexing ? "ENANTIOMORPH+INDEXING" : "ENANTIOMORPH") : (took_indexing ? "INDEXING" : "NONE"), - {"NONE", "INDEXING", "ENANTIOMORPH", "ENANTIOMORPH+INDEXING"})); - Add(s, KeyBool("MODEL_ENANTIOMORPH_ADOPTED", mv.adopted_model_enantiomorph)); + {"NONE", "INDEXING", "ENANTIOMORPH", "ENANTIOMORPH+INDEXING", "MR_ENANTIOMORPH"})); + Add(s, KeyBool("MODEL_ENANTIOMORPH_ADOPTED", mv.adopted_model_enantiomorph || mr_hand)); Add(s, KeyText("MODEL_INDEXING_OPERATOR", mv.indexing_op.triplet())); if (!(mv.change_of_basis_op == gemmi::Op::identity())) { Add(s, KeyText("MODEL_CHANGE_OF_BASIS", mv.change_of_basis_op.triplet())); @@ -2424,6 +2443,12 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, } if (!mv.maps_prefix.empty()) Add(s, KeyText("MAPS_PREFIX", mv.maps_prefix)); + // The model's choice of another group or of the twin hypothesis' subgroup, by either route: files of + // their own, never the answer's. + if (!result.model_prior_file_prefix.empty()) { + Add(s, KeyText("MODEL_CHOICE_FILE_PREFIX", result.model_prior_file_prefix)); + Add(s, KeyText("MODEL_CHOICE_SPACE_GROUP", result.model_prior_file_space_group)); + } // The twin hypothesis (section 4): the model scored against that merge too, and which of // the two the maps and the placed model were made from. The keys above are the answer's. if (result.twin_subgroup_merge && result.twin_subgroup_merge->model_validation) { @@ -2536,12 +2561,12 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, " MODEL_FIT says whether these data accept the model: the same model was refitted from\n" " {} random orientations, and MODEL_FIT_SIGMA ({:+.2f}) is how far the real fit sits\n" " above that null.", mv.null_replicates, mv.r_work_sigma))); - else - Add(s, Prose(" MODEL_FIT= NOT_TESTED is not a failed test: the model asserted neither another hand\n" - " nor another indexing, so there was nothing to arbitrate and the null was not built.\n" - " The R-factors, the maps and the rigid-body placement are exactly what they would\n" - " have been.")); - if (mv.fit_tested && !mv.model_fits) + if (by_mr) + Add(s, Prose(" The model was placed by molecular replacement (MODEL_PLACEMENT, and the MR_ keys below):\n" + " the R-factors, the maps and CC(model, data) above are of that placement, validated\n" + " against the merge it was chosen on with nothing left to decide. MODEL_FIT, where\n" + " present, is the rigid body's verdict on the model as it was read, never the search's.")); + if (mv.fit_tested && !mv.model_fits && !by_mr) Add(s, Prose(fmt::format( " The model was tried and REJECTED - R-work {:.4f} against {:.4f} +- {:.4f} for the same\n" " model in random orientations, which is no better than chance. Nothing downstream moved:\n" @@ -2571,123 +2596,118 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, "the reflections have not been reindexed, which would have hidden which of the two it is", mv.anomalous_deepest_site, mv.anomalous_deepest_sigma)); } - doc.sections.push_back(std::move(s)); - } - - // ------------------------------------------------- MOLECULAR REPLACEMENT (--mr, experimental) - if (result.model_prior.has_value()) { - const auto &mp = *result.model_prior; - ReportSection s; - s.title = "MOLECULAR REPLACEMENT (EXPERIMENTAL)"; - Add(s, KeyText("MR_MODEL_FILE", mp.model_path)); - if (!mp.ok) { - Add(s, KeyEnum("MR", "NOT_PERFORMED", {"PERFORMED", "NOT_PERFORMED"})); - Add(s, KeyText("MR_REASON", mp.failure_reason)); - } else { - Add(s, KeyEnum("MR", "PERFORMED", {"PERFORMED", "NOT_PERFORMED"})); - Add(s, KeyEnum("MR_MODEL_KIND", mp.predicted_model ? "PREDICTED" : "EXPERIMENTAL", - {"PREDICTED", "EXPERIMENTAL"})); - Add(s, KeyInt("MR_RESIDUES_READ", mp.residues_read)); - Add(s, KeyInt("MR_RESIDUES_KEPT", mp.residues_kept)); - Add(s, KeyInt("MR_PAE_DOMAINS", mp.pae_domains)); - Add(s, KeyText("MR_SEARCH_MODEL", mp.search_model)); - Add(s, KeyReal("MR_SEARCH_D_MIN", mp.d_min_search, "{:.2f}")); - Add(s, KeyReal("MR_SCORE_D_MIN", mp.d_min_score, "{:.2f}")); - Add(s, KeyText("MR_DECISION", mp.decision)); - if (mp.chosen >= 0) - { - Add(s, KeyText("MR_SPACE_GROUP", mp.candidates[mp.chosen].space_group.xhm())); - Add(s, KeyInt("MR_COPIES", mp.candidates[mp.chosen].copies)); - } - if (mp.hand_check.performed) { - const auto &h = mp.hand_check; - Add(s, KeyEnum("MR_HAND_CHECK", h.adopted ? "ADOPTED" : "VETOED", {"ADOPTED", "VETOED"})); - Add(s, KeyReal("MR_HAND_FIT_SIGMA", h.fit_sigma, "{:+.2f}")); - Add(s, KeyBool("MR_ANOMALOUS_INVERTED", h.anomalous_inverted)); - Add(s, KeyInt("MR_ANOMALOUS_PAIRS", h.anomalous_pairs)); - if (std::isfinite(h.anomalous_scatterer_mean_sigma)) - Add(s, KeyReal("MR_ANOMALOUS_SCATTERER_MEAN_SIGMA", h.anomalous_scatterer_mean_sigma, - "{:+.2f}")); - if (h.anomalous_inverted) - Add(s, KeyText("MR_ANOMALOUS_DEEPEST_SITE", - fmt::format("{} {:.1f}", h.deepest_site, h.deepest_sigma))); - } - if (!result.model_prior_file_prefix.empty()) - Add(s, KeyText("MR_FILE_PREFIX", result.model_prior_file_prefix)); - Add(s, KeyEnum("MR_SEARCH_DEVICE", mp.gpu_search ? "GPU" : "CPU", {"GPU", "CPU"})); - Add(s, KeyReal("MR_SECONDS", mp.t_model + mp.t_transform + mp.t_rotation + mp.t_translation - + mp.t_rescore + mp.t_null + mp.t_domains + mp.hand_check.seconds, "{:.2f}")); + // Molecular replacement, where it ran: for a prediction, with --mr, or for a model the rigid body did + // not place. In this section, beside the validation, since it is the same model's. + if (result.model_prior.has_value()) { + const auto &mp = *result.model_prior; Add(s, Blank()); - ReportEntry e; - e.kind = ReportEntry::Kind::Table; - e.table.columns = {"MERGE", "SPACE_GROUP", "VM", "COPIES", "COPIES_LIKELY", "RF_Z", "TF_Z", "LLG_Z", "CLASH", - "R_WORK", "R_FREE", "R_FREE_TWIN", "NULL_SIGMA", "VERDICT"}; - e.table.text_header = " MERGE SPACE_GROUP VM COPIES LIKELY RF_Z TF_Z LLG_Z CLASH R_WORK R_FREE R_FREE_TWIN NULL_SIGMA VERDICT\n" - " ------------------- ----------- ----- ------ ------ ----- ----- ----- ----- ------ ------ ----------- ---------- ----------"; - for (const auto &c : mp.candidates) { - const std::string merge = c.merge == 0 ? "MERGE" : "TWIN_SUBGROUP_MERGE"; - const auto r = [](double v) { return std::isfinite(v) ? fmt::format("{:.3f}", v) : std::string("-"); }; - e.table.text_rows.push_back(fmt::format( - " {:<19s} {:<11s} {:5.2f} {:6d} {:6d} {:5.1f} {:5.1f} {:>5s} {:4.1f}% {:>6s} {:>6s} {:>11s} {:>10s} {}", - merge, c.space_group.xhm(), c.vm, c.copies, c.copies_likely, c.rf_z, c.tf_z, - std::isfinite(c.llg_z) ? fmt::format("{:.1f}", c.llg_z) : std::string("-"), - 100 * c.clash_fraction, r(c.r_work), - r(c.r_free), r(c.r_free_twin), std::isfinite(c.null_sigma) ? fmt::format("{:+.1f}", c.null_sigma) : std::string("-"), - c.verdict)); - std::vector row; - row.push_back(KeyText("", merge).value); - row.push_back(KeyText("", c.space_group.xhm()).value); - row.push_back(KeyReal("", c.vm, "{:.2f}").value); - row.push_back(KeyInt("", c.copies).value); - row.push_back(KeyInt("", c.copies_likely).value); - row.push_back(KeyReal("", c.rf_z, "{:.1f}").value); - row.push_back(KeyReal("", c.tf_z, "{:.1f}").value); - row.push_back(KeyReal("", c.llg_z, "{:.1f}").value); - row.push_back(KeyReal("", c.clash_fraction, "{:.3f}").value); - row.push_back(KeyReal("", c.r_work, "{:.3f}").value); - row.push_back(KeyReal("", c.r_free, "{:.3f}").value); - row.push_back(KeyReal("", c.r_free_twin, "{:.3f}").value); - row.push_back(KeyReal("", c.null_sigma, "{:+.1f}").value); - row.push_back(KeyText("", c.verdict).value); - e.table.cells.push_back(std::move(row)); + if (!mp.ok) { + Add(s, KeyEnum("MR", "NOT_PERFORMED", {"PERFORMED", "NOT_PERFORMED"})); + Add(s, KeyText("MR_REASON", mp.failure_reason)); + } else { + Add(s, KeyEnum("MR", "PERFORMED", {"PERFORMED", "NOT_PERFORMED"})); + Add(s, KeyEnum("MR_MODEL_KIND", mp.predicted_model ? "PREDICTED" : "EXPERIMENTAL", + {"PREDICTED", "EXPERIMENTAL"})); + Add(s, KeyInt("MR_RESIDUES_READ", mp.residues_read)); + Add(s, KeyInt("MR_RESIDUES_KEPT", mp.residues_kept)); + Add(s, KeyInt("MR_PAE_DOMAINS", mp.pae_domains)); + Add(s, KeyText("MR_SEARCH_MODEL", mp.search_model)); + Add(s, KeyReal("MR_SEARCH_D_MIN", mp.d_min_search, "{:.2f}")); + Add(s, KeyReal("MR_SCORE_D_MIN", mp.d_min_score, "{:.2f}")); + Add(s, KeyText("MR_DECISION", mp.decision)); + if (mp.chosen >= 0) + { + Add(s, KeyText("MR_SPACE_GROUP", mp.candidates[mp.chosen].space_group.xhm())); + Add(s, KeyInt("MR_COPIES", mp.candidates[mp.chosen].copies)); + } + if (mp.hand_check.performed) { + const auto &h = mp.hand_check; + Add(s, KeyEnum("MR_HAND_CHECK", h.adopted ? "ADOPTED" : "VETOED", {"ADOPTED", "VETOED"})); + Add(s, KeyReal("MR_HAND_FIT_SIGMA", h.fit_sigma, "{:+.2f}")); + Add(s, KeyBool("MR_ANOMALOUS_INVERTED", h.anomalous_inverted)); + Add(s, KeyInt("MR_ANOMALOUS_PAIRS", h.anomalous_pairs)); + if (std::isfinite(h.anomalous_scatterer_mean_sigma)) + Add(s, KeyReal("MR_ANOMALOUS_SCATTERER_MEAN_SIGMA", h.anomalous_scatterer_mean_sigma, + "{:+.2f}")); + if (h.anomalous_inverted) + Add(s, KeyText("MR_ANOMALOUS_DEEPEST_SITE", + fmt::format("{} {:.1f}", h.deepest_site, h.deepest_sigma))); + } + Add(s, KeyEnum("MR_SEARCH_DEVICE", mp.gpu_search ? "GPU" : "CPU", {"GPU", "CPU"})); + Add(s, KeyReal("MR_SECONDS", mp.t_model + mp.t_transform + mp.t_rotation + mp.t_translation + + mp.t_rescore + mp.t_null + mp.t_domains + mp.hand_check.seconds, "{:.2f}")); + Add(s, Blank()); + ReportEntry e; + e.kind = ReportEntry::Kind::Table; + e.table.columns = {"MERGE", "SPACE_GROUP", "VM", "COPIES", "COPIES_LIKELY", "RF_Z", "TF_Z", "LLG_Z", "CLASH", + "R_WORK", "R_FREE", "R_FREE_TWIN", "NULL_SIGMA", "VERDICT"}; + e.table.text_header = " MERGE SPACE_GROUP VM COPIES LIKELY RF_Z TF_Z LLG_Z CLASH R_WORK R_FREE R_FREE_TWIN NULL_SIGMA VERDICT\n" + " ------------------- ----------- ----- ------ ------ ----- ----- ----- ----- ------ ------ ----------- ---------- ----------"; + for (const auto &c : mp.candidates) { + const std::string merge = c.merge == 0 ? "MERGE" : "TWIN_SUBGROUP_MERGE"; + const auto r = [](double v) { return std::isfinite(v) ? fmt::format("{:.3f}", v) : std::string("-"); }; + e.table.text_rows.push_back(fmt::format( + " {:<19s} {:<11s} {:5.2f} {:6d} {:6d} {:5.1f} {:5.1f} {:>5s} {:4.1f}% {:>6s} {:>6s} {:>11s} {:>10s} {}", + merge, c.space_group.xhm(), c.vm, c.copies, c.copies_likely, c.rf_z, c.tf_z, + std::isfinite(c.llg_z) ? fmt::format("{:.1f}", c.llg_z) : std::string("-"), + 100 * c.clash_fraction, r(c.r_work), + r(c.r_free), r(c.r_free_twin), std::isfinite(c.null_sigma) ? fmt::format("{:+.1f}", c.null_sigma) : std::string("-"), + c.verdict)); + std::vector row; + row.push_back(KeyText("", merge).value); + row.push_back(KeyText("", c.space_group.xhm()).value); + row.push_back(KeyReal("", c.vm, "{:.2f}").value); + row.push_back(KeyInt("", c.copies).value); + row.push_back(KeyInt("", c.copies_likely).value); + row.push_back(KeyReal("", c.rf_z, "{:.1f}").value); + row.push_back(KeyReal("", c.tf_z, "{:.1f}").value); + row.push_back(KeyReal("", c.llg_z, "{:.1f}").value); + row.push_back(KeyReal("", c.clash_fraction, "{:.3f}").value); + row.push_back(KeyReal("", c.r_work, "{:.3f}").value); + row.push_back(KeyReal("", c.r_free, "{:.3f}").value); + row.push_back(KeyReal("", c.r_free_twin, "{:.3f}").value); + row.push_back(KeyReal("", c.null_sigma, "{:+.1f}").value); + row.push_back(KeyText("", c.verdict).value); + e.table.cells.push_back(std::move(row)); + } + Add(s, std::move(e)); + Add(s, Blank()); + Add(s, Prose(fmt::format( + " The model - a prediction trimmed to pLDDT >= 70, or an experimental structure\n" + " (MR_MODEL_KIND) - was searched for in every space group each merge\n" + " cannot tell apart (rotation function on a 6 deg sampling, translation function, {:.1f}-{:.1f} A),\n" + " and the best placements taken through the rigid body; then further copies, one at a time, up to\n" + " the largest copy number with a Matthews probability of 1% (COPIES_LIKELY: the most probable),\n" + " each with the translation function beside the copies already placed, its peaks ranked by the\n" + " likelihood gain of the copy (sigma_A from the pLDDT's coordinate error - 0.67 A for an\n" + " experimental structure - its scale refined), kept\n" + " where the copy stands out, packs and R-free drops by 0.01. R_WORK and R_FREE are of all COPIES,\n" + " to MR_SCORE_D_MIN; TF_Z is the highest of the copies'; LLG_Z the highest Z-score of a copy's refined\n" + " likelihood gain over the same copy at random positions, where its TF_Z was below 8. A copy stands\n" + " out where TF_Z or LLG_Z is 8 or more. PLACED needs a copy that stands out, at most 5% of the\n" + " C-alpha atoms clashing with another copy or a symmetry mate, one copy fitting the cell\n" + " (MATTHEWS: Vm below 1.5) and, for the candidate the choice falls on, an R-work at least 3 sd\n" + " better than the same copies refitted from random orientations (NULL where it is not; the next\n" + " candidate is then taken). The model chooses within a merge where its placed group leads every\n" + " other that packs by 0.03 in R-free, and the twin hypothesis' merge where no group of the\n" + " answer's merge can hold it at all, or where, on that merge and against its data as twinned by\n" + " its law, the model searched for in the subgroup leads the chosen placement read as a crystal\n" + " there by 0.03 in R-free (R_FREE_TWIN). A group one copy cannot fit is checked for\n" + " packing only, not refined. Where the AlphaFold DB's predicted aligned error sits beside the\n" + " model, the prediction is cut into rigid domains (MR_PAE_DOMAINS), residues in none left out;\n" + " where the whole model is placed nowhere, each domain is searched by itself, and the domains\n" + " placed have to agree (MR_SEARCH_MODEL: what the decision was made with).\n" + " The answer's files take the enantiomorph from it - only where model validation of the placed\n" + " model accepts it as it accepts a --model's: a fit over its null and an anomalous map at its\n" + " atoms that is not inverted (ENANTIOMORPH_VETOED otherwise) - and nothing else; another group it\n" + " chooses is written as {}.mtz/.cif with the placed model.", + mp.d_max_search, mp.d_min_search, + result.model_prior_file_prefix.empty() ? std::string("_mr") : result.model_prior_file_prefix))); + Add(s, Prose(fmt::format( + " Seconds: model {:.2f}, transform {:.2f}, rotation {:.2f}, translation {:.2f}, rigid body and\n" + " R-free {:.2f}, null {:.2f}, PAE domains {:.2f}, hand check {:.3f}.", mp.t_model, mp.t_transform, + mp.t_rotation, mp.t_translation, mp.t_rescore, mp.t_null, mp.t_domains, mp.hand_check.seconds), true)); } - Add(s, std::move(e)); - Add(s, Blank()); - Add(s, Prose(fmt::format( - " The model - a prediction trimmed to pLDDT >= 70, or an experimental structure\n" - " (MR_MODEL_KIND) - was searched for in every space group each merge\n" - " cannot tell apart (rotation function on a 6 deg sampling, translation function, {:.1f}-{:.1f} A),\n" - " and the best placements taken through the rigid body; then further copies, one at a time, up to\n" - " the largest copy number with a Matthews probability of 1% (COPIES_LIKELY: the most probable),\n" - " each with the translation function beside the copies already placed, its peaks ranked by the\n" - " likelihood gain of the copy (sigma_A from the pLDDT's coordinate error - 0.67 A for an\n" - " experimental structure - its scale refined), kept\n" - " where the copy stands out, packs and R-free drops by 0.01. R_WORK and R_FREE are of all COPIES,\n" - " to MR_SCORE_D_MIN; TF_Z is the highest of the copies'; LLG_Z the highest Z-score of a copy's refined\n" - " likelihood gain over the same copy at random positions, where its TF_Z was below 8. A copy stands\n" - " out where TF_Z or LLG_Z is 8 or more. PLACED needs a copy that stands out, at most 5% of the\n" - " C-alpha atoms clashing with another copy or a symmetry mate, one copy fitting the cell\n" - " (MATTHEWS: Vm below 1.5) and, for the candidate the choice falls on, an R-work at least 3 sd\n" - " better than the same copies refitted from random orientations (NULL where it is not; the next\n" - " candidate is then taken). The model chooses within a merge where its placed group leads every\n" - " other that packs by 0.03 in R-free, and the twin hypothesis' merge where no group of the\n" - " answer's merge can hold it at all, or where, on that merge and against its data as twinned by\n" - " its law, the model searched for in the subgroup leads the chosen placement read as a crystal\n" - " there by 0.03 in R-free (R_FREE_TWIN). A group one copy cannot fit is checked for\n" - " packing only, not refined. Where the AlphaFold DB's predicted aligned error sits beside the\n" - " model, the prediction is cut into rigid domains (MR_PAE_DOMAINS), residues in none left out;\n" - " where the whole model is placed nowhere, each domain is searched by itself, and the domains\n" - " placed have to agree (MR_SEARCH_MODEL: what the decision was made with).\n" - " The answer's files take the enantiomorph from it - only where model validation of the placed\n" - " model accepts it as it accepts a --model's: a fit over its null and an anomalous map at its\n" - " atoms that is not inverted (ENANTIOMORPH_VETOED otherwise) - and nothing else; another group it\n" - " chooses is written as {}.mtz/.cif with the placed model.", - mp.d_max_search, mp.d_min_search, - result.model_prior_file_prefix.empty() ? std::string("_mr") : result.model_prior_file_prefix))); - Add(s, Prose(fmt::format( - " Seconds: model {:.2f}, transform {:.2f}, rotation {:.2f}, translation {:.2f}, rigid body and\n" - " R-free {:.2f}, null {:.2f}, PAE domains {:.2f}, hand check {:.3f}.", mp.t_model, mp.t_transform, - mp.t_rotation, mp.t_translation, mp.t_rescore, mp.t_null, mp.t_domains, mp.hand_check.seconds), true)); } doc.sections.push_back(std::move(s)); } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 2d6571c8c..203418910 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -200,12 +200,10 @@ struct ProcessConfig { // Fo-Fc maps (--model). std::string model_path; - // Experimental (--mr): a model - a prediction (AlphaFold, pLDDT in the B column) or an experimental - // structure, told apart from the file - placed by molecular replacement, as many copies as the cell holds, - // in every space group the merge cannot tell apart and in the twin hypothesis' subgroup. It never changes - // the answer's files beyond the enantiomorph label; a different group it chooses is written beside them - // as _mr. - std::string mr_model_path; + // --mr: the model is searched for by molecular replacement straight away, without the rigid-body + // placement first. Without it a --model that is a prediction, or an experimental structure the rigid body + // did not place, is searched for anyway. + bool force_mr = false; // --finalist-ledger: report-only symmetry evidence table (and the twin-immune zones), folded from // a merge of its own in P1 that only this flag makes. @@ -573,10 +571,15 @@ struct ProcessResult { // that was asked to validate looking like one that was not. Unset when no model was asked for. std::optional model_validation; - // What --mr found, when one was given; failure_reason says why where it could not run. + // How the --model was placed: RIGID_BODY (it fits against its null in the answer's merge or the twin + // hypothesis' one), MR (molecular replacement placed it) or NOT_PLACED. Empty without a model. + std::string model_placement; + // What molecular replacement found, where it ran; failure_reason says why where it could not. std::optional model_prior; - // The file prefix the model's choice was written under (_mr), empty where nothing was. + // The file prefix the model's choice of another group or the twin hypothesis' subgroup was written under + // (_mr), and that group; empty where nothing was. std::string model_prior_file_prefix; + std::string model_prior_file_space_group; // Space group used for the final re-scale/merge and written to the master file: determined by the // search when the user did not fix one, otherwise the fixed group. The whole group, so it keeps diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 6855c612f..83d61b1bf 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -81,9 +81,7 @@ std::string RugnuxCommandLine(const ProcessConfig &config, add("--reference-column", config.reference_column); } if (!config.model_path.empty()) - add("--model", config.model_path); - if (!config.mr_model_path.empty()) - add("--mr", config.mr_model_path); + add(config.force_mr ? "--mr" : "--model", config.model_path); } if (config.start_image != 0) add("-s", std::to_string(config.start_image)); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index eb55df05b..46747403c 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1425,8 +1425,19 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { std::string twin_overruled_reason; std::optional twin_choice; std::string twin_choice_reason; - // The twin hypothesis' merge made only for --mr to test: written, and reported, only where it chooses it. + // The twin hypothesis' merge made only for the model to test: written, and reported, only where it chooses it. bool twin_for_mr_only = false; + // The --model's kind, from the file alone: a prediction has no placement of its own and goes straight to + // molecular replacement; an experimental structure's group, where it has fewer rotations than the merge's, + // asks for the twin hypothesis' merge. + bool model_predicted = false; + size_t model_rotations = 0; + if (!config_.model_path.empty() && write_files && !geometry_prepass) { + const SearchModel m = PrepareSearchModel(config_.model_path); + model_predicted = m.failure_reason.empty() && m.predicted; + if (const gemmi::SpaceGroup *sg = m.st.find_spacegroup(); sg != nullptr && m.st.cell.is_crystal()) + model_rotations = sg->operations().sym_ops.size(); + } int twin_choice_group_number = 0; // the adopted group's, when the choice was made // A canonical pass the geometry walk is going to replace stops here, where the answer is // already in (GeometryWalkTakesRound); it is asked again below, before the analyses. Not with a @@ -2876,9 +2887,14 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // too, with the law - so a refinement program can take the twin hypothesis without a re-run. // The answer itself is not touched by any of this. const bool pre_twin_like = result.pre_promotion_twinning && result.pre_promotion_twinning->twinning_suspected; - // A predicted model can test the twin hypothesis itself (--mr's twin comparison), so with one the - // merge in the indexed lattice's group is made wherever a metric re-ask promoted from it. - const bool model_tests_twin = !config_.mr_model_path.empty() && twin_metric_subgroup.has_value(); + // A model can test the twin hypothesis itself (the twin validation, molecular replacement's twin + // comparison), so with one the merge in the subgroup is made wherever a metric re-ask promoted from the + // indexed lattice's group, or the model is in a group with fewer rotations than the adopted one. + const bool model_tests_twin = + !config_.model_path.empty() + && (twin_metric_subgroup.has_value() + || (!model_predicted && model_rotations > 0 + && model_rotations < result.space_group->operations().sym_ops.size())); if (!geometry_prepass && result.space_group && sg_search.point_group_order > 1 && (twin_overruled_subgroup || pre_twin_like || model_tests_twin)) { twin_choice = ChooseTwinSubgroup(*result.space_group, sg_search, @@ -2891,8 +2907,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { : pre_twin_like ? fmt::format("the merge read twin-like before the space-group search ({})", TwinningVerdictLine(*result.pre_promotion_twinning)) - : std::string("it was promoted from the indexed lattice's group, which the predicted model " - "can test"); + : std::string("the model can test it"); if (twin_choice) logger.Info("Twin hypothesis: {} was adopted, but {} - the merge in {} is written too, " "with twin law {}", result.space_group->xhm(), twin_choice_reason, @@ -3723,18 +3738,43 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { adopted.xhm(), c.subgroup.xhm(), t.twin_law, fraction.fraction, fraction.pairs, t.file_prefix); }; - if (result.consensus_cell && write_files && !config_.model_path.empty() && !geometry_prepass - && !superseded) { + // The model's choice of another group, or of the twin hypothesis' subgroup: never the answer's files, + // which take from a model at most the enantiomorph and the indexing, but files of their own beside + // them, _mr, from the merge it was chosen on, with the model as placed there. + struct ModelPriorFile { + std::vector merged; + UnitCell cell; + gemmi::SpaceGroup space_group; + std::optional law; + TwinLawFraction fraction; + std::shared_ptr placed; + }; + std::optional model_prior_file; + // The coordinates that go with the maps beside the answer's files, written once the model is done. + struct PlacedModelFile { + std::shared_ptr placed; + UnitCell cell; + gemmi::SpaceGroup space_group; + }; + std::optional placed_model_out; + const bool model_run = result.consensus_cell && write_files && !config_.model_path.empty() && !geometry_prepass + && !superseded; + // The merge's own shell bounds, so CC(model, data) is reported on the shells CC1/2 was. + std::vector report_shell_d_min; + for (const auto &sh : sm.statistics.shells) + report_shell_d_min.push_back(sh.d_min); + // Placed by the rigid body: the model fits against its null in the answer's merge or in the twin + // hypothesis' one. Where it is not - or is a prediction, or --mr asked for the search - it is searched + // for by molecular replacement below. + bool placed_by_rigid_body = false; + bool twin_chosen_by_model = false; + if (model_run && !model_predicted && !config_.force_mr) { phase("Validating against model"); Step("Model validation"); TimingMark mark(logger, "model validation"); const auto data_sg = experiment_.GetGemmiSpaceGroup(); // With a reference MTZ the merohedral indexing was already resolved against it (rotation // merge / stills scaling), so trust that; only probe indexing by R-free when model-only. - // The merge's own shell bounds, so CC(model, data) is reported on the shells CC1/2 was. - std::vector report_shell_d_min; - for (const auto &sh : sm.statistics.shells) - report_shell_d_min.push_back(sh.d_min); // The validation in the model's setting further down depends on this one only through the // setting and the indexing it settles, and both are forecast well before this one has finished // (ModelFrameForecast). So it is started on that forecast, beside this one - on the next card @@ -3784,6 +3824,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { } ModelValidationSchedule schedule; schedule.write_gate = gate; + schedule.fit_test = false; // the verdict is the first validation's (KeepModelVerdict) ahead.merged = std::move(merged); return ValidateAgainstModel(ahead.merged, cell, config_.model_path, config_.output_prefix, ahead.log, in, /*probe_indexing_ambiguity=*/false, @@ -3819,9 +3860,6 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // it. Without this the run ends successfully with no R-free, no maps and nothing in the // report - indistinguishable from a run that was never given --model at all. result.model_validation = validation; - if (!validation.failure_reason.empty()) - result.warnings.push_back({PathologyCode::MODEL_NOT_VALIDATED, - "Model validation did not run: " + validation.failure_reason}); // The ambiguity warning above was issued before the model was read; take it back only // where the model actually settled the indexing. if (validation.indexing_decided) @@ -3893,12 +3931,14 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { if (!(to_model == gemmi::Op::identity())) relabel_output(to_model.inverse(), "MODEL"); const auto sg_now = experiment_.GetGemmiSpaceGroup(); + ModelValidationSchedule kept_verdict; + kept_verdict.fit_test = false; remade = ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path, config_.output_prefix, logger, sg_now ? &*sg_now : nullptr, /*probe_indexing_ambiguity=*/false, static_cast(config_.nthreads), - experiment_.GetWavelength_A(), report_shell_d_min); + experiment_.GetWavelength_A(), report_shell_d_min, kept_verdict); } if (remade.ok) { KeepModelVerdict(remade, validation); @@ -3909,6 +3949,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { "first left its maps to it - no maps were written", remade.failure_reason); } } + placed_by_rigid_body = validation.ok && validation.model_fits; // The coordinates that go with the maps just written, in the frame the reflection files // beside them carry. Written whenever the validation ran, not only where the rigid-body // step moved the model: it is also re-fractionalized into the data cell and may be @@ -3920,7 +3961,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // subgroup's merge, plainly and as twinned by its law (the twin target, at a scanned fraction), // against the answer's merge as above. The maps and the placed model come from whichever merge // the model fits better on R-free; the answer's files are the answer's either way - the model - // chooses the maps, never the merge (see docs/RUGNUX_REPORT.md). + // chooses the maps, never the merge (see docs/RUGNUX_REPORT.md). A merge made only for the model to + // test is chosen only where the model fits against its null there too; otherwise it is not written. prepare_twin_file(); bool maps_from_twin = false; if (twin_file && validation.ok) { @@ -3928,17 +3970,20 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { const auto twin_r_free = [](const ModelValidationResult &v) { return v.twin_target ? std::min(v.r_free, v.r_free_twin) : v.r_free; }; - ModelValidationSchedule twin_schedule; - twin_schedule.maps_superseded = [&](const ModelValidationResult &v) { - return !(twin_r_free(v) < answer_r_free); + const auto maps_from = [&](const ModelValidationResult &v) { + return twin_r_free(v) < answer_r_free && (!twin_for_mr_only || v.model_fits); }; + ModelValidationSchedule twin_schedule; + twin_schedule.maps_superseded = [&](const ModelValidationResult &v) { return !maps_from(v); }; phase("Validating against model (twin hypothesis)"); auto twin_validation = ValidateAgainstModel( twin_file->merged, twin_file->cell, config_.model_path, config_.output_prefix, logger, &twin_file->space_group, /*probe_indexing_ambiguity=*/config_.reference_data.empty(), static_cast(config_.nthreads), experiment_.GetWavelength_A(), {}, twin_schedule, &twin_file->law); - maps_from_twin = twin_validation.ok && twin_r_free(twin_validation) < answer_r_free; + maps_from_twin = twin_validation.ok && maps_from(twin_validation); + const bool placed_in_twin = twin_validation.ok && twin_validation.model_fits; + placed_by_rigid_body = placed_by_rigid_body || placed_in_twin; // The subgroup file follows its own validation's indexing, as the answer follows its own: // the file, the R-factors and the maps then describe one indexing. if (twin_validation.ok) { @@ -3964,50 +4009,41 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { : std::string(), maps_from_twin ? "the merge in " + twin_file->space_group.xhm() : "the answer's merge"); if (maps_from_twin && twin_validation.placed_model) - placed_model_file = std::async(std::launch::async, - [placed = twin_validation.placed_model, cell = twin_file->cell, - sg = twin_file->space_group, prefix = config_.output_prefix, &logger] { - WritePlacedModel(*placed, cell, sg, prefix, logger); - }); + placed_model_out = PlacedModelFile{twin_validation.placed_model, twin_file->cell, + twin_file->space_group}; + // The model placed there and fitting better there chooses the subgroup: its files, _mr. + if (maps_from_twin && placed_in_twin) { + twin_chosen_by_model = true; + model_prior_file = ModelPriorFile{twin_file->merged, twin_file->cell, twin_file->space_group, + twin_file->law, twin_file->fraction, twin_validation.placed_model}; + result.model_prior_file_prefix = config_.output_prefix + "_mr"; + } result.twin_subgroup_merge->maps_from_this_merge = maps_from_twin; result.twin_subgroup_merge->model_validation = std::move(twin_validation); } if (validation.placed_model && !maps_from_twin) - placed_model_file = std::async(std::launch::async, - [placed = validation.placed_model, cell = *result.consensus_cell, - sg = experiment_.GetSpaceGroupOrP1(), prefix = config_.output_prefix, &logger] { - WritePlacedModel(*placed, cell, sg, prefix, logger); - }); + placed_model_out = PlacedModelFile{validation.placed_model, *result.consensus_cell, + experiment_.GetSpaceGroupOrP1()}; } prepare_twin_file(); - // --mr (experimental): the model placed by molecular replacement, as many copies as the cell holds, in the - // groups the answer's merge cannot tell apart and in the twin hypothesis' subgroup. It takes the - // enantiomorph label into the answer, as --model does, and nothing else: a different group it - // chooses goes to files of its own beside the answer's, _mr, from the merge it was chosen on. - struct ModelPriorFile { - std::vector merged; - UnitCell cell; - gemmi::SpaceGroup space_group; - std::optional law; - TwinLawFraction fraction; - std::shared_ptr placed; - }; - std::optional model_prior_file; - if (result.consensus_cell && write_files && !config_.mr_model_path.empty() && !geometry_prepass - && !superseded) { - phase("Placing the predicted model"); + // Molecular replacement: for a prediction, which has no placement of its own; with --mr; and for an + // experimental structure the rigid body did not place (or could not use: no crystal cell). The model is + // placed, as many copies as the cell holds, in the groups the answer's merge cannot tell apart and in the + // twin hypothesis' subgroup. It takes the enantiomorph label into the answer - only past the anomalous + // veto - and nothing else, never an indexing: a different group or the subgroup it chooses goes to the + // files of its own, _mr, from the merge it was chosen on. It never makes the rigid body's + // verdict (MODEL_FIT) either; its own null is what placed it. + if (model_run && (model_predicted || config_.force_mr || !placed_by_rigid_body)) { + phase("Molecular replacement"); Step("Molecular replacement"); TimingMark mark(logger, "model prior"); std::vector merges{{"MERGE", &sm.merged, *result.consensus_cell, experiment_.GetSpaceGroupOrP1()}}; if (twin_file) merges.push_back({"TWIN_SUBGROUP_MERGE", &twin_file->merged, twin_file->cell, twin_file->space_group, twin_file->law}); - MolecularReplacementResult mr = RunMolecularReplacement(config_.mr_model_path, merges, logger, + MolecularReplacementResult mr = RunMolecularReplacement(config_.model_path, merges, logger, static_cast(config_.nthreads)); - if (!mr.ok) - result.warnings.push_back({PathologyCode::MODEL_NOT_VALIDATED, - "The --mr model could not be searched for: " + mr.failure_reason}); if (mr.chosen >= 0) { const MRCandidate &c = mr.candidates[mr.chosen]; if (mr.decision == "ENANTIOMORPH") @@ -4036,15 +4072,64 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { } else if (mr.decision == "TWIN_SUBGROUP" && twin_file) { model_prior_file = ModelPriorFile{twin_file->merged, twin_file->cell, c.space_group, twin_file->law, twin_file->fraction, c.placed}; + twin_chosen_by_model = true; + } + result.model_prior_file_prefix = model_prior_file ? config_.output_prefix + "_mr" : std::string(); + // The maps, the R-factors and CC(model, data): the placement validated against the merge it was + // chosen on, with nothing left for that validation to decide. The rigid body's verdict stays. + const ModelPriorFile *own = model_prior_file ? &*model_prior_file : nullptr; + const gemmi::SpaceGroup data_sg = own ? own->space_group : experiment_.GetSpaceGroupOrP1(); + const UnitCell cell = own ? own->cell : *result.consensus_cell; + ModelValidationSchedule placed_elsewhere; + placed_elsewhere.fit_test = false; + phase("Validating against model (molecular replacement)"); + ModelValidationResult v = ValidateAgainstModel( + own ? own->merged : sm.merged, cell, config_.model_path, + own ? result.model_prior_file_prefix : config_.output_prefix, logger, &data_sg, + /*probe_indexing_ambiguity=*/false, static_cast(config_.nthreads), + experiment_.GetWavelength_A(), report_shell_d_min, placed_elsewhere, nullptr, c.placed.get()); + if (v.ok) { + if (result.model_validation) { + const ModelValidationResult &rb = *result.model_validation; + v.fit_tested = rb.fit_tested; + v.model_fits = rb.model_fits; + v.null_r_work_mean = rb.null_r_work_mean; + v.null_r_work_sd = rb.null_r_work_sd; + v.r_work_sigma = rb.r_work_sigma; + v.null_replicates = rb.null_replicates; + } + if (own) + model_prior_file->placed = v.placed_model; + else + placed_model_out = PlacedModelFile{v.placed_model, cell, data_sg}; + result.model_validation = std::move(v); } } - if (model_prior_file) - result.model_prior_file_prefix = config_.output_prefix + "_mr"; - if (twin_for_mr_only && mr.decision != "TWIN_SUBGROUP") { + result.model_prior = std::move(mr); + } + if (model_run) { + result.model_placement = placed_by_rigid_body ? "RIGID_BODY" + : result.model_prior && result.model_prior->chosen >= 0 ? "MR" : "NOT_PLACED"; + // A model that was asked for and could not be used has to say so where anyone will see it. Without + // this the run ends successfully with no R-free, no maps and nothing in the report - indistinguishable + // from a run that was never given --model at all. A model that was used and not placed is a result. + std::string why = result.model_validation ? result.model_validation->failure_reason : std::string(); + if (why.empty() && result.model_prior && !result.model_prior->ok) + why = "the model could not be searched for: " + result.model_prior->failure_reason; + if (!(result.model_validation && result.model_validation->ok) && !why.empty()) + result.warnings.push_back({PathologyCode::MODEL_NOT_VALIDATED, "Model validation did not run: " + why}); + // The twin hypothesis' merge made only for the model to test: written, and reported, only where it chose it. + if (twin_for_mr_only && !twin_chosen_by_model) { twin_file.reset(); result.twin_subgroup_merge.reset(); } - result.model_prior = std::move(mr); + if (model_prior_file) + result.model_prior_file_space_group = model_prior_file->space_group.xhm(); + if (placed_model_out) + placed_model_file = std::async(std::launch::async, [f = *placed_model_out, prefix = config_.output_prefix, + &logger] { + WritePlacedModel(*f.placed, f.cell, f.space_group, prefix, logger); + }); } // The unmerged MTZ (written with the other output files, see there) is built from here on, diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 285e561a3..ab0f8e056 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -201,8 +201,8 @@ void print_usage() { std::cout << " --scaling-iterations Cap on the per-frame scaling iterations; the loop stops earlier once the scales settle (default: 100)" << std::endl; std::cout << " -z, --reference Reference data of the same crystal form: an MTZ, or a PDB structure-factor mmCIF (e.g. 1abc-sf.cif, first merged block with a usable column; R-free set from _refln.status f/o, else pdbx_r_free_flag). Gzipped or not, recognised by content. --reference-mtz is the same option" << std::endl; std::cout << " --reference-column