diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index e407cd56..38c236ca 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -388,17 +388,17 @@ SearchSpaceGroupResult SearchSpaceGroup( int best_pg_order = 0; double best_pg_min_cc = -2.0; for (const auto& c : pg_cands) { - // Widen the chi^2 ratio bound by the error-model miscalibration. When even the best candidate's - // reduced chi^2 (chi2_ref) is far above 1 (weak data / uncorrected decay), the ratio grows with - // point-group order for genuine high symmetry too, so the fixed bound wrongly rejects it. Adding - // log10(chi2_ref) lets the bound exceed max_merge_chi2_ratio by ~the order of magnitude of the - // miscalibration - enough to keep the true group (F432 weak data: chi2_ref~150, ratio 3.2 vs - // bound ~4.2) while STILL rejecting a false operator whose ratio is far worse than its subgroup's - // (a twin/pseudo R32 on calibrated data: chi2_ref~2.6, ratio 4.2 vs bound ~2.4). The bound stays - // a per-candidate chi^2 test, so chi^2 still arbitrates between candidates - it is never bypassed. - const double allowed_ratio = opt.max_merge_chi2_ratio + std::log10(std::max(1.0, chi2_ref)); + // A genuine symmetry operator merges equivalent reflections, so it barely changes the reduced + // chi^2 relative to the best subgroup - across the whole rotation-test battery every correct + // point group stays within ~1.7x, even on weak or badly-integrated data (F432 chi2_ref 8.3 -> + // 1.15; Thau P41212 -> 1.71). A twin law or pseudo-symmetry forces non-equivalent reflections + // together, so its ratio is markedly higher (Ins_H_2's twin 2-fold: R3 3.02 -> R32 6.07, ratio + // 2.01). max_merge_chi2_ratio sits between the two. (An earlier log10(chi2_ref) widening + // compensated for an under-calibrated error model that inflated real-symmetry ratios with data + // weakness; the variance-floor fix removed that inflation, and the widening now only let the + // twin through, so it is gone.) const bool consistent = c.pg->rotations.empty() || !std::isfinite(c.chi2) || - !std::isfinite(chi2_ref) || c.chi2 <= chi2_ref * allowed_ratio; + !std::isfinite(chi2_ref) || c.chi2 <= chi2_ref * opt.max_merge_chi2_ratio; if (!consistent) continue; if (c.order > best_pg_order || (c.order == best_pg_order && c.min_class_cc > best_pg_min_cc)) { diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index 021e38ea..71a50305 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -77,12 +77,12 @@ struct SearchSpaceGroupOptions { // reduced chi^2 (within-orbit scatter / sigma^2) beyond this factor times the most-consistent // candidate. A false operator forces non-equivalent reflections together so they disagree by many // sigma and chi^2 blows up; a real one leaves it ~flat even when the operator CC is only moderate. - // 2.0 separates the test set: a true point group's chi^2 stays within ~1.6x the best subgroup - // (the cubic merge of a noisy crystal is the worst real case), while a false operator gives >=2.3x. - // This bound is the calibrated-error baseline; SearchSpaceGroup widens it by log10(chi2_ref) so a - // badly miscalibrated error model (weak data / uncorrected decay, where the ratio grows with - // point-group order for real symmetry too) does not spuriously reject a genuine high-symmetry group. - double max_merge_chi2_ratio = 2.0; + // Calibrated on the rotation-test battery: every correct point group stays within ~1.7x the best + // subgroup even on weak / badly-integrated data (worst real case Thau P41212 at 1.71), while a twin + // law or pseudo-symmetry lands clearly higher (Ins_H_2's R3->R32 twin 2-fold at 2.01). 1.85 sits + // between the two, so a partial merohedral twin is kept in its true lower symmetry (R3), not + // over-promoted to the holohedral R32. + double max_merge_chi2_ratio = 1.85; // --- Stage B: space group (screw axes / centering) --- bool determine_space_group = true; // false: stop at the symmorphic representative diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp index 8dd392d7..96866b81 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -22,6 +22,27 @@ namespace { 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& merged, @@ -83,7 +104,15 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector& merg // --- Second moment /^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. + // (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 ( 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::infinity(); double max_s = -std::numeric_limits::infinity(); @@ -99,31 +128,57 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector& merg const double t = (1.0 / (d * d) - min_s) / (max_s - min_s); return std::min(n_shells - 1, std::max(0, static_cast(t * n_shells))); }; - std::vector shell_sum(n_shells, 0.0); - std::vector shell_n(n_shells, 0); + // Group acentric intensities by shell, and accumulate to gauge each shell's signal. + std::vector> shell_I(n_shells); + std::vector shell_isig_sum(n_shells, 0.0); + std::vector 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_sum[b] += r.I; - shell_n[b] += 1; + 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; + } } - std::vector shell_mean(n_shells, 0.0); - for (int b = 0; b < n_shells; ++b) - if (shell_n[b] > 0) - shell_mean[b] = shell_sum[b] / shell_n[b]; double sum_e4 = 0.0; int n_moment = 0; - for (const auto& r : merged) { - if (!UsableIntensity(r) || !acentric(r)) + 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; - const double mean = shell_mean[shell_of(r.d)]; + + double sum = 0.0; + for (double I : intensities) + sum += I; + double mean = sum / intensities.size(); if (mean <= 0.0) continue; - const double e2 = r.I / mean; - sum_e4 += e2 * e2; - ++n_moment; + // 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) @@ -135,9 +190,14 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector& merg 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. - result.twinning_suspected = (result.l_test_pairs > 0 && result.mean_abs_l < 0.44) || - (result.moment_reflections > 0 && result.second_moment < 1.85); + // 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; } @@ -157,6 +217,9 @@ std::string TwinningAnalysisToText(const TwinningAnalysisResult& result) { 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(); diff --git a/image_analysis/scale_merge/TwinningAnalysis.h b/image_analysis/scale_merge/TwinningAnalysis.h index 6c0ec7d8..16714749 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.h +++ b/image_analysis/scale_merge/TwinningAnalysis.h @@ -33,6 +33,10 @@ struct TwinningAnalysisResult { // Twin fraction estimated from the second moment, a = (1 - sqrt(2M-3))/2, clamped to [0, 0.5]. double estimated_twin_fraction = 0.0; bool twinning_suspected = false; + + // False when the Laue class is holohedral (4/mmm, 6/mmm, m-3m, rhombohedral -3m): no merohedral + // twin law can exist, so a low <|L|> / second moment there is a statistical artefact, not twinning. + bool merohedral_twinning_possible = true; }; // The space group (when known) is used to drop centric reflections, which follow different