diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 64064a11d..91ca925f4 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -172,8 +172,8 @@ axis preserves the class of a half-integer pseudo-translation, which is what a p That robustness does not extend to a pseudo-translation which is not half-integer, and the partner steps are restricted when one is detected - see the translational-pseudo-symmetry note below. -**The twin fraction of the twin hypothesis** follows Yeates: the intensity a twinned crystal records -is (1−α) I(h) + α I(Th), and the disagreement +**The twin fraction of the twin hypothesis, and the twinned model it is scored with**, follow Yeates: +the intensity a twinned crystal records is (1−α) I(h) + α I(Th), and the disagreement H = |I(h) − I(Th)|/(I(h) + I(Th)) of the pairs a law relates is uniform on [0, 1−2α], so the median gives the fraction. T. O. Yeates, "Detecting and overcoming crystal twinning" (1997), Methods Enzymol. 276, 344-358 [doi:10.1016/S0076-6879(97)76068-3](https://doi.org/10.1016/S0076-6879(97)76068-3). diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index b592e9dc5..937391670 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -47,7 +47,7 @@ The methods draw on, and in places reimplement, solutions from: - P. Gasparotto, L. Barba, H.-C. Stadler et al., "TORO Indexer: a PyTorch-based indexing algorithm for kilohertz serial crystallography", *J. Appl. Cryst.* **57** (2024), 931-944 (the algorithm of the `ffbidx` fast-feedback indexer, §4). - I. Křivý & B. Gruber, "A unified algorithm for determining the reduced (Niggli) cell", *Acta Cryst.* **A32** (1976), 297-298, and International Tables for Crystallography Vol. A, Table 9.2.5.1 (the Niggli reduction and the lattice-character table of §5.3/§6). - J. E. Padilla & T. O. Yeates, "A statistic for local intensity differences: robustness to anisotropy and pseudo-centering and utility for detecting twinning", *Acta Cryst.* **D59** (2003), 1124-1130 (the L test, §13.2). -- T. O. Yeates, "Detecting and overcoming crystal twinning", *Methods Enzymol.* **276** (1997), 344-358 (the $H$ statistic behind the twin fraction of the twin hypothesis, §13.2). +- T. O. Yeates, "Detecting and overcoming crystal twinning", *Methods Enzymol.* **276** (1997), 344-358 (the twinned intensity model and the $H$ statistic behind the twin fraction of the twin hypothesis, §13.2, §14.9). - A. J. C. Wilson, "The probability distribution of X-ray intensities", *Acta Cryst.* **2** (1949), 318-321, and E. R. Howells, D. C. Phillips & D. Rogers, "The probability distribution of X-ray intensities. II. Experimental investigation and the X-ray detection of centres of symmetry", *Acta Cryst.* **3** (1950), 210-214 (the acentric and centric intensity distributions behind the centre-of-symmetry statistics of §13.1 and the Wilson outlier test of §13.3). - W. H. Baur & D. Kassner, "The perils of Cc: comparing the frequencies of falsely assigned space groups with their general population", *Acta Cryst.* **B48** (1992), 356-369 (writing the centrosymmetric group where the absences cannot decide the centre, §13.1). - R. J. Read, P. D. Adams & A. J. McCoy, "Intensity statistics in the presence of translational noncrystallographic symmetry", *Acta Cryst.* **D69** (2013), 176-183 (the native-Patterson detection of translational pseudo-symmetry, and the intensity modulation it produces, which the axial-zone screw-absence test scores against). diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index c6812bf84..ae3571ea6 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -64,7 +64,7 @@ angle < 0.44$) and the averaging closed a third or more of the gap from the unme A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$ — 0.500 and 0.333 untwinned, 0.375 and 0.200 for a perfect twin) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ (2.0 for untwinned acentric data, 1.5 for a perfect twin) are written to the merged mmCIF as a twinning diagnostic. Both are taken on intensities divided by their resolution-shell mean, over the shells whose $\langle I/\sigma\rangle$ reaches 1, with Wilson outliers rejected; the $L$-test pairs are therefore compared on one scale even where two index steps span a steep fall-off, as in phenix.xtriage and ctruncate, and the selection is by shell, never by the individual reflection's $I/\sigma$, which would cut the weak tail and bias $\langle|L|\rangle$ down. The twin fraction is quoted from the statistic that carries the verdict — the $L$-test unless the call rests on the second moment alone — and a second moment is not turned into a fraction under a detected pseudo-translation, which inflates it. A twin law exists only where the **lattice** has a rotation the group lacks, so whether one is possible is read off the cell, not off the name of the Laue class: the lattice's two-folds are found as for the free set (Le Page, within 3° obliquity — GEMMI's `find_twin_laws`), and where none lies outside the group no twin is called. That covers the holohedral classes ($4/mmm$ on a tetragonal lattice, $6/mmm$, $m\bar{3}m$, $\bar{3}m$ on a rhombohedral one) and equally $mmm$ or $2/m$ on a lattice with no special metric, while a pseudo-merohedral $2/m$ or $mmm$ lattice — a monoclinic $\beta$ near 90°, a C-centred cell with $b \approx a\sqrt{3}$ — admits a law and stays eligible. The $L$-test is still read where no law is possible, for a different question: merging $I(h)$ with $I(Th)$ under a false operator $T$ gives $(I(h)+I(Th))/2$ whatever the twin fraction, which has the perfect-twin distribution, so $\langle|L|\rangle$ below 0.42 there is reported as **an adopted operator averaging unequal intensities** — the space group is too high, or a twin law was absorbed into the point group — a warning, never a change to the space group. A genuine operator leaves the untwinned 0.5. Reflections that overlap along a very long axis narrow the distribution the same way, which is one reason it stays a warning. The same numbers measured on the P1 merge the space-group search was given, before any point group was adopted, are reported beside them (`_BEFORE_SEARCH`), with that merge's own pseudo-translation declared and the reflections its lattice centring extinguishes left out. Between that bar and the twin band's 0.44 the verdict is `NO_TWIN_LAW`: the intensities read twin-like, but the lattice admits no law for the group, so they are not a twin in it — the group may be too high (a subgroup's twin law absorbed, or non-crystallographic symmetry taken for real), or something other than twinning narrowed the distribution. It is reported as such rather than as a twin that "needs a twin law" no refinement program could be given. -**The twin hypothesis.** Where the adopted group may have absorbed a twin law — a twin gate of the search refused it and was overruled (§13.1), or a group above $P1$ was adopted on a merge that read twin-like before the search — the same observations are also merged in a subgroup and written beside the answer as `_twin_.mtz` / `.cif`, in the subgroup's standard setting, to the answer's resolution and with the answer's free set (each reflection takes the flag of the reflection of the answer's merge its orbit falls in, so twin mates share it). The subgroup is the one a decision named — the group a refused promotion was adopted over, or the indexed lattice's group a metric re-ask promoted from — and otherwise, among the index-2 subgroups of the adopted group that the search confirmed as point-group finalists, the one whose lost two-fold agrees worst with the intensities (lowest operator CC): the most likely twin law. One file is written; the other subgroups that qualified are named in the report. The law is in the mmCIF's `_pdbx_reflns_twin` with its fraction, from Yeates' $H$ on the pairs it relates in the subgroup merge — $H = |I(h)-I(Th)|/(I(h)+I(Th))$ is uniform on $[0, 1-2\alpha]$, so $\alpha = 1/2 - \mathrm{median}\,H$, a lower bound since measurement error only adds to $H$. Nothing the answer is made of depends on it: it is a second merge from the same scaled observations, and the answer's files are byte for byte what they are without it. Refinement programs want twinned data merged in the true lower group, never detwinned, which is why the file is the merge and not a correction of the answer; a near-perfect twin of a subgroup and the higher group predict the same intensities, so which of the two is the crystal is a question for a model. +**The twin hypothesis.** Where the adopted group may have absorbed a twin law — a twin gate of the search refused it and was overruled (§13.1), or a group above $P1$ was adopted on a merge that read twin-like before the search — the same observations are also merged in a subgroup and written beside the answer as `_twin_.mtz` / `.cif`, in the subgroup's standard setting, to the answer's resolution and with the answer's free set (each reflection takes the flag of the reflection of the answer's merge its orbit falls in, so twin mates share it). The subgroup is the one a decision named — the group a refused promotion was adopted over, or the indexed lattice's group a metric re-ask promoted from — and otherwise, among the index-2 subgroups of the adopted group that the search confirmed as point-group finalists, the one whose lost two-fold agrees worst with the intensities (lowest operator CC): the most likely twin law. One file is written; the other subgroups that qualified are named in the report. The law is in the mmCIF's `_pdbx_reflns_twin` with its fraction, from Yeates' $H$ on the pairs it relates in the subgroup merge — $H = |I(h)-I(Th)|/(I(h)+I(Th))$ is uniform on $[0, 1-2\alpha]$, so $\alpha = 1/2 - \mathrm{median}\,H$, a lower bound since measurement error only adds to $H$. Nothing the answer is made of depends on it: it is a second merge from the same scaled observations, and the answer's files are byte for byte what they are without it. Refinement programs want twinned data merged in the true lower group, never detwinned, which is why the file is the merge and not a correction of the answer; a near-perfect twin of a subgroup and the higher group predict the same intensities, so which of the two is the crystal is a question for a model (§14.9). Beside them, on rotation data with `--finalist-ledger`, the run reports **twin-immune zone evidence** for the operators the adopted point group adds over each of its index-2 subgroups, read on a full-resolution merge of the observations in P1 that the option makes for the purpose. Under a twin law $T$, $I_\mathrm{obs}(h)=(1-\alpha)I(h)+\alpha I(Th)$; a reflection whose twin mate is itself up to the subgroup and Friedel is untouched at every $\alpha$, and those are exactly the reflections centric in the group but acentric in the subgroup. They read centric ($\langle|E^2-1|\rangle = 0.968$) if the added operators are real and acentric (0.736) if they are a twin law or a pseudo-symmetry — the one intensity statistic that still separates the two at $\alpha = 0.5$, where every operator statistic reads "real". Each zone is normalised against its own mean in resolution bins (and within the two phase classes of a detected pseudo-translation), read only in shells with $\langle I/\sigma\rangle \ge 5$ because noise inflates every class towards centric, and reported with $n$, a standard error and the centric-over-acentric Wilson log-likelihood ratio in nats (both densities convolved with each reflection's measurement error), beside the acentric control. It is read absolutely, never as the difference to the control: a perfect twin's control (0.541) plus noise inflation would otherwise look like true symmetry. An absolute reading is only as good as the normalisation, though, and the acentric control certifies it: an acentric population reads $-0.130$ nats per reflection when the normalisation is right, and whatever the control reads above that is the normalisation's — anisotropy, a pseudo-translation, a pseudo-centring, noise all inflate every class towards centric alike — so the zone, normalised the same way, carries the same per reflection and the **calibrated** evidence has it taken off (a twinned control reads *below* the expectation, so nothing is taken off a twin). The centric side is read **twinned at the fraction the lattice's other operations show** by their own correlation — the strongest CC of a lattice rotation outside the group, relative to the mean CC of the group's own operators, inverted through $\rho = 2\alpha(1-\alpha)/((1-\alpha)^2+\alpha^2)$ — because a twin by another law reaches a genuine zone exactly as it reaches the control: a genuine 321 crystal twinned by a 622 law read its 2-folds at −206 nats against the untwinned centric density and +209 at that law's fraction of 0.07. The operators' own twin law cannot reach their zone, so the acentric side stays untwinned, and the control's expectation is taken at the same fraction. Measured, a 6/m crystal with a 67 Ų anisotropy read its control at $+0.10$ nats per reflection and the zone of a refused 622 the same, so the zone's $+124$ nats were the normalisation's; calibrated it reads $-115$, and genuine promotions keep $+200$ and above. The calibrated zone is what a refused promotion is decided on (above). It is an in-house method: published practice (Yeates' $H$-test, xtriage) excludes these reflections as uninformative about the twin fraction, and no published test was found that uses them positively; the one precedent here is a trigonal crystal whose twin-immune zone read centric and whose higher group was then confirmed by refinement. @@ -273,3 +273,9 @@ The rotation and translation actually taken are reported (`RIGID_BODY_ROTATION_D The coordinates **as placed** are written as `_model.cif`, and as `_model.pdb` where the PDB format can hold the cell (the fragment-screening tools the file feeds — PanDDA, and the pair dimple produces — take a PDB beside the MTZ), so there is a coordinate file that describes the maps: the input's chains, residues, ligands, waters, B-factors, occupancies and anisotropic $U$s, moved, in the same cell and space group as `.mtz`. Taking the frame from the written reflections rather than from the input model matters — §14.6 may have relabelled them to the model's enantiomorph, which is neither the data's original group nor, necessarily, the model's. It is written **whenever the maps are**, not only where the rigid-body step was committed: the model is re-fractionalized into the data cell and may be relabelled whatever the placement decided, so an unmoved model is still not the input file; and a model the null of §14.5 rejected is scored, placed and mapped like any other, which is precisely the case where the density is worth looking at. +### 14.9 The twin hypothesis: which merge the maps come from + +Where the run wrote a twin hypothesis (§13.2), the model is scored against that merge too: as merged, and **as twinned** by its law — $|F_\mathrm{twin}(h)|^2 = (1-\alpha)|F_\mathrm{model}(h)|^2 + \alpha|F_\mathrm{model}(Th)|^2$ from the placed, scaled model, one overall scale refitted at each $\alpha$, $\alpha$ scanned from 0 to 0.5 in steps of 0.01 on the working set and R-free read at the chosen one. The placement, the anisotropic $B$ and the bulk solvent are those of the untwinned fit (§14.2, §14.8), so this is the same few-parameter comparison as §14.3, not a twin refinement. Against the answer's merge the model is scored as in §14.3, the subgroup model's structure factors simply read at each of the higher group's reflections. + +The maps and the placed model are then made from whichever merge the model fits better on R-free — the lower of the subgroup's plain and twinned R-free against the answer's — and the report says which (`MAPS_FROM=`) and with what R-free for each. The map coefficients are those of §14.4 against the chosen merge's amplitudes, not detwinned. **This chooses the maps and nothing else**: the space group, the answer's reflection files and every number in the report outside the twin keys are what the run decided without the model, because a model can be wrong and the reflection files must not depend on one (§14.5). The structure factors come from the same GPU engine as §14.1 where there is a GPU. + diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 72a683ba6..0ac0a01f1 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -170,7 +170,15 @@ PDB: PanDDA's per-dataset input is `.pdb` beside `.mtz`, and dimple pair. (The PDB is skipped, with a log line, for a cell its fixed-width format cannot hold; the mmCIF is unconditional.) That is the file to open with the maps, and the one to hand to REFMAC5 or `phenix.refine`; its starting R-free is the `R_FREE=` the report quotes. -(`.cif` is the merged **reflections** — hence the separate name.) Both are written whenever the +(`.cif` is the merged **reflections** — hence the separate name.) + +Where the run also wrote a **twin hypothesis** — `_twin_.mtz`/`.cif`, the same data merged +in a subgroup with a twin law, written when the adopted space group may have absorbed one (see +[the twin hypothesis](RUGNUX_REPORT.md#the-twin-hypothesis-twin_subgroup_)) — the model is scored +against that merge as well, plainly and as twinned by the law, and the maps and `_model.cif` +come from whichever of the two merges the model fits better on R-free. The report says which +(`MAPS_FROM=`) and gives the R-free of both. The model only chooses the maps: `.mtz` and the +space group are the run's own either way. Both are written whenever the maps are, including for a model the data rejected: the rejection is a result, and it is exactly the case where someone wants to look at the model in the density. diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 8634ca78b..fd76f8b88 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -861,6 +861,14 @@ as usual. | `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 | +Where a twin hypothesis was written (above), the model is scored against that merge too - as merged, +and as twinned by its law at a fraction scanned on the working set - and **the maps and the placed model +come from whichever merge the model fits better on R-free**. `MAPS_FROM=` says which (`MERGE` or +`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 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 verdict (`MODEL_FIT`, `MODEL_FIT_SIGMA`, `MODEL_DECISIONS_TAKEN` and the two decisions). diff --git a/docs/RUGNUX_TUTORIAL.md b/docs/RUGNUX_TUTORIAL.md index 89eac3cfe..67988771c 100644 --- a/docs/RUGNUX_TUTORIAL.md +++ b/docs/RUGNUX_TUTORIAL.md @@ -287,7 +287,8 @@ Absorption is corrected only by the empirical surface a rotation run fits and ke [Validating against a model](RUGNUX_ADVANCED.md#validating-against-a-model-rugnux---model). - `_twin_.mtz` / `.cif` — only where the chosen space group may have absorbed a twin law: the same data merged in a subgroup ``, with the twin law and its fraction, same resolution and - free set as `.mtz`, for refinement with that law. See + free set as `.mtz`, for refinement with that law. With `--model`, the maps come from + whichever of the two merges the model fits better. See [The twin hypothesis](RUGNUX_REPORT.md#the-twin-hypothesis-twin_subgroup_). - `_report.txt` — the **results report**: what the run determined, in a form both a person and a beamline script can read. Always written, next to the files above. See diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index e355be145..ab8797727 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -423,6 +423,7 @@ ModelValidationResult Validate(const std::vector &merged, double wavelength_A, const std::vector &report_shell_d_min, const ModelValidationSchedule &schedule, + const gemmi::Op *twin_law, bool rigid_body_gpu) { ModelValidationResult result; result.model_path = model_path; @@ -1228,6 +1229,61 @@ ModelValidationResult Validate(const std::vector &merged, result.b_sol = best.b_sol; result.k_overall = best.k_overall; + // --- the twin target, where a twin law was given --- + // The intensity a twinned crystal records at h is (1-a) I(h) + a I(Th), so the model is compared + // against the data as twinned: |F_twin(h)|^2 = (1-a)|F_model(h)|^2 + a|F_model(Th)|^2 from the + // placed, scaled model, with one overall scale refitted at each a. The fraction is scanned in steps + // of 0.01 and chosen on the working set; R-free is read at it. Nothing else is refitted against the + // twinned target - the placement, the anisotropic B and the bulk solvent are those of the untwinned + // fit - so this is the same few-parameter comparison as R above, not a twin refinement. + // Following Yeates (1997) Methods Enzymol. 276, 344-358 + if (twin_law != nullptr) { + std::unordered_map model_i; + model_i.reserve(fmodel.v.size()); + for (const auto &hv : fmodel.v) + model_i[hkl_key(hv.hkl)] = std::norm(hv.value); + struct TwinTerm { double fo, i1, i2; bool free; }; + std::vector twin_terms; + for (const auto &hv : fmodel.v) { + const auto it = obs_by_hkl.find(hkl_key(hv.hkl)); + if (it == obs_by_hkl.end()) continue; + const gemmi::Miller mate = asu.to_asu(twin_law->apply_to_hkl(hv.hkl), gops).first; + const auto m = model_i.find(hkl_key(mate)); + if (m == model_i.end()) continue; + twin_terms.push_back({it->second.F, std::norm(hv.value), m->second, it->second.free}); + } + double best_r_work = 2.0; + for (int step = 0; step <= 50 && !twin_terms.empty(); ++step) { + const double a = 0.01 * step; + double fo_fc = 0, fc_fc = 0; + for (const auto &t : twin_terms) + if (!t.free) { + const double fc = std::sqrt((1.0 - a) * t.i1 + a * t.i2); + fo_fc += t.fo * fc; + fc_fc += fc * fc; + } + const double k = fc_fc > 0 ? fo_fc / fc_fc : 1.0; + double num_w = 0, den_w = 0, num_f = 0, den_f = 0; + for (const auto &t : twin_terms) { + const double diff = std::fabs(t.fo - k * std::sqrt((1.0 - a) * t.i1 + a * t.i2)); + if (t.free) { num_f += diff; den_f += t.fo; } + else { num_w += diff; den_w += t.fo; } + } + if (den_w > 0 && den_f > 0 && num_w / den_w < best_r_work) { + best_r_work = num_w / den_w; + result.r_work_twin = best_r_work; + result.r_free_twin = num_f / den_f; + result.twin_fraction = a; + } + } + result.twin_target = std::isfinite(result.r_free_twin); + if (result.twin_target) + logger.Info("Model validation: against the data as twinned by {} - fraction {:.2f} (scanned on the " + "working set) - R-work {:.4f}, R-free {:.4f}; untwinned {:.4f} / {:.4f}", + twin_law->as_hkl().triplet('h'), result.twin_fraction, result.r_work_twin, + result.r_free_twin, result.r_work, result.r_free); + } + // --- CC(model, data) by resolution shell --- // The correlation of the merged intensities with |F_model|^2, binned on the merge table's own // shells so the two tables line up row for row. Nearly free: the scaled Fmodel is already here, @@ -1655,14 +1711,16 @@ ModelValidationResult ValidateAgainstModel(const std::vector & size_t nthreads, double wavelength_A, const std::vector &report_shell_d_min, - const ModelValidationSchedule &schedule) { + const ModelValidationSchedule &schedule, + const gemmi::Op *twin_law) { #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, // so the reflection files are those of a run without a model. try { return Validate(merged, cell, model_path, output_prefix, logger, data_space_group, - probe_indexing_ambiguity, nthreads, wavelength_A, report_shell_d_min, schedule, true); + probe_indexing_ambiguity, nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, + true); } catch (const RigidBodyGPUFailure &e) { cuda_clear_error(); ModelValidationResult failed; @@ -1673,7 +1731,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, false); + nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, false); #endif } diff --git a/image_analysis/structure_refinement/ModelValidation.h b/image_analysis/structure_refinement/ModelValidation.h index c89feb567..adc686bf2 100644 --- a/image_analysis/structure_refinement/ModelValidation.h +++ b/image_analysis/structure_refinement/ModelValidation.h @@ -184,6 +184,14 @@ struct ModelValidationResult { bool adopted_model_enantiomorph = false; int model_space_group_number = 0; + // Against the data as twinned by the twin law ValidateAgainstModel was given (none given: not + // made): R-work and R-free with |F_twin|^2 = (1-a)|F_model(h)|^2 + a|F_model(Th)|^2 and one overall + // scale, at the fraction a (0 to 0.5) that gives the lowest R-work. + bool twin_target = false; + double twin_fraction = NAN; + double r_work_twin = NAN; + double r_free_twin = NAN; + // The model as it was scored: re-fractionalized into the data cell and moved by the rigid-body // step, everything else exactly as it was read. Null where the validation did not run. The frame // it is written in is only settled by AdoptModelFrame below, so it is the caller that puts it on @@ -224,6 +232,9 @@ struct ModelValidationResult { // nthreads is what the null's replicates and the candidate changes of basis run on - each is // independent of the others, so they run at once. The answer does not depend on it. // +// 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). +// // 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. @@ -261,7 +272,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector & size_t nthreads = 1, double wavelength_A = 0.0, const std::vector &report_shell_d_min = {}, - const ModelValidationSchedule &schedule = {}); + const ModelValidationSchedule &schedule = {}, + const gemmi::Op *twin_law = 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 497d3150b..d8f50ed55 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -2392,6 +2392,33 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, } if (!mv.maps_prefix.empty()) Add(s, KeyText("MAPS_PREFIX", mv.maps_prefix)); + // 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) { + const auto &t = *result.twin_subgroup_merge; + const auto &tv = *t.model_validation; + Add(s, KeyEnum("MAPS_FROM", t.maps_from_this_merge ? "TWIN_SUBGROUP_MERGE" : "MERGE", + {"MERGE", "TWIN_SUBGROUP_MERGE"})); + if (tv.ok) { + Add(s, KeyReal("TWIN_SUBGROUP_R_WORK", tv.r_work, "{:.4f}")); + Add(s, KeyReal("TWIN_SUBGROUP_R_FREE", tv.r_free, "{:.4f}")); + if (tv.twin_target) { + Add(s, KeyReal("TWIN_SUBGROUP_R_WORK_TWINNED", tv.r_work_twin, "{:.4f}")); + Add(s, KeyReal("TWIN_SUBGROUP_R_FREE_TWINNED", tv.r_free_twin, "{:.4f}")); + Add(s, KeyReal("TWIN_SUBGROUP_MODEL_TWIN_FRACTION", tv.twin_fraction, "{:.2f}")); + } + } + Add(s, Blank()); + Add(s, Prose(fmt::format( + " The model was also scored against the twin hypothesis' merge in {} ({}.mtz): as merged,\n" + " and as twinned by {} - |F_twin|^2 = (1-a)|F_model(h)|^2 + a|F_model(Th)|^2, the\n" + " fraction a scanned on the working set. The maps and the placed model come from {},\n" + " the merge the model fits better on R-free ({:.4f} against {:.4f}). This chooses the\n" + " maps only: the answer's reflection files and every number above are the answer's.", + t.space_group, t.file_prefix, t.twin_law, + t.maps_from_this_merge ? "the merge in " + t.space_group : "the answer's merge", + tv.twin_target ? std::min(tv.r_free, tv.r_free_twin) : tv.r_free, mv.r_free))); + } if (!mv.cc_model_shells.empty()) { Add(s, Blank()); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index c6eb42c66..782347db0 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -491,6 +491,10 @@ struct ProcessResult { double twin_fraction = NAN; // EstimateTwinLawFraction: a lower bound int twin_fraction_pairs = 0; std::vector other_subgroups; // other confirmed subgroups that qualified + // --model: the model scored against this merge (plain and as twinned by the law), and whether + // the maps and the placed model come from it rather than from the answer's merge. + std::optional model_validation; + bool maps_from_this_merge = false; }; std::optional twin_subgroup_merge; diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 12641df51..4dde9ed34 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -3738,7 +3738,61 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // file". A rejected model is scored, placed and mapped like any other, and that is the case // where someone most wants to see the model in the density, so it gets its file too. // Written beside the reflection files below: nothing after this reads it. - if (validation.placed_model) + // The model asked about the twin hypothesis too, where one was written: scored against the + // 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). + prepare_twin_file(); + bool maps_from_twin = false; + if (twin_file && validation.ok) { + const double answer_r_free = validation.r_free; + 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); + }; + 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; + // 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) { + const gemmi::Op op = twin_validation.indexing_op; + twin_file->space_group = *AdoptModelFrame(twin_validation, twin_file->merged, + twin_file->space_group, + experiment_.GetScalingSettings().GetMergeFriedel(), + logger); + if (!(op == gemmi::Op::identity())) + twin_file->law = op.inverse().combine(twin_file->law).combine(op); + result.twin_subgroup_merge->twin_law = twin_file->law.as_hkl().triplet(); + result.twin_subgroup_merge->space_group = twin_file->space_group.xhm(); + } + logger.Info("Model validation, twin hypothesis: R-free {:.4f} against the answer's merge in {}, " + "{:.4f} against the merge in {}{} - the maps and the placed model come from {}", + answer_r_free, experiment_.GetSpaceGroupOrP1().xhm(), twin_validation.r_free, + twin_file->space_group.xhm(), + twin_validation.twin_target + ? fmt::format(" and {:.4f} as twinned at a fraction of {:.2f}", + twin_validation.r_free_twin, twin_validation.twin_fraction) + : 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); + }); + 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] {