Space-group search: judge a screw axis against its own axial row
A reflection the group predicts absent counted as a violation when
I/sigma > 3 AND E^2 = I/<I>(shell) > 0.3. Neither half survives contact
with real data:
* merged sigma is floored at b|I|, so merged I/sigma saturates at ISa
for nearly every reflection - the I/sigma half is an on/off switch
keyed on ISa vs 3, not a per-reflection test. On one crystal the
absent class read <I/s> 4.10 against 3.73 for the present class while
being genuinely extinct;
* <I>(shell) decays with resolution while a systematically-absent
reflection keeps a small NON-decaying residual (background / profile
leakage), so absent reflections drift over an absolute E^2 cut at high
resolution. That cost a tetragonal 42_12 crystal its 4_1: 18 of its 47
absent 00l crossed the cut, all beyond 3.7 A, at absolute intensities
identical to the low-resolution ones correctly judged absent, while
their l=4n row-mates sat 20-60x higher at the same resolution.
A screw extinguishes only the reflections that lie ON its axis, so the
fair yardstick is the rest of that same row. The threshold is now
0.3 * max(1, median E^2 of the reflection's own row), the row being the
gcd-reduced reciprocal-space direction and the control class the same-row
reflections the group predicts present. Floored at 1, so it only ever
relaxes: a screw can be recovered by it, never lost.
Per row, not pooled. A 4_1 along c and a 2_1 along a are separate
conditions with separate controls; pooling let the weak a/b rows (median
E^2 ~0.5) set the threshold for a strong c row (8.4) and the rescue never
fired.
The candidate table now reports the screw evidence (median E^2 of the
absent class and of its rows) - the <I/s> columns are the centering
evidence and say nothing about screws, for the sigma-floor reason above.
Rotation battery, 33 crystals: 31 decisions bit-identical, the 42_12
crystal recovers its 4_1 (0 violations, row E^2 8.4 vs absent 0.12), and
one crystal with a long axis and heavy 00l overlap moves to a 4_1 group at
exactly 10.0% violations - marginal, and its sister crystal of the same
form sits at 13.3% and does not move. Real screws now span 0-9.3%
violations, so max_absent_violation_fraction cannot be tightened below
0.10 without risking a genuine one.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
@@ -190,4 +190,48 @@ TEST_CASE("SearchSpaceGroup finds a screw axis despite under-estimated sigmas on
|
||||
REQUIRE(result.best_space_group.has_value());
|
||||
CHECK(result.best_space_group->short_name() == "P2");
|
||||
}
|
||||
}
|
||||
|
||||
// Regression: the E^2 gate above compares a reflection to the mean of its RESOLUTION SHELL, which
|
||||
// falls off with resolution, while a systematically-absent reflection keeps a small non-decaying
|
||||
// residual (background / profile leakage). On a crystal whose axial rows are much stronger than an
|
||||
// average reflection, that turns the high-resolution residuals into screw-axis violations and the
|
||||
// screw is lost, although the reflections beside them in the same row are tens of times stronger.
|
||||
// A tetragonal 42_12 case failed exactly this way (18 of 47 absent 00l over the cut, all beyond
|
||||
// 3.7 A, at 1-2% of the l=4n reflections next to them). The threshold is therefore taken relative to
|
||||
// the axial row the screw constrains, not to the shell.
|
||||
TEST_CASE("SearchSpaceGroup finds a screw axis whose absent class is weak only within its own row") {
|
||||
const gemmi::SpaceGroup& sg = gemmi::get_spacegroup_by_name("P 43 21 2");
|
||||
auto merged = GenerateMergedReflectionsForSpaceGroup(sg, 12);
|
||||
|
||||
// Axial rows 40x stronger than a general reflection, and an absent class carrying ~2% of its own
|
||||
// row - but half of a general reflection, so a threshold set against the shell calls every one of
|
||||
// them a violation while a threshold set against the row calls none.
|
||||
const gemmi::GroupOps gops = sg.operations();
|
||||
int absent_on_axis = 0;
|
||||
for (auto& r : merged) {
|
||||
const gemmi::Op::Miller hkl{{r.h, r.k, r.l}};
|
||||
if (gops.epsilon_factor_without_centering(hkl) <= 1)
|
||||
continue;
|
||||
if (gops.is_systematically_absent(hkl)) {
|
||||
r.I = 300.0f;
|
||||
r.sigma = 1.0f;
|
||||
++absent_on_axis;
|
||||
} else {
|
||||
r.I *= 40.0f;
|
||||
}
|
||||
}
|
||||
REQUIRE(absent_on_axis >= 8);
|
||||
|
||||
SearchSpaceGroupOptions opt;
|
||||
opt.merge_friedel = true;
|
||||
|
||||
const auto result = SearchSpaceGroup(merged, opt);
|
||||
INFO(SearchSpaceGroupResultToText(result));
|
||||
REQUIRE(result.best_space_group.has_value());
|
||||
// P4_1 2_1 2 and P4_3 2_1 2 are enantiomorphs and indistinguishable from intensities.
|
||||
std::vector<std::string> accepted{result.best_space_group->short_name()};
|
||||
for (const auto& alt : result.alternatives)
|
||||
accepted.push_back(alt.short_name());
|
||||
CHECK(std::find(accepted.begin(), accepted.end(), "P43212") != accepted.end());
|
||||
}
|
||||
Reference in New Issue
Block a user