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

rugnux: measure and report diffraction anisotropy

Browse Source 
rugnux now says whether a dataset's fall-off is direction-dependent, and by
how much. It corrects nothing and truncates nothing: no intensity is changed,
no reflection is dropped on a directional criterion, and the written files do
not depend on direction at all.

Two quantities, because they are not the same thing. The anisotropic deltaB is
the range of the principal components of the anisotropy tensor - a rate of
fall-off. The diffraction limit along each principal direction is where
<I/sigma(I)> in a 20 degree cone falls through 2 - where signal actually runs
out. One battery case has only 0.28 A between its directional limits and a 58x
ratio in cone <I/sigma>, so reporting either alone would miss it.

The tensor is a Laue-constrained deviatoric ADP tensor fitted on INTENSITIES
with no positivity cut, by weighted Gauss-Newton over 12 shells x 60
directions with a free constant per shell. Fitting amplitudes after a
positivity cut, which is what xtriage and ctruncate do, destroys about 40% of
the measured anisotropy - the cut keeps only the positive noise excursions in
whichever direction has died, and that is the direction carrying the signal.
Against the same 38 merged files rugnux reads 1.24x xtriage's eigenvalue
spread and 1.61x ctruncate's; on strong near-isotropic data all three agree to
a few percent, and they diverge exactly where a direction has died.

The verdict is gated three ways - not detected, detected, or cannot determine
- against the dataset's own systematic floor, measured in the tensor
directions its Laue symmetry forbids. The floor cannot be measured on merged
reflections, which have exact Laue symmetry by construction, so the floor is
taken from the unmerged observations and the verdict is "cannot determine"
without them. Triclinic has no forbidden subspace and always returns cannot
determine. A cubic crystal returns exactly zero, because that is its symmetry
and not a measurement.

A second axis reports the resolution signature: a genuine Debye-Waller
fall-off is linear through the origin in s^2, and a deficit that is flat is
something else. Magnitude alone had promoted a crystal that is 68% not a
Debye-Waller B into the top five of this battery; it now reads not detected
with the caution attached.

Following Sheriff & Hendrickson (1987) Acta Cryst. A43, 118-121 for the tensor
and Popov & Bourenkov (2003) Acta Cryst. D59, 1145-1153 for the estimator.

The directional limits are written as jfjoch_ local mmCIF items rather than
_reflns.pdbx_aniso_diffraction_limit_*, whose dictionary definition is
explicitly the ellipsoid fitted to a diffraction cut-off surface - a
construction rugnux does not perform. The generic anisotropic B tensor items
are written.

Changes no existing number; only REPORT_VERSION moves, 1 to 2.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P
This commit is contained in:
leonarski_f
2026-08-25 19:13:16 +02:00
co-authored by Claude Opus 5
parent 0f5a6e41d6
commit 80e10f51b3
13 changed files with 1842 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/ACKNOWLEDGEMENT.md
+15
View File
@@ -101,6 +101,21 @@ enumerated from a cell. D. W. Moreau, H. Atakisi and R. E. Thorne, "Ice in biomo
cryocrystallography" (2021), Acta Cryst. D77, 540-554
[doi:10.1107/S2059798321001170](https://doi.org/10.1107/S2059798321001170).
**Diffraction anisotropy** — the description of the overall fall-off by a single anisotropic
displacement tensor, its symmetry constraints, and the fact that only its deviatoric part is
determined (the isotropic part being degenerate with the overall scale) are Sheriff and Hendrickson's.
The estimator fits that tensor to the observed intensity distribution, taking sigma(I) into account,
in the sense of Popov and Bourenkov. The directional diffraction limits - <I/sigma(I)> in a cone about
each principal direction, and the reporting of the anisotropic deltaB as the range of the principal
components - follow AIMLESS. rugnux reports these; it corrects no intensity and removes no reflection
on a directional criterion. S. Sheriff and W. A. Hendrickson, "Description of overall anisotropy in
diffraction from macromolecular crystals" (1987), Acta Cryst. A43, 118-121
[doi:10.1107/S010876738709977X](https://doi.org/10.1107/S010876738709977X); A. N. Popov and
G. P. Bourenkov, "Choice of data-collection parameters based on statistic modelling" (2003), Acta
Cryst. D59, 1145-1153 [doi:10.1107/S0907444903008163](https://doi.org/10.1107/S0907444903008163);
P. R. Evans and G. N. Murshudov, "How good are my data and what is the resolution?" (2013), Acta
Cryst. D69, 1204-1214 [doi:10.1107/S0907444913000061](https://doi.org/10.1107/S0907444913000061).
**Data-quality statistics** follow the established conventions rather than any one program: R_meas
and R_pim, CC1/2 and CC\*, and the reporting of I/sigma(I). K. Diederichs and P. A. Karplus, "Improved
R-factors for diffraction data analysis in macromolecular crystallography" (1997), Nat. Struct. Biol.
docs/CHANGELOG.md
+1
View File
@@ -5,6 +5,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as
* The PCIe driver DKMS package builds for the kernel it is being installed for instead of the running one, so a module built while a kernel update is being applied loads after the reboot.
* The PCIe driver builds on RHEL 9.5 and later, and on their CentOS Stream, Rocky and AlmaLinux equivalents, where the `vm_flags` kernel interface was backported into the 5.14 kernel.
* rugnux reports diffraction anisotropy: the anisotropic deltaB and the diffraction limit along each
* rugnux: `--mode scale` reports the Wilson B-factor estimate instead of `WILSON_B= nan`.
* rugnux: `--export-unmerged` writes the integrated observations as `<prefix>_unmerged.mtz`, an unmerged MTZ readable by aimless, pointless, careless and `iotbx.merging_statistics`, in `--mode mx` and `--mode scale` alike; each rotation reflection's partials are summed into one full, and `--export-unmerged-partials` writes one row per image instead. Lattice-centring absences are not written; screw and glide absences are.
docs/CPU_DATA_ANALYSIS.md
+27 -1
View File
@@ -37,6 +37,8 @@ The methods draw on, and in places reimplement, solutions from:
- P. Evans, "Scaling and assessment of data quality", *Acta Cryst.* **D62** (2006), 72-82, and P. R. Evans, *Acta Cryst.* **D67** (2011), 282-292 (POINTLESS: operator-by-operator point-group scoring, and the axial-zone screw-absence test).
- A. G. W. Leslie & H. R. Powell, "Processing diffraction data with MOSFLM" (2007), NATO Science Series II **245**, 41-51 (post-refinement practice: what is refined per image and what over a wedge).
- D. W. Moreau, H. Atakisi & R. E. Thorne, "Ice in biomolecular cryocrystallography", *Acta Cryst.* **D77** (2021), 540-554 (measured hexagonal-ice ring positions, used by the ice-ring score, the ice flagging and the ice calibrant).
- S. Sheriff & W. A. Hendrickson, "Description of overall anisotropy in diffraction from macromolecular crystals", *Acta Cryst.* **A43** (1987), 118-121, and A. N. Popov & G. P. Bourenkov, *Acta Cryst.* **D59** (2003), 1145-1153 (the overall anisotropic B tensor, its symmetry constraints, and its estimation from the observed intensities).
- P. R. Evans & G. N. Murshudov, "How good are my data and what is the resolution?", *Acta Cryst.* **D69** (2013), 1204-1214 (AIMLESS: the anisotropic deltaB as the range of the principal components, and diffraction limits from a cone about each principal direction).
- K. Diederichs & P. A. Karplus, *Nat. Struct. Biol.* **4** (1997), 269-275, and P. A. Karplus & K. Diederichs, *Science* **336** (2012), 1030-1033 (R_meas / R_pim, CC1/2 and CC\*).
- IUCr Commission on Crystallographic Nomenclature, "Statistical descriptors in crystallography", *Acta Cryst.* **A45** (1989), 63-75, and *Acta Cryst.* **A51** (1995), 565-569 (uncertainty conventions).
@@ -958,7 +960,31 @@ Merging applies an optional per-observation median-based $N\sigma$ cut (`--rejec
By default the reported/written high-resolution limit is trimmed where $\mathrm{CC}_{1/2}$ falls off: a logistic is fitted to $\mathrm{CC}_{1/2}(s)$, and the limit is set **one reported-shell width past** the point where the fit crosses 0.30 — deliberately "one shell too far", so weak-but-real data below the crossing are kept rather than discarded. The extension is measured over the range that is actually kept, not the full measured range, so a detector reaching far past where the crystal diffracts cannot inflate it. `--scaling-high-resolution` overrides the limit and `--resolution-cutoff off` disables it.
### 13.5 Practical notes and limitations
### 13.5 Diffraction anisotropy
How fast the intensity falls off with resolution can depend on direction. rugnux measures that, reports it, and does nothing else with it: no intensity is corrected, no reflection is removed on a directional criterion, and the merged data and the written files do not depend on direction at all.
**The tensor.** A deviatoric anisotropic displacement tensor is fitted to the merged intensities as
$$\ln \langle I(\mathbf{s})\rangle = c(\text{shell}) - \tfrac{1}{2}\,\mathbf{s}^\mathsf{T} B\, \mathbf{s},\qquad \mathbf{s} = \text{reciprocal-space vector},\ |\mathbf{s}| = 1/d$$
with one free constant per resolution shell, so every isotropic feature — the Wilson curve, an ice ring, a noise floor, a scaling error — is absorbed exactly and only the $\ell = 2$ angular part drives the tensor. For an isotropic $B$ this reduces to the ordinary Wilson plot, so $B$ here is the ordinary crystallographic ($B = 8\pi^2 U$) $B$, directly comparable with phenix.xtriage's `B_cart`, ctruncate's anisotropic $B$ eigenvalues and AIMLESS's anisotropic $\Delta B$. Only the deviatoric part is fitted: the isotropic part is degenerate with the overall scale. The tensor is constrained to the directions the Laue class allows — five free deviatoric parameters in triclinic, three in monoclinic, two in orthorhombic, one in tetragonal, trigonal and hexagonal, and **none at all in cubic**, where symmetry forces $\Delta B$ to be exactly zero.
The fit is on **intensities, with no positivity cut**. Fitting amplitudes, or dropping non-positive intensities as an amplitude-based tool must, loses roughly 40% of the signal: in a direction that has died half the merged intensities are negative, so a positivity cut keeps only the positive noise excursions and flattens the fall-off exactly where the anisotropy is largest.
**Two different quantities are reported, and they are not interchangeable.** $\Delta B$ (the range of the principal components) is a *rate*; the diffraction limit along each principal direction — where $\langle I/\sigma(I)\rangle$ in a 20° cone about that direction falls through 2 — is where the signal actually runs out. A crystal can have a large $\Delta B$ and almost no spread in directional limit, or the reverse.
**The resolution signature.** A genuine Debye–Waller $B$ makes the directional deficit a straight line through the origin in $s^2$. The per-shell $\ell = 2$ amplitude is therefore fitted against $s^2$ and the curve is classified: *linear* (a real $B$), *flat* (a deficit that does not follow $\exp(-\tfrac12 \mathbf{s}^\mathsf{T} B \mathbf{s})$ at all, so the fitted $\Delta B$ describes the data with the wrong functional form and may be an **under**-estimate), or *convex* (a deficit that grows faster than $s^2$, which a $B$ cannot do). The verdict is re-derived at 8 and at 16 shells, and reported as undetermined if it moves.
**The verdict, and what it is measured against.** Whether an anisotropy is real is not decided against a counting-statistics error bar. Real data carry systematic error far larger than counting error, and gating on the latter reports anisotropy on datasets that have none. Instead the data set measures its own systematic error: in the tensor directions the Laue class *forbids*, the true tensor is exactly zero whatever the crystal is, so whatever is measured there is systematic. That measurement needs the unmerged observations — a merge has exact Laue symmetry by construction, and the forbidden directions are identically zero in it — so it is made on the scaled, rocking-curve-assembled observations. The counting part is subtracted, the counting error of the directions actually being tested is added back, and the ratio of $\Delta B$ to the result is banded: below 2 not established, 2–3.5 marginal, above 3.5 established, above 5 strong.
The report says **NOT DETECTED**, **DETECTED**, or **CANNOT DETERMINE**, and the third is a real answer rather than an evasion. It is returned when the Laue class is triclinic (no forbidden direction exists, so there is no internal measurement of the systematic error and no substitute for it), when the observed rotation range is under about 90° (a lab-fixed systematic then reaches several tensor directions instead of one), when the merged data are at the noise floor, when the scale model carried no dose term (an uncorrected dose ramp manufactures anisotropy that no significance test can see through), or when no unmerged observations were available. The smallest $\Delta B$ that could have been established on the data set is reported with the verdict; it is set by the systematic error rather than by counting, so it does **not** improve with more reflections or a longer exposure.
A too-high space-group assignment is the one failure mode that is silent: real anisotropy is then pushed into the directions used to measure the systematic error, which inflates the floor and biases the answer towards reporting none. A caution says so wherever the Laue class leaves a single free direction.
Everything lands in `<prefix>_report.txt` section 9 (`ANISOTROPY_*` keys), in the printed statistics, and in the merged mmCIF: the eigen-decomposition of the tensor as the standard `_reflns.pdbx_aniso_B_tensor_*` items (relative to the weakest direction, since only the deviatoric part is determined), and the directional limits, the shape and the verdict under the `_reflns.jfjoch_aniso_*` local prefix.
### 13.6 Practical notes and limitations
- **Bragg integration is profile-fitted by default** (per-shell Gaussian profile, Kabsch extraction; §9.3), with plain box summation available as a fallback (`--integrator boxsum`). The profiles are built per frame from that frame's strong spots, which suits fast-feedback and serial/streaming use; a profile shared across many frames (as in full offline workflows) is not currently formed.
- **Space-group symmetry** beyond centering absences is not enforced during prediction/integration unless the space group is supplied and used downstream.
image_analysis/WriteReflections.cpp
+42
View File
@@ -274,6 +274,48 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
if (std::isfinite(statistics.radiation_damage_delta_b))
out << "_reflns.jfjoch_radiation_damage_relative_B " << Fmt(statistics.radiation_damage_delta_b, 2)
<< " # relative-B first->last over the run (A^2); + = high-res fades with dose\n";
// Diffraction anisotropy. The eigen-decomposition of the anisotropy tensor has standard PDBx
// items; the eigenvalues there must be non-negative, and only the deviatoric part of the tensor
// is determined at all (its isotropic part is degenerate with the overall scale), so they are
// written relative to the weakest direction - eigenvalue_3 is 0 by construction and
// eigenvalue_1 is the anisotropic deltaB. The eigenvectors are in the PDB orthogonalisation
// convention, which is the one gemmi (and hence rugnux) uses throughout.
//
// The directional diffraction LIMITS are deliberately NOT written as
// _reflns.pdbx_aniso_diffraction_limit_*: the dictionary defines those as the semi-axes of an
// ellipsoid fitted to a diffraction cut-off surface, which is a different construction from the
// one below and one rugnux does not perform - it cuts nothing on a directional criterion. They
// go under the jfjoch local prefix with their own definition instead.
const auto &an = statistics.anisotropy;
if (an.n_reflections > 0 && an.n_cells > 0 && std::isfinite(an.delta_b)) {
out << "_reflns.pdbx_orthogonalization_convention pdb\n";
for (int i = 0; i < 3; ++i) {
out << "_reflns.pdbx_aniso_B_tensor_eigenvalue_" << (i + 1) << " "
<< Fmt(an.eigenvalue[i] - an.eigenvalue[2], 2) << "\n";
for (int j = 0; j < 3; ++j)
out << "_reflns.pdbx_aniso_B_tensor_eigenvector_" << (i + 1) << "_ortho[" << (j + 1)
<< "] " << Fmt(an.eigenvector[i][j], 4) << "\n";
}
out << "_reflns.jfjoch_aniso_delta_B " << Fmt(an.delta_b, 2)
<< " # range of the principal components (A^2), fitted on intensities\n";
if (std::isfinite(an.delta_b_linear))
out << "_reflns.jfjoch_aniso_delta_B_linear " << Fmt(an.delta_b_linear, 2)
<< " # deltaB implied by the s^2 slope alone: what a Debye-Waller B accounts for\n";
out << "_reflns.jfjoch_aniso_shape " << AnisotropyShapeCode(an.shape)
<< " # resolution signature of the directional deficit\n";
if (std::isfinite(an.floor))
out << "_reflns.jfjoch_aniso_floor " << Fmt(an.floor, 3)
<< " # deltaB this data set's own systematic error could manufacture (A^2)\n";
if (std::isfinite(an.significance))
out << "_reflns.jfjoch_aniso_significance " << Fmt(an.significance, 2)
<< " # deltaB(linear) / floor\n";
out << "_reflns.jfjoch_aniso_verdict " << AnisotropyVerdictCode(an.verdict) << "\n";
for (int i = 0; i < 3; ++i)
if (std::isfinite(an.d_min_axis[i]))
out << "_reflns.jfjoch_aniso_d_min_" << (i + 1) << " "
<< Fmt(an.d_min_axis[i], 2)
<< " # <I/sigma(I)> = 2 in a 20 deg cone about eigenvector " << (i + 1) << "\n";
}
out << "#\n";
// Per-batch relative-B curve (the radiation-damage monitor, rotation): one relative Debye-Waller B
image_analysis/scale_merge/AnisotropyAnalysis.cpp
+1357
View File
File diff suppressed because it is too large Load Diff
image_analysis/scale_merge/AnisotropyAnalysis.h
+154
View File
@@ -0,0 +1,154 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
#include <string>
#include <vector>
#include "../../common/Reflection.h"
#include "../IntegrationOutcome.h"
#include "gemmi/symmetry.hpp"
#include "gemmi/unitcell.hpp"
// Diffraction-anisotropy diagnostic. It DESCRIBES and it REPORTS: no intensity is corrected, no
// reflection is removed, and the written reflection file does not depend on direction in any way.
// Two quantities are reported, and they are not the same thing:
//
// - the anisotropic deltaB, the range of the principal components of the anisotropy tensor (the
// AIMLESS quantity), which is a RATE of fall-off;
// - the diffraction limit along each principal direction, from <I/sigma(I)> in a cone about that
// direction (the AIMLESS construction), which is where the signal actually runs out.
//
// A crystal can have a large deltaB and almost no spread in directional limit, or the reverse.
//
// CONVENTION. The tensor is fitted as ln <I(s)> = c(shell) - 1/2 s^T B s, s = 1/d in A^-1, B in A^2.
// For an isotropic B this is the ordinary Wilson plot, so B here is the ordinary crystallographic
// (structure-factor-level, B = 8 pi^2 U) B - the same scale as phenix.xtriage's B_cart, ctruncate's
// anisotropic B eigenvalues and AIMLESS's anisotropic deltaB, and directly comparable with all three.
// Only the DEVIATORIC part is fitted: the isotropic part is degenerate with the overall scale and with
// the Wilson curve, and one free constant per resolution shell absorbs it exactly.
//
// The tensor is fitted on INTENSITIES, with no positivity cut. Fitting amplitudes, or dropping
// non-positive intensities as the amplitude-based tools must, loses roughly 40% of the signal: in a
// direction that has died, half the merged intensities are negative, so a positivity cut keeps only
// the positive noise excursions and flattens the fall-off exactly where the anisotropy is largest.
//
// TWO INPUTS, because they answer different halves of the question. The tensor, its resolution
// signature and the directional limits are measured on the MERGED reflections. The error bar the
// verdict is gated on cannot be: merged data has exact Laue symmetry by construction, so the
// symmetry-forbidden tensor directions - where the true tensor is zero whatever the crystal is, and
// which are therefore this dataset's own measurement of its systematic error - are identically zero
// there. That measurement needs the UNMERGED, scaled observations, which still carry the differences
// between symmetry mates. Without them the verdict is "cannot determine" and says so.
//
// Method credits are at the algorithms in the .cpp.
enum class AnisotropyShape {
Undetermined, // too few shells, or the verdict moved when the binning changed
Linear, // the deficit follows exp(-1/2 s^T B s): a Debye-Waller B
Flat, // the deficit does not follow exp(-1/2 s^T B s) - see the note on Flat below
Convex // the deficit grows faster than s^2
};
enum class AnisotropyVerdict {
CannotDetermine, // this dataset cannot answer the question - `refusal` says why
NotDetected, // anisotropy is not established above this dataset's own systematic error
Detected // anisotropy is established; `band` grades it
};
const char *AnisotropyShapeCode(AnisotropyShape shape); // e.g. "LINEAR"
const char *AnisotropyVerdictCode(AnisotropyVerdict verdict); // e.g. "NOT_DETECTED"
struct AnisotropyResult {
// "Could not decide" is a counter, not a flag: every renderer tests n_reflections and falls silent
// when it is zero (the TwinningAnalysis idiom). n_free_parameters is the deviatoric degrees of
// freedom the Laue class allows - 5 triclinic, 3 monoclinic, 2 orthorhombic, 1 tetragonal/trigonal/
// hexagonal, 0 cubic - so 0 means "no anisotropy is possible here", not "nothing was measured", and
// a cubic crystal reports n_reflections > 0 with n_cells = 0 and deltaB exactly zero.
int n_reflections = 0;
int n_cells = 0;
int n_free_parameters = 0;
int n_observations = 0; // unmerged observations the systematic-error floor was measured on
// --- the tensor ---------------------------------------------------------------------------
double delta_b = NAN; // B_max - B_min of the deviatoric tensor, A^2
double eigenvalue[3] = {NAN, NAN, NAN}; // descending; deviatoric, so they sum to zero
double eigenvector[3][3] = {{NAN, NAN, NAN}, {NAN, NAN, NAN}, {NAN, NAN, NAN}}; // Cartesian rows
double fold_weakening = NAN; // exp(delta_b / 2 d_min^2): strong/weak intensity ratio at d_min
// --- the diffraction limits (model-free) --------------------------------------------------
// Highest resolution at which <I/sigma(I)> in a 20 deg cone about principal axis n is still 2.0.
double d_min_axis[3] = {NAN, NAN, NAN};
double d_min_spread = NAN; // max - min of the above
// --- the resolution signature of the deficit ----------------------------------------------
// The l=2 directional deficit A of each shell against s^2, fitted as A = c0 + c1 s^2. A genuine
// Debye-Waller B gives a straight line through the origin (c0 = 0, c1 = |B_dev|).
int shape_shells = 0;
double shape_intercept = NAN; // c0, dimensionless
double shape_intercept_z = NAN; // c0 / sigma(c0)
double shape_slope = NAN; // c1, A^2
double shape_curvature_z = NAN; // significance of an s^4 term
double shape_curvature_share = NAN; // fraction of the deficit that term carries at the limit
double shape_residual = NAN; // chi2/dof of the through-origin line
AnisotropyShape shape = AnisotropyShape::Undetermined;
bool shape_stable = false; // the verdict survived rebinning at 8 and 16 shells
// deltaB implied by the s^2 slope alone - what a Debye-Waller B accounts for - and the rest.
// The first can exceed delta_b when the fitted curve passes below the origin.
double delta_b_linear = NAN;
double delta_b_flat = NAN;
// --- the gate -----------------------------------------------------------------------------
// Everything here is measured from this dataset alone. sigma_systematic is the per-component
// systematic error scale left in the symmetry-FORBIDDEN tensor directions of the unmerged
// observations after their own counting noise has been taken out; floor is what that systematic,
// plus counting noise, would manufacture in the symmetry-ALLOWED directions; significance is
// delta_b_linear / floor.
double sigma_systematic = NAN;
double floor = NAN;
double significance = NAN;
double detection_limit = NAN; // the smallest delta_b this dataset could establish, A^2
AnisotropyVerdict verdict = AnisotropyVerdict::CannotDetermine;
std::string band; // "not established" / "marginal" / "established" / "strong"
std::string refusal; // why, when the verdict is CannotDetermine
std::vector<std::string> cautions; // things that bias the verdict, each measured
};
// One scaled, unmerged observation: the Miller index it was measured at (NOT reduced to the
// asymmetric unit - the whole point is that symmetry mates are separate measurements) and its
// intensity on the merge's own scale.
struct AnisotropyObservation {
int32_t h = 0, k = 0, l = 0;
float I = NAN;
float sigma = NAN;
float d = NAN;
};
// Build those from the per-image integration outcomes, on the scale rugnux itself applied. On rotation
// data each reflection's partials are first assembled into one full, as the 3D combine and the
// unmerged export do; on stills each reflection is already a whole measurement. Images with no fitted
// per-image scale are left out, and so is an event that caught less than min_partiality of its rocking
// curve - the combine does not merge one either.
std::vector<AnisotropyObservation> ScaledObservations(const std::vector<IntegrationOutcome> &outcomes,
bool rotation, double min_partiality = 0.5);
// What the caller knows about the run and the merge that the reflections alone do not say.
struct AnisotropyRunInfo {
// Rotation actually covered by the processed images, in degrees (NOT the header's nominal sweep -
// one battery dataset spans 179 deg against a 360 deg header). NaN for stills / not known.
double observed_rotation_deg = NAN;
// Whether the scale model that produced these intensities carried a dose (decay) term. Mandatory
// on rotation data: see the note at the refusal in the .cpp.
bool dose_term_in_scale_model = true;
// The run's measured radiation damage (MergeStatistics::radiation_damage_delta_b), for the caution.
double radiation_damage_relative_b = NAN;
};
AnisotropyResult AnalyzeAnisotropy(const std::vector<MergedReflection> &merged,
const std::vector<AnisotropyObservation> &unmerged,
const gemmi::UnitCell &cell,
const gemmi::SpaceGroup *space_group,
const AnisotropyRunInfo &run = {});
std::string AnisotropyToText(const AnisotropyResult &result);
image_analysis/scale_merge/CMakeLists.txt
+2
View File
@@ -3,6 +3,8 @@ ADD_LIBRARY(JFJochScaleMerge
SearchSpaceGroup.h
TwinningAnalysis.cpp
TwinningAnalysis.h
AnisotropyAnalysis.cpp
AnisotropyAnalysis.h
Merge.cpp
Merge.h
ScaleOnTheFly.cpp
image_analysis/scale_merge/Merge.h
+7
View File
@@ -14,6 +14,7 @@
#include "../../common/Reflection.h"
#include "../IntegrationOutcome.h"
#include "AnisotropyAnalysis.h"
#include "HKLKey.h"
struct MergeStatisticsShell {
@@ -104,6 +105,12 @@ struct MergeStatistics {
// Stretches of the sweep over which the crystal delivered much less than the rest of the run
// (MeasureSweepQuality). Report-only - no observation is dropped because of it.
SweepQuality sweep_quality;
// Diffraction anisotropy (AnalyzeAnisotropy): the anisotropy tensor, the diffraction limit along
// each of its principal directions, and whether either is established above this dataset's own
// systematic error. Report-only - nothing is corrected and no reflection is removed. Empty
// (n_reflections = 0) when the diagnostic did not run.
AnisotropyResult anisotropy;
};
rugnux/ResultReport.cpp
+45 -3
View File
@@ -1,6 +1,7 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <algorithm>
#include <fstream>
#include <sstream>
@@ -8,6 +9,7 @@
#include "../common/GitInfo.h"
#include "../common/time_utc.h"
#include "../image_analysis/scale_merge/AnisotropyAnalysis.h"
#include "../image_analysis/scale_merge/Merge.h"
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
@@ -17,7 +19,7 @@
namespace {
// The version of this file format. Bumped when a key is renamed or removed, a table column moves,
// or a reason code changes meaning - a consumer can gate on it.
constexpr int REPORT_VERSION = 1;
constexpr int REPORT_VERSION = 2;
const char *BANNER = " ******************************************************************************";
@@ -263,11 +265,51 @@ std::string RenderResultReport(const std::string &output_prefix,
}
os << " ----------- ----------- --------- -------- -------------------- -------- ------ ------ --------\n";
// --------------------------------------------------------------- 9. WARNINGS
// ---------------------------------------------------------- 9. DIFFRACTION ANISOTROPY
const auto &an = result.merge_statistics.anisotropy;
if (merged && an.n_reflections > 0) {
Section(os, "9. DIFFRACTION ANISOTROPY");
os << " How much the fall-off depends on direction, and whether that is established above this\n"
<< " data set's own systematic error. Nothing here corrects an intensity or removes a\n"
<< " reflection: the merged data and the written files do not depend on direction at all.\n"
<< " ANISOTROPY_DELTA_B is the range of the principal components of the anisotropy tensor,\n"
<< " on the ordinary crystallographic B scale (the same scale as phenix.xtriage's B_cart and\n"
<< " ctruncate's anisotropic B), fitted on intensities with nothing dropped.\n\n";
Key(os, "ANISOTROPY_VERDICT", AnisotropyVerdictCode(an.verdict));
Key(os, "ANISOTROPY_FREE_DIRECTIONS", an.n_free_parameters);
Key(os, "ANISOTROPY_DELTA_B", fmt::format("{:.2f}", an.delta_b));
Key(os, "ANISOTROPY_DELTA_B_LINEAR", fmt::format("{:.2f}", an.delta_b_linear));
Key(os, "ANISOTROPY_PRINCIPAL_B", fmt::format("{:.2f} {:.2f} {:.2f}",
an.eigenvalue[0] - an.eigenvalue[2],
an.eigenvalue[1] - an.eigenvalue[2], 0.0));
Key(os, "ANISOTROPY_FOLD_WEAKENING", fmt::format("{:.1f}", an.fold_weakening));
Key(os, "ANISOTROPY_D_MIN_PRINCIPAL", fmt::format("{:.2f} {:.2f} {:.2f}", an.d_min_axis[0],
an.d_min_axis[1], an.d_min_axis[2]));
Key(os, "ANISOTROPY_D_MIN_SPREAD", fmt::format("{:.2f}", an.d_min_spread));
Key(os, "ANISOTROPY_SHAPE", AnisotropyShapeCode(an.shape));
Key(os, "ANISOTROPY_SHAPE_INTERCEPT", fmt::format("{:.3f}", an.shape_intercept));
Key(os, "ANISOTROPY_SHAPE_INTERCEPT_Z", fmt::format("{:.1f}", an.shape_intercept_z));
Key(os, "ANISOTROPY_SHAPE_SLOPE", fmt::format("{:.2f}", an.shape_slope));
Key(os, "ANISOTROPY_SHAPE_RESIDUAL", fmt::format("{:.1f}", an.shape_residual));
Key(os, "ANISOTROPY_SIGMA_SYSTEMATIC", fmt::format("{:.3f}", an.sigma_systematic));
Key(os, "ANISOTROPY_FLOOR", fmt::format("{:.3f}", an.floor));
Key(os, "ANISOTROPY_SIGNIFICANCE", fmt::format("{:.2f}", an.significance));
Key(os, "ANISOTROPY_DETECTION_LIMIT", fmt::format("{:.2f}", an.detection_limit));
os << "\n" << AnisotropyToText(an) << "\n";
if (an.verdict == AnisotropyVerdict::Detected && an.d_min_spread > 0.5)
warnings.emplace_back(fmt::format(
"Diffraction is anisotropic (deltaB {:.1f} A^2; the diffraction limit runs from "
"{:.2f} to {:.2f} A depending on direction) - refinement and map interpretation "
"should allow for it; no intensity has been corrected for it here",
an.delta_b, *std::max_element(an.d_min_axis, an.d_min_axis + 3),
*std::min_element(an.d_min_axis, an.d_min_axis + 3)));
}
// --------------------------------------------------------------- 10. WARNINGS
if (result.cancelled)
warnings.emplace_back(fmt::format("Processing was cancelled after {} images - this report "
"describes an incomplete run", result.images_processed));
Section(os, "9. WARNINGS");
Section(os, "10. WARNINGS");
os << " Everything that needs a person's attention, one line each, marked so a script can find\n"
<< " them with a single grep for \"WARNING:\".\n\n";
Key(os, "WARNING_COUNT", warnings.size());
rugnux/Rugnux.cpp
+22
View File
@@ -47,6 +47,7 @@
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
#include "../image_analysis/geom_refinement/PostRefine.h"
#include "../image_analysis/lattice_search/LatticeSearch.h"
#include "../image_analysis/scale_merge/AnisotropyAnalysis.h"
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
#include "../image_analysis/scale_merge/HKLKey.h"
#include "../image_analysis/scale_merge/ScaleOnTheFly.h"
@@ -2815,6 +2816,27 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
if (std::isfinite(wilson.b) && wilson.b > 0.0)
logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells)",
wilson.b, wilson.correlation, wilson.n_shells);
// Diffraction anisotropy. Report-only, like the twinning and Wilson analyses above: it
// describes the data and corrects nothing. The tensor and its resolution signature come from
// the merged intensities; the error bar the verdict is gated on has to come from the
// unmerged observations, because a merge has exact Laue symmetry by construction and the
// tensor directions the symmetry forbids - the only place a dataset measures its own
// systematic error - are identically zero in it.
if (result.consensus_cell) {
AnisotropyRunInfo aniso_run;
if (sm.statistics.sweep_quality.measured && sm.statistics.sweep_quality.sweep_deg > 0.0f)
aniso_run.observed_rotation_deg = sm.statistics.sweep_quality.sweep_deg;
aniso_run.dose_term_in_scale_model =
experiment_.GetScalingSettings().GetCorrectionSurfaces();
aniso_run.radiation_damage_relative_b = sm.statistics.radiation_damage_delta_b;
sm.statistics.anisotropy = AnalyzeAnisotropy(
sm.merged,
ScaledObservations(indexer->GetIntegrationOutcome(),
experiment_.IsRotationIndexing()),
*result.consensus_cell, twin_sg, aniso_run);
stats_text << AnisotropyToText(sm.statistics.anisotropy) << "\n";
}
}
stats_text << sm.statistics;
rugnux/rugnux_cli.cpp
+16
View File
@@ -36,6 +36,7 @@
#include "../image_analysis/scale_merge/StillsPartialityRefine.h"
#include "../image_analysis/scale_merge/RotationScaleMerge.h"
#include "../image_analysis/scale_merge/ResolutionCutoff.h"
#include "../image_analysis/scale_merge/AnisotropyAnalysis.h"
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
#include "Rugnux.h"
@@ -1552,6 +1553,21 @@ static int RunRugnux(int argc, char **argv) {
const auto twinning = AnalyzeTwinning(merged_reflections, twin_sg);
std::cout << std::endl << TwinningAnalysisToText(twinning) << std::endl;
// Diffraction anisotropy, exactly as the full pipeline computes it (Rugnux.cpp) - report-only,
// nothing is corrected. The stored reflections have just been re-scaled above, so they carry
// this merge's own per-image scale.
if (experiment.GetUnitCell()) {
AnisotropyRunInfo aniso_run;
if (merged_statistics.sweep_quality.measured && merged_statistics.sweep_quality.sweep_deg > 0.0f)
aniso_run.observed_rotation_deg = merged_statistics.sweep_quality.sweep_deg;
aniso_run.dose_term_in_scale_model = experiment.GetScalingSettings().GetCorrectionSurfaces();
aniso_run.radiation_damage_relative_b = merged_statistics.radiation_damage_delta_b;
merged_statistics.anisotropy = AnalyzeAnisotropy(merged_reflections,
ScaledObservations(reflections, is_rotation),
*experiment.GetUnitCell(), twin_sg, aniso_run);
std::cout << AnisotropyToText(merged_statistics.anisotropy) << std::endl;
}
// Unmerged observations (--export-unmerged), from the integrated observations rather than the
// merged ones: the partiality and the per-image scale are left for the reading program, which
// fits a scale model of its own.
tests/AnisotropyAnalysisTest.cpp
+153
View File
@@ -0,0 +1,153 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <catch2/catch_all.hpp>
#include <cmath>
#include <vector>
#include "../image_analysis/scale_merge/AnisotropyAnalysis.h"
#include "gemmi/symmetry.hpp"
#include "gemmi/unitcell.hpp"
namespace {
gemmi::UnitCell Cell(double a, double b, double c, double al, double be, double ga) {
gemmi::UnitCell out;
out.set(a, b, c, al, be, ga);
return out;
}
// A synthetic merge: Wilson-distributed intensities with an isotropic fall-off and a known
// deviatoric anisotropy on top, over every hkl inside the resolution limit. The "structure
// factor" is a hash of the index, so the set is reproduced bit for bit.
std::vector<MergedReflection> SyntheticMerge(const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg,
double d_min, double b_iso,
const gemmi::SMat33<double> &b_dev) {
const gemmi::GroupOps ops = sg->operations();
const gemmi::ReciprocalAsu asu(sg);
std::vector<MergedReflection> out;
const int hmax = static_cast<int>(cell.a / d_min) + 1;
const int kmax = static_cast<int>(cell.b / d_min) + 1;
const int lmax = static_cast<int>(cell.c / d_min) + 1;
for (int h = -hmax; h <= hmax; ++h)
for (int k = -kmax; k <= kmax; ++k)
for (int l = -lmax; l <= lmax; ++l) {
const gemmi::Op::Miller hkl{{h, k, l}};
if ((h == 0 && k == 0 && l == 0) || !asu.is_in(hkl) || ops.is_systematically_absent(hkl))
continue;
const double d2 = cell.calculate_1_d2(hkl);
if (d2 <= 0 || 1.0 / std::sqrt(d2) < d_min)
continue;
// Wilson draw from a hash of the index: deterministic, and spanning a realistic range.
const uint32_t seed = static_cast<uint32_t>(h * 73856093 ^ k * 19349663 ^ l * 83492791);
const double u = ((seed * 2654435761u) >> 8) / static_cast<double>(1 << 24);
const double wilson = -std::log(std::max(1e-6, u));
// The tensor is applied in the Cartesian frame, s = F^T h.
const gemmi::Vec3 s = cell.frac.mat.left_multiply(gemmi::Vec3(h, k, l));
const double aniso = -0.5 * b_dev.r_u_r(s);
MergedReflection r;
r.h = h; r.k = k; r.l = l;
r.d = static_cast<float>(1.0 / std::sqrt(d2));
r.I = static_cast<float>(1000.0 * wilson * std::exp(-0.5 * b_iso * d2 + aniso));
r.sigma = static_cast<float>(0.02 * std::fabs(r.I) + 1.0);
out.push_back(r);
}
return out;
}
}
// The number of free deviatoric anisotropy parameters is fixed by the Laue class alone. This is the
// self-test of the constraint basis: 5 / 3 / 2 / 1 / 1 / 0 for triclinic / monoclinic / orthorhombic /
// tetragonal / trigonal-hexagonal / cubic, and nothing else is possible.
TEST_CASE("Anisotropy free-parameter count", "[anisotropy]") {
struct Case {
const char *space_group;
gemmi::UnitCell cell;
int free_directions;
};
const std::vector<Case> cases{
{"P 1", Cell(51, 62, 73, 84.0, 95.0, 103.0), 5},
{"P 1 21 1", Cell(51, 62, 73, 90.0, 95.0, 90.0), 3},
{"C 1 2 1", Cell(91, 62, 73, 90.0, 105.0, 90.0), 3},
{"P 21 21 21", Cell(51, 62, 73, 90.0, 90.0, 90.0), 2},
{"I 2 2 2", Cell(51, 62, 73, 90.0, 90.0, 90.0), 2},
{"P 43 21 2", Cell(79, 79, 38, 90.0, 90.0, 90.0), 1}, // lysozyme, the field's test specimen
{"P 31 2 1", Cell(62, 62, 91, 90.0, 90.0, 120.0), 1},
{"R 3 :H", Cell(78, 78, 33, 90.0, 90.0, 120.0), 1},
{"P 63", Cell(62, 62, 91, 90.0, 90.0, 120.0), 1},
{"I 2 3", Cell(78, 78, 78, 90.0, 90.0, 90.0), 0},
{"F 4 3 2", Cell(78, 78, 78, 90.0, 90.0, 90.0), 0},
};
for (const auto &c : cases) {
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(c.space_group);
REQUIRE(sg != nullptr);
std::vector<MergedReflection> merged(1);
merged[0].h = 1; merged[0].k = 0; merged[0].l = 0;
merged[0].I = 1.0f; merged[0].sigma = 1.0f; merged[0].d = 10.0f;
const auto result = AnalyzeAnisotropy(merged, {}, c.cell, sg);
INFO(c.space_group);
CHECK(result.n_free_parameters == c.free_directions);
}
}
// A refined cell need not obey its space group's metric constraints exactly, and rugnux writes the
// unconstrained refined cell. An angle a hundredth of a degree off 90 must not create an extra
// anisotropy direction.
TEST_CASE("Anisotropy free-parameter count with an off-metric refined cell", "[anisotropy]") {
std::vector<MergedReflection> merged(1);
merged[0].h = 1; merged[0].k = 0; merged[0].l = 0;
merged[0].I = 1.0f; merged[0].sigma = 1.0f; merged[0].d = 10.0f;
CHECK(AnalyzeAnisotropy(merged, {}, Cell(51, 62, 73, 90.02, 105.0, 89.97),
gemmi::find_spacegroup_by_name("C 1 2 1")).n_free_parameters == 3);
CHECK(AnalyzeAnisotropy(merged, {}, Cell(51.0, 62.0, 73.0, 89.98, 90.03, 90.01),
gemmi::find_spacegroup_by_name("I 2 2 2")).n_free_parameters == 2);
CHECK(AnalyzeAnisotropy(merged, {}, Cell(78.01, 77.99, 78.02, 90.01, 89.99, 90.0),
gemmi::find_spacegroup_by_name("I 2 3")).n_free_parameters == 0);
}
// A cubic crystal has no free deviatoric parameter, so its deltaB is exactly zero by symmetry - not
// small, not measured, zero - and the verdict is a statement about symmetry rather than about data.
TEST_CASE("Anisotropy is exactly zero in a cubic Laue class", "[anisotropy]") {
const gemmi::UnitCell cell = Cell(78, 78, 78, 90, 90, 90);
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("I 2 3");
const auto merged = SyntheticMerge(cell, sg, 2.5, 20.0, {8.0, -4.0, -4.0, 0.0, 0.0, 0.0});
REQUIRE(merged.size() > 1000);
const auto result = AnalyzeAnisotropy(merged, {}, cell, sg);
CHECK(result.n_free_parameters == 0);
CHECK(result.delta_b == 0.0);
CHECK(result.verdict == AnisotropyVerdict::NotDetected);
}
// The tensor itself: put a known deviatoric B into a tetragonal merge and read it back. The
// tetragonal Laue class leaves one free direction, along c*, and its magnitude is what deltaB means.
TEST_CASE("Anisotropy tensor is recovered from a synthetic merge", "[anisotropy]") {
const gemmi::UnitCell cell = Cell(79, 79, 38, 90, 90, 90);
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 43 21 2");
// Uniaxial about c, deltaB = B_zz - B_xx = 15 A^2.
const gemmi::SMat33<double> b_dev{-5.0, -5.0, 10.0, 0.0, 0.0, 0.0};
const auto merged = SyntheticMerge(cell, sg, 2.0, 20.0, b_dev);
REQUIRE(merged.size() > 2000);
const auto result = AnalyzeAnisotropy(merged, {}, cell, sg);
REQUIRE(result.n_free_parameters == 1);
REQUIRE(result.n_cells > 0);
CHECK(result.delta_b == Catch::Approx(15.0).margin(1.5));
// c* is the weak direction here, so the largest principal value points along z.
CHECK(std::fabs(result.eigenvector[0][2]) == Catch::Approx(1.0).margin(0.05));
// A pure Debye-Waller fall-off is a straight line through the origin in s^2.
CHECK(result.shape == AnisotropyShape::Linear);
// With no unmerged observations the systematic-error scale cannot be measured, and the verdict
// says so rather than falling back on a counting-statistics error bar.
CHECK(result.verdict == AnisotropyVerdict::CannotDetermine);
CHECK_FALSE(result.refusal.empty());
}
// An isotropic merge must not produce a tensor, whatever the Laue class allows.
TEST_CASE("Anisotropy of an isotropic synthetic merge is small", "[anisotropy]") {
const gemmi::UnitCell cell = Cell(79, 79, 38, 90, 90, 90);
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 43 21 2");
const auto merged = SyntheticMerge(cell, sg, 2.0, 20.0, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
REQUIRE(merged.size() > 2000);
const auto result = AnalyzeAnisotropy(merged, {}, cell, sg);
REQUIRE(result.n_cells > 0);
CHECK(std::fabs(result.delta_b) < 1.0);
}
tests/CMakeLists.txt
+1
View File
@@ -89,6 +89,7 @@ ADD_EXECUTABLE(jfjoch_test
SyntheticMergedReflections.h
XDSPluginTest.cpp
MergeScaleTest.cpp
AnisotropyAnalysisTest.cpp
RfreeFlagsTest.cpp
FrenchWilsonTest.cpp
ReindexAmbiguityTest.cpp