leonarski_fandClaude Opus 4.8
4aab8078d6
SearchSpaceGroup: drop the log10(chi2_ref) widening of the point-group chi^2-ratio bound and tighten max_merge_chi2_ratio 2.0 -> 1.85. The variance-floor fix removed the error-model miscalibration the widening compensated for, so the widening now only let a partial merohedral twin through (Ins_H_2 R3->R32 twin 2-fold: ratio 2.01). Every genuine high symmetry across the rotation-test battery stays within ~1.7x (worst real case Thau P41212 at 1.71), so 1.85 keeps R3 in its true lower symmetry. TwinningAnalysis: make the <I^2>/<I>^2 second moment robust - skip noise-only shells (<I/sigma> < 1) and reject Wilson outliers (E^2 > 8) with one shell-mean re-iteration, so a single strong reflection in a collapsed-mean shell no longer dominates the moment. Add MerohedralTwinningPossible: in a holohedral Laue class (4/mmm, 6/mmm, m-3m, rhombohedral -3m) no twin law exists, so a low <|L|> there is a statistical artefact and is no longer flagged as twinning. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
227 lines
11 KiB
C++
227 lines
11 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "TwinningAnalysis.h"
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <cstdint>
|
|
#include <iomanip>
|
|
#include <limits>
|
|
#include <sstream>
|
|
#include <unordered_map>
|
|
#include <vector>
|
|
|
|
namespace {
|
|
int64_t PackHKL(int h, int k, int l) {
|
|
constexpr int64_t bias = 1 << 20; // indices assumed within +/- 2^20
|
|
return ((h + bias) << 42) | ((k + bias) << 21) | (l + bias);
|
|
}
|
|
|
|
bool UsableIntensity(const MergedReflection& r) {
|
|
return std::isfinite(r.I) && std::isfinite(r.d) && r.d > 0.0;
|
|
}
|
|
|
|
// Merohedral twinning needs a twin law - a lattice symmetry operation that is not a symmetry of the
|
|
// crystal - which exists only when the Laue class is a proper subgroup of the lattice holohedry.
|
|
// The holohedral high-symmetry Laue classes (4/mmm, 6/mmm, m-3m, and -3m on a rhombohedral lattice)
|
|
// admit no such operation, so twinning is geometrically impossible and the intensity statistics
|
|
// cannot be indicating it. Low-symmetry classes stay eligible because pseudo-merohedral twinning
|
|
// through an accidental metric specialisation cannot be excluded from the symmetry alone.
|
|
bool MerohedralTwinningPossible(const gemmi::SpaceGroup* sg) {
|
|
if (!sg)
|
|
return true; // P1 / unknown symmetry: cannot rule twinning out
|
|
switch (sg->laue_class()) {
|
|
case gemmi::Laue::L4mmm: // 4/mmm - tetragonal holohedry
|
|
case gemmi::Laue::L6mmm: // 6/mmm - hexagonal holohedry
|
|
case gemmi::Laue::Lm3m: // m-3m - cubic holohedry
|
|
return false;
|
|
case gemmi::Laue::L3m: // -3m is holohedral on a rhombohedral (R) lattice, but a
|
|
return sg->hm[0] != 'R'; // hexagonal-P 32/3m crystal can still twin towards 6/mmm
|
|
default:
|
|
return true;
|
|
}
|
|
}
|
|
}
|
|
|
|
TwinningAnalysisResult AnalyzeTwinning(const std::vector<MergedReflection>& merged,
|
|
const gemmi::SpaceGroup* space_group, int resolution_shells) {
|
|
TwinningAnalysisResult result;
|
|
if (merged.empty())
|
|
return result;
|
|
|
|
// Centric reflections follow different statistics and must be excluded. In P1 (no space group)
|
|
// none are centric.
|
|
const gemmi::GroupOps gops = space_group ? space_group->operations() : gemmi::GroupOps{};
|
|
auto acentric = [&](const MergedReflection& r) {
|
|
return !space_group || !gops.is_reflection_centric(gemmi::Op::Miller{{r.h, r.k, r.l}});
|
|
};
|
|
|
|
// --- L-test ---
|
|
// Pair each reflection with a symmetry-independent neighbour two steps away along an axis (the
|
|
// step of 2 keeps the partner local in resolution while avoiding the reflection itself). The
|
|
// merged reflections are unique in the asymmetric unit, so any other merged reflection is
|
|
// genuinely non-equivalent - exactly the pairing the L-test wants. Only acentric reflections
|
|
// with positive intensity enter, which also keeps L = (I1-I2)/(I1+I2) bounded in [-1, 1].
|
|
std::unordered_map<int64_t, double> intensity;
|
|
intensity.reserve(merged.size() * 2);
|
|
for (const auto& r : merged)
|
|
if (UsableIntensity(r) && r.I > 0.0 && acentric(r))
|
|
intensity.emplace(PackHKL(r.h, r.k, r.l), r.I);
|
|
|
|
// Step to a nearby, symmetry-independent partner. The axis step of 2 preserves the reflection
|
|
// condition of P/I/C/F/A/B lattices (parity-based), but violates R-centring (-h+k+l = 0 mod 3,
|
|
// as 2 != 0 mod 3) -> the partner is systematically absent and rhombohedral crystals yield zero
|
|
// pairs. The diagonal (1,1,0)/(1,1,3) steps preserve the mod-3 condition in BOTH obverse and
|
|
// reverse settings; they are tried only when the axis steps find no present partner, so P/I/C/F
|
|
// behaviour is unchanged (first present partner wins).
|
|
const std::array<std::array<int, 3>, 5> steps{
|
|
{{2, 0, 0}, {0, 2, 0}, {0, 0, 2}, {1, 1, 0}, {1, 1, 3}}};
|
|
|
|
double sum_abs_l = 0.0, sum_l2 = 0.0;
|
|
int n_pairs = 0;
|
|
for (const auto& [key, i1] : intensity) {
|
|
const int h = static_cast<int>((key >> 42) & 0x1FFFFF) - (1 << 20);
|
|
const int k = static_cast<int>((key >> 21) & 0x1FFFFF) - (1 << 20);
|
|
const int l = static_cast<int>(key & 0x1FFFFF) - (1 << 20);
|
|
for (const auto& s : steps) {
|
|
const auto it = intensity.find(PackHKL(h + s[0], k + s[1], l + s[2]));
|
|
if (it == intensity.end())
|
|
continue;
|
|
const double lstat = (i1 - it->second) / (i1 + it->second);
|
|
sum_abs_l += std::fabs(lstat);
|
|
sum_l2 += lstat * lstat;
|
|
++n_pairs;
|
|
break; // one neighbour per reflection
|
|
}
|
|
}
|
|
result.l_test_pairs = n_pairs;
|
|
if (n_pairs > 0) {
|
|
result.mean_abs_l = sum_abs_l / n_pairs;
|
|
result.mean_l_squared = sum_l2 / n_pairs;
|
|
}
|
|
|
|
// --- Second moment <I^2>/<I>^2 of acentric intensities, normalised per resolution shell ---
|
|
// Binning by 1/d^2 removes the resolution fall-off, so the moment is 2.0 (untwinned) or 1.5
|
|
// (perfect twin) regardless of the overall B-factor. The moment divides by the *square* of the
|
|
// shell-mean intensity, so it is not robust: on weak or mis-integrated data a shell mean can
|
|
// collapse to the noise floor and one outlier reflection then dominates (I/mean)^2 (a single
|
|
// I=158 in a mean~1 shell contributed 78% of a whole dataset's value). To keep this a twinning
|
|
// indicator rather than a data-quality artefact - as phenix.xtriage does - we skip noise-only
|
|
// shells (<I/sigma> below 1) and reject Wilson outliers (E^2 above 8, ~exp(-8) upper tail) with
|
|
// one shell-mean re-iteration so the outlier does not corrupt the normalising mean either.
|
|
constexpr double min_shell_isig = 1.0; // shells below this are noise, not signal
|
|
constexpr double wilson_outlier_e2 = 8.0; // reject improbably strong reflections (P ~ e^-8)
|
|
int n_shells = std::max(1, resolution_shells);
|
|
double min_s = std::numeric_limits<double>::infinity();
|
|
double max_s = -std::numeric_limits<double>::infinity();
|
|
for (const auto& r : merged) {
|
|
if (!UsableIntensity(r) || !acentric(r))
|
|
continue;
|
|
const double s = 1.0 / (r.d * r.d);
|
|
min_s = std::min(min_s, s);
|
|
max_s = std::max(max_s, s);
|
|
}
|
|
if (std::isfinite(min_s) && max_s > min_s) {
|
|
auto shell_of = [&](double d) {
|
|
const double t = (1.0 / (d * d) - min_s) / (max_s - min_s);
|
|
return std::min(n_shells - 1, std::max(0, static_cast<int>(t * n_shells)));
|
|
};
|
|
// Group acentric intensities by shell, and accumulate <I/sigma> to gauge each shell's signal.
|
|
std::vector<std::vector<double>> shell_I(n_shells);
|
|
std::vector<double> shell_isig_sum(n_shells, 0.0);
|
|
std::vector<int> shell_isig_n(n_shells, 0);
|
|
for (const auto& r : merged) {
|
|
if (!UsableIntensity(r) || !acentric(r))
|
|
continue;
|
|
const int b = shell_of(r.d);
|
|
shell_I[b].push_back(r.I);
|
|
if (std::isfinite(r.sigma) && r.sigma > 0.0) {
|
|
shell_isig_sum[b] += r.I / r.sigma;
|
|
shell_isig_n[b] += 1;
|
|
}
|
|
}
|
|
|
|
double sum_e4 = 0.0;
|
|
int n_moment = 0;
|
|
for (int b = 0; b < n_shells; ++b) {
|
|
const auto& intensities = shell_I[b];
|
|
if (intensities.empty() || shell_isig_n[b] == 0
|
|
|| shell_isig_sum[b] / shell_isig_n[b] < min_shell_isig)
|
|
continue;
|
|
|
|
double sum = 0.0;
|
|
for (double I : intensities)
|
|
sum += I;
|
|
double mean = sum / intensities.size();
|
|
if (mean <= 0.0)
|
|
continue;
|
|
// Re-fit the mean over the reflections that pass the outlier cut, so the outlier does not
|
|
// inflate the very mean it is measured against.
|
|
sum = 0.0;
|
|
int n_kept = 0;
|
|
for (double I : intensities)
|
|
if (I / mean <= wilson_outlier_e2) {
|
|
sum += I;
|
|
++n_kept;
|
|
}
|
|
if (n_kept == 0)
|
|
continue;
|
|
mean = sum / n_kept;
|
|
if (mean <= 0.0)
|
|
continue;
|
|
|
|
for (double I : intensities) {
|
|
const double e2 = I / mean;
|
|
if (e2 > wilson_outlier_e2)
|
|
continue;
|
|
sum_e4 += e2 * e2;
|
|
++n_moment;
|
|
}
|
|
}
|
|
result.moment_reflections = n_moment;
|
|
if (n_moment > 0)
|
|
result.second_moment = sum_e4 / n_moment;
|
|
}
|
|
|
|
// Twin fraction from the second moment M = 2(1 - a + a^2): a = (1 - sqrt(2M-3))/2.
|
|
if (result.second_moment > 0.0) {
|
|
const double m = std::clamp(result.second_moment, 1.5, 2.0);
|
|
result.estimated_twin_fraction = (1.0 - std::sqrt(std::max(0.0, 2.0 * m - 3.0))) / 2.0;
|
|
}
|
|
// Either indicator dropping clearly below its untwinned value is suspicious - but only where a twin
|
|
// law can actually exist. In a holohedral Laue class (e.g. lysozyme's 422) no merohedral twinning is
|
|
// possible, so a low <|L|> is a statistical artefact (correlated near-neighbours) rather than a
|
|
// twin, and must not be flagged.
|
|
result.merohedral_twinning_possible = MerohedralTwinningPossible(space_group);
|
|
result.twinning_suspected = result.merohedral_twinning_possible &&
|
|
((result.l_test_pairs > 0 && result.mean_abs_l < 0.44) ||
|
|
(result.moment_reflections > 0 && result.second_moment < 1.85));
|
|
return result;
|
|
}
|
|
|
|
std::string TwinningAnalysisToText(const TwinningAnalysisResult& result) {
|
|
std::ostringstream os;
|
|
os << std::fixed << std::setprecision(3);
|
|
os << "Twinning analysis\n";
|
|
if (result.l_test_pairs > 0)
|
|
os << " L-test (Padilla-Yeates): <|L|> = " << result.mean_abs_l
|
|
<< ", <L^2> = " << result.mean_l_squared
|
|
<< " [untwinned 0.500 / 0.333, perfect twin 0.375 / 0.200; " << result.l_test_pairs
|
|
<< " pairs]\n";
|
|
if (result.moment_reflections > 0)
|
|
os << " Second moment <I^2>/<I>^2 = " << result.second_moment
|
|
<< " [untwinned 2.00, perfect twin 1.50]\n";
|
|
if (result.twinning_suspected)
|
|
os << " => Twinning suspected (estimated twin fraction ~"
|
|
<< result.estimated_twin_fraction << "). Statistics flag the presence of twinning, not\n"
|
|
<< " the twin law; confirm with a dedicated twin-law analysis.\n";
|
|
else if (!result.merohedral_twinning_possible)
|
|
os << " => No twinning: the Laue class is holohedral, so no merohedral twin law exists\n"
|
|
<< " (any <|L|> below 0.5 here is a statistical artefact, not twinning).\n";
|
|
else
|
|
os << " => No twinning indicated.\n";
|
|
return os.str();
|
|
}
|