leonarski_fandClaude Opus 5.5
18263bccac
The refused group's representative can be given on the primitive lattice, which named primitive siblings (P222 beside an adopted I41). The refused group is taken with the adopted group's centring before its index-2 subgroups are listed. Report-only on every set the subset battery holds; a sibling-group model of a centred crystal is now matched. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37
527 lines
29 KiB
C++
527 lines
29 KiB
C++
// 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 <algorithm>
|
|
#include <array>
|
|
#include <string>
|
|
#include <cmath>
|
|
#include <map>
|
|
#include <tuple>
|
|
#include <numbers>
|
|
#include <vector>
|
|
|
|
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
|
|
#include "SyntheticMergedReflections.h"
|
|
|
|
namespace {
|
|
// A translational pseudo-symmetry u puts the factor |1 + exp(2 pi i h.u)|^2 on every intensity.
|
|
// Modelled here with a depth f so the weak class is suppressed rather than extinguished, which is
|
|
// what a real pair of copies with different orientations gives.
|
|
std::vector<MergedReflection> WithPseudoTranslation(std::vector<MergedReflection> merged,
|
|
const std::array<double, 3> &u, double f) {
|
|
for (auto &r : merged) {
|
|
const double phase = 2.0 * std::numbers::pi * (r.h * u[0] + r.k * u[1] + r.l * u[2]);
|
|
r.I = static_cast<float>(r.I * (1.0 + f * std::cos(phase)));
|
|
r.sigma = static_cast<float>(0.02 * std::fabs(r.I) + 1.0);
|
|
}
|
|
return merged;
|
|
}
|
|
|
|
std::vector<MergedReflection> Crystal(double twin_fraction) {
|
|
jfjoch_test::SyntheticMergeParams p;
|
|
p.true_space_group = "P 1 2 1";
|
|
p.twin_supergroup = "P 2 2 2";
|
|
p.twin_fraction = twin_fraction;
|
|
p.d_min_A = 3.0;
|
|
return jfjoch_test::GenerateSyntheticMerged(p);
|
|
}
|
|
}
|
|
|
|
// The L-test pairs a reflection with a partner two steps along an axis, and that choice is
|
|
// load-bearing for a reason the code did not state until now: an even step preserves the class of a
|
|
// HALF-INTEGER pseudo-translation, so the commonest tNCS leaves <|L|> alone by construction. A
|
|
// one-third pseudo-translation along the same axis does not, and it moves <|L|> by enough to change
|
|
// the verdict. This is the regression test for both halves of that.
|
|
TEST_CASE("L-test partner steps and a pseudo-translation", "[twinning][tncs]") {
|
|
const auto clean = Crystal(0.0);
|
|
const double l_clean = AnalyzeTwinning(clean, nullptr).mean_abs_l;
|
|
REQUIRE(l_clean == Catch::Approx(0.5).margin(0.02));
|
|
|
|
SECTION("a half-integer pseudo-translation does not move it") {
|
|
const auto modulated = WithPseudoTranslation(clean, {0.5, 0.0, 0.0}, 0.8);
|
|
const auto r = AnalyzeTwinning(modulated, nullptr);
|
|
CHECK(r.mean_abs_l == Catch::Approx(l_clean).margin(0.005));
|
|
}
|
|
|
|
SECTION("a one-third pseudo-translation moves it a long way up") {
|
|
const auto modulated = WithPseudoTranslation(clean, {1.0 / 3.0, 0.0, 0.0}, 0.8);
|
|
const auto r = AnalyzeTwinning(modulated, nullptr);
|
|
CHECK(r.mean_abs_l > l_clean + 0.03);
|
|
CHECK(r.mean_abs_l > 0.50); // into the "contradicts a twin" branch
|
|
}
|
|
|
|
SECTION("declaring the vector repairs it") {
|
|
const std::array<double, 3> u{1.0 / 3.0, 0.0, 0.0};
|
|
const std::vector<std::array<double, 3>> uv{u};
|
|
const auto modulated = WithPseudoTranslation(clean, u, 0.8);
|
|
const auto r = AnalyzeTwinning(modulated, nullptr, 20, &uv);
|
|
CHECK(r.l_test_tncs_step_restricted);
|
|
CHECK_FALSE(r.l_test_contaminated_by_tncs);
|
|
// The control is the SAME crystal without the pseudo-translation measured with the SAME
|
|
// restricted steps: <|L|> differs slightly between step directions (they are different
|
|
// distances in reciprocal space), so comparing against the unrestricted number would be
|
|
// measuring that instead of the repair.
|
|
const auto control = AnalyzeTwinning(clean, nullptr, 20, &uv);
|
|
REQUIRE(control.l_test_tncs_step_restricted);
|
|
CHECK(r.mean_abs_l == Catch::Approx(control.mean_abs_l).margin(0.005));
|
|
CHECK(r.l_test_pairs > 0.9 * AnalyzeTwinning(modulated, nullptr).l_test_pairs);
|
|
}
|
|
|
|
SECTION("a vector no step can preserve is declared unreadable, not repaired") {
|
|
const std::array<double, 3> u{1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0};
|
|
const std::vector<std::array<double, 3>> uv{u};
|
|
const auto modulated = WithPseudoTranslation(clean, u, 0.8);
|
|
const auto r = AnalyzeTwinning(modulated, nullptr, 20, &uv);
|
|
CHECK(r.l_test_contaminated_by_tncs);
|
|
CHECK_FALSE(r.l_test_tncs_step_restricted);
|
|
CHECK(r.l_test_pairs > 0); // still reported, just not read
|
|
}
|
|
|
|
SECTION("declaring a half-integer vector changes nothing") {
|
|
const std::array<double, 3> u{0.5, 0.0, 0.0};
|
|
const std::vector<std::array<double, 3>> uv{u};
|
|
const auto modulated = WithPseudoTranslation(clean, u, 0.8);
|
|
const auto with = AnalyzeTwinning(modulated, nullptr, 20, &uv);
|
|
const auto without = AnalyzeTwinning(modulated, nullptr);
|
|
CHECK(with.mean_abs_l == Catch::Approx(without.mean_abs_l).margin(0.002));
|
|
CHECK_FALSE(with.l_test_contaminated_by_tncs);
|
|
}
|
|
}
|
|
|
|
// The defect the repair exists for: a twinned crystal that also carries a one-third
|
|
// pseudo-translation loses its twin call entirely, because the pseudo-symmetry pushes <|L|> up out
|
|
// of the twinned range AND pushes the second moment up out of it at the same time.
|
|
TEST_CASE("A pseudo-translation can hide a twin, and declaring it restores the call",
|
|
"[twinning][tncs]") {
|
|
const auto twin = Crystal(0.5);
|
|
const auto baseline = AnalyzeTwinning(twin, nullptr);
|
|
REQUIRE(baseline.mean_abs_l < 0.44); // an unambiguous twin when nothing masks it
|
|
REQUIRE(baseline.twinning_suspected);
|
|
|
|
const std::array<double, 3> u{1.0 / 3.0, 0.0, 0.0};
|
|
const std::vector<std::array<double, 3>> uv{u};
|
|
const auto masked = WithPseudoTranslation(twin, u, 0.8);
|
|
|
|
const auto undeclared = AnalyzeTwinning(masked, nullptr);
|
|
CHECK(undeclared.mean_abs_l > baseline.mean_abs_l + 0.05);
|
|
CHECK_FALSE(undeclared.twinning_suspected); // the twin call is lost
|
|
|
|
const auto declared = AnalyzeTwinning(masked, nullptr, 20, &uv);
|
|
// Against the same twin without the pseudo-translation, measured with the same restricted steps.
|
|
const auto control = AnalyzeTwinning(twin, nullptr, 20, &uv);
|
|
CHECK(declared.mean_abs_l == Catch::Approx(control.mean_abs_l).margin(0.005));
|
|
CHECK(declared.mean_abs_l < 0.44);
|
|
CHECK(declared.twinning_suspected); // and restored
|
|
}
|
|
|
|
namespace {
|
|
// The merge of a synthetic crystal in `group`: one row per reflection of that group's asymmetric
|
|
// unit, the average of the P1 intensities of its orbit - which is what merging under the group
|
|
// does, whether or not its operators are real.
|
|
std::vector<MergedReflection> MergedIn(const std::vector<MergedReflection> &p1, const char *group) {
|
|
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(group);
|
|
const gemmi::GroupOps gops = sg->operations();
|
|
const gemmi::ReciprocalAsu asu(sg);
|
|
std::map<std::array<int, 3>, std::pair<MergedReflection, int>> sum;
|
|
for (const auto &r : p1) {
|
|
const auto key = asu.to_asu(gemmi::Op::Miller{{r.h, r.k, r.l}}, gops).first;
|
|
auto [it, fresh] = sum.try_emplace({key[0], key[1], key[2]}, r, 0);
|
|
if (!fresh)
|
|
it->second.first.I += r.I;
|
|
it->second.second += 1;
|
|
}
|
|
std::vector<MergedReflection> out;
|
|
for (auto &[key, v] : sum) {
|
|
MergedReflection r = v.first;
|
|
r.h = key[0];
|
|
r.k = key[1];
|
|
r.l = key[2];
|
|
r.I /= static_cast<float>(v.second);
|
|
out.push_back(r);
|
|
}
|
|
return out;
|
|
}
|
|
|
|
std::vector<MergedReflection> Tetragonal(const char *true_group, double twin_fraction) {
|
|
jfjoch_test::SyntheticMergeParams p;
|
|
p.true_space_group = true_group;
|
|
p.twin_supergroup = "P 4 2 2";
|
|
p.twin_fraction = twin_fraction;
|
|
p.d_min_A = 2.5;
|
|
return jfjoch_test::GenerateSyntheticMerged(p);
|
|
}
|
|
}
|
|
|
|
// The L-test compares two intensities as if they had the same expected value, and two index steps are
|
|
// not the same resolution: on a small cell with a steep fall-off the raw pair differs systematically,
|
|
// which reads as an untwinned crystal being "more untwinned" than 0.5. The shell normalisation removes
|
|
// that, and the fraction quoted is the L-test's own.
|
|
TEST_CASE("L-test on shell-normalised intensities", "[twinning]") {
|
|
jfjoch_test::SyntheticMergeParams p;
|
|
p.true_space_group = "P 1";
|
|
p.twin_supergroup = "P 1";
|
|
p.wilson_b_A2 = 60.0;
|
|
p.d_min_A = 2.5;
|
|
const auto r = AnalyzeTwinning(jfjoch_test::GenerateSyntheticMerged(p), nullptr);
|
|
INFO("<|L|> " << r.mean_abs_l);
|
|
CHECK(r.mean_abs_l == Catch::Approx(0.5).margin(0.02));
|
|
CHECK_FALSE(r.twinning_suspected);
|
|
|
|
const auto twin = AnalyzeTwinning(Crystal(0.2), nullptr);
|
|
CHECK(twin.twin_fraction_source == TwinFractionSource::LTest);
|
|
CHECK(twin.estimated_twin_fraction == Catch::Approx(0.2).margin(0.05));
|
|
}
|
|
|
|
// A P1 merge of a centred lattice holds the reflections the centring extinguishes. Told the centring,
|
|
// the test leaves them out and reads the crystal; not told, it pairs present with absent ones.
|
|
TEST_CASE("L-test on a P1 merge of a centred lattice", "[twinning]") {
|
|
jfjoch_test::SyntheticMergeParams p;
|
|
p.true_space_group = "R 3 :H";
|
|
p.twin_supergroup = "R 3 :H";
|
|
p.d_min_A = 3.0;
|
|
const auto merged = jfjoch_test::GenerateSyntheticMerged(p);
|
|
const auto told = AnalyzeTwinning(merged, nullptr, 20, nullptr, 'R');
|
|
CHECK(told.mean_abs_l == Catch::Approx(0.5).margin(0.02));
|
|
const auto not_told = AnalyzeTwinning(merged, nullptr);
|
|
CHECK(std::fabs(not_told.mean_abs_l - 0.5) > 0.05);
|
|
}
|
|
|
|
// A genuine holohedral crystal: no twin law exists, and its operators read as real.
|
|
TEST_CASE("No twin is called in a genuine holohedral Laue class", "[twinning]") {
|
|
const auto merged = MergedIn(Tetragonal("P 4 2 2", 0.0), "P 4 2 2");
|
|
const gemmi::UnitCell tetragonal(47, 47, 63, 90, 90, 90);
|
|
const auto r = AnalyzeTwinning(merged, gemmi::find_spacegroup_by_name("P 4 2 2"), 20, nullptr, 'P', &tetragonal);
|
|
CHECK_FALSE(r.merohedral_twinning_possible);
|
|
CHECK(r.mean_abs_l == Catch::Approx(0.5).margin(0.03));
|
|
CHECK_FALSE(r.twinning_suspected);
|
|
CHECK_FALSE(r.adopted_operators_suspect);
|
|
}
|
|
|
|
// Merging under a false operator averages I(h) with I(Th), and (I(h) + I(Th))/2 has the perfect-twin
|
|
// distribution whatever the twin fraction - so a subgroup crystal merged in its lattice's holohedry
|
|
// reads <|L|> ~0.375 at every fraction, and the holohedral class must say so rather than "no twin
|
|
// law exists".
|
|
TEST_CASE("A holohedral merge under a false operator is called suspect", "[twinning]") {
|
|
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 4 2 2");
|
|
const gemmi::UnitCell tetragonal(47, 47, 63, 90, 90, 90);
|
|
for (double alpha : {0.0, 0.2, 0.5}) {
|
|
const auto r = AnalyzeTwinning(MergedIn(Tetragonal("P 4", alpha), "P 4 2 2"), sg, 20, nullptr, 'P', &tetragonal);
|
|
CAPTURE(alpha);
|
|
CHECK(r.mean_abs_l == Catch::Approx(0.375).margin(0.02));
|
|
CHECK(r.adopted_operators_suspect);
|
|
CHECK(r.twin_fraction_source == TwinFractionSource::None); // lost in the merge
|
|
CHECK(TwinningVerdictLine(r).rfind("SYMMETRY SUSPECT", 0) == 0);
|
|
}
|
|
}
|
|
|
|
// Whether a twin law can exist is read off the lattice, not off a list of Laue classes: a group whose
|
|
// lattice has a rotation the group lacks admits one, whatever the class is called.
|
|
TEST_CASE("A twin law is possible exactly where the lattice has a rotation the group lacks", "[twinning]") {
|
|
const auto possible = [](const char *group, const gemmi::UnitCell &cell) {
|
|
return TwinLawPossible(cell, gemmi::get_spacegroup_by_name(group));
|
|
};
|
|
// Merohedral: a tetragonal lattice under a group of lower Laue class than its holohedry.
|
|
CHECK(possible("P 4", gemmi::UnitCell(47, 47, 63, 90, 90, 90)));
|
|
CHECK_FALSE(possible("P 4 2 2", gemmi::UnitCell(47, 47, 63, 90, 90, 90)));
|
|
// Orthorhombic: no law on a general metric - mmm is the lattice's own holohedry - and a law where
|
|
// the C-centred cell is pseudo-hexagonal (b = a sqrt(3)).
|
|
CHECK_FALSE(possible("C 2 2 21", gemmi::UnitCell(56, 157, 143, 90, 90, 90)));
|
|
CHECK(possible("C 2 2 21", gemmi::UnitCell(50, 50 * std::sqrt(3.0), 143, 90, 90, 90)));
|
|
// Monoclinic: none on a general metric, one where beta is 90 (pseudo-merohedral).
|
|
CHECK_FALSE(possible("P 1 21 1", gemmi::UnitCell(75, 94, 99, 90, 93.1, 90)));
|
|
CHECK(possible("P 1 21 1", gemmi::UnitCell(55, 108, 55.3, 90, 90.0, 90)));
|
|
// ... and within the obliquity it is read at, a degree off still counts.
|
|
CHECK(possible("P 1 21 1", gemmi::UnitCell(54, 142, 84, 90, 91.0, 90)));
|
|
}
|
|
|
|
// A twin-like <|L|> in a group the lattice admits no twin law for is not a twin in that group: it is
|
|
// reported as such, and not as "needs a twin law". Where the lattice does admit one, it is a twin.
|
|
TEST_CASE("Twin-like statistics where no twin law is possible", "[twinning]") {
|
|
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 2 2 2");
|
|
const auto merged = MergedIn(Tetragonal("P 2 2 2", 0.15), "P 2 2 2");
|
|
const gemmi::UnitCell orthorhombic(47, 52, 63, 90, 90, 90);
|
|
const gemmi::UnitCell tetragonal(47, 47, 63, 90, 90, 90);
|
|
|
|
const auto none = AnalyzeTwinning(merged, sg, 20, nullptr, 'P', &orthorhombic);
|
|
INFO("<|L|> " << none.mean_abs_l);
|
|
REQUIRE(none.mean_abs_l > 0.42);
|
|
REQUIRE(none.mean_abs_l < 0.44);
|
|
CHECK_FALSE(none.merohedral_twinning_possible);
|
|
CHECK_FALSE(none.twinning_suspected);
|
|
CHECK_FALSE(none.adopted_operators_suspect);
|
|
CHECK(none.twin_like_without_twin_law);
|
|
CHECK(TwinningVerdictLine(none).rfind("TWIN-LIKE, NO TWIN LAW", 0) == 0);
|
|
|
|
const auto law = AnalyzeTwinning(merged, sg, 20, nullptr, 'P', &tetragonal);
|
|
CHECK(law.merohedral_twinning_possible);
|
|
CHECK(law.twinning_suspected);
|
|
CHECK_FALSE(law.twin_like_without_twin_law);
|
|
|
|
// Without a cell nothing can be ruled out.
|
|
CHECK(AnalyzeTwinning(merged, sg).merohedral_twinning_possible);
|
|
}
|
|
|
|
namespace {
|
|
// A P1 merge whose centric reflections are centric: SyntheticMergedReflections draws every
|
|
// intensity from the acentric distribution, which is all the searches it was written for need, but
|
|
// the zone test reads exactly the difference. E^2 is chi^2 with one degree of freedom for a
|
|
// reflection centric in the true group (a squared Box-Muller normal), exponential otherwise, keyed
|
|
// on the asymmetric-unit index so the true group's symmetry is exact. `anisotropy_b` is an extra
|
|
// Debye-Waller B along c* alone (A^2) - the same in both twin domains, a twin law being a lattice
|
|
// operation - and `scale_jitter` a log-normal factor of that width on every observed intensity,
|
|
// a nuisance with no direction that no normalisation removes.
|
|
// `twin_super` names the group whose extra operator is the twin law.
|
|
std::vector<MergedReflection> WilsonMerge(const char *true_group, const char *twin_super,
|
|
const gemmi::UnitCell &cell, double twin_fraction,
|
|
double anisotropy_b = 0.0, double scale_jitter = 0.0) {
|
|
const gemmi::SpaceGroup &sub = gemmi::get_spacegroup_by_name(true_group);
|
|
const gemmi::SpaceGroup &super = gemmi::get_spacegroup_by_name(twin_super);
|
|
const gemmi::Op twin = jfjoch_test::TwinLaw(sub, super);
|
|
const gemmi::GroupOps gops = sub.operations();
|
|
const gemmi::ReciprocalAsu rasu(&sub);
|
|
const int hmax = static_cast<int>(cell.a / 2.5) + 1, kmax = static_cast<int>(cell.b / 2.5) + 1,
|
|
lmax = static_cast<int>(cell.c / 2.5) + 1;
|
|
auto true_intensity = [&](const gemmi::Op::Miller &hkl) {
|
|
const auto asu = rasu.to_asu(hkl, gops).first;
|
|
const double u1 = jfjoch_test::detail::UniformFromHkl(asu);
|
|
const double u2 = jfjoch_test::detail::UniformFromHkl({{asu[0] + 300, asu[1] + 300, asu[2] + 300}});
|
|
const double c = std::cos(2.0 * std::numbers::pi * u2);
|
|
const double e2 = gops.is_reflection_centric(hkl) ? -2.0 * std::log(u1) * c * c : -std::log(u1);
|
|
const double s_c = hkl[2] / cell.c;
|
|
return 1000.0 * e2 * std::exp(-20.0 / (2.0 * std::pow(cell.calculate_d(hkl), 2)))
|
|
* std::exp(-anisotropy_b * s_c * s_c / 2.0);
|
|
};
|
|
std::vector<MergedReflection> out;
|
|
for (int h = -hmax; h <= hmax; ++h)
|
|
for (int k = -kmax; k <= kmax; ++k)
|
|
for (int l = -lmax; l <= lmax; ++l) {
|
|
if (std::make_tuple(h, k, l) <= std::make_tuple(-h, -k, -l))
|
|
continue;
|
|
const gemmi::Op::Miller hkl{{h, k, l}};
|
|
const double d = cell.calculate_d(hkl);
|
|
if (d < 2.5)
|
|
continue;
|
|
MergedReflection r;
|
|
r.h = h;
|
|
r.k = k;
|
|
r.l = l;
|
|
r.d = static_cast<float>(d);
|
|
const double u3 = jfjoch_test::detail::UniformFromHkl({{h + 700, k + 700, l + 700}});
|
|
const double u4 = jfjoch_test::detail::UniformFromHkl({{h + 900, k + 900, l + 900}});
|
|
const double jitter = std::exp(scale_jitter * std::sqrt(-2.0 * std::log(u3))
|
|
* std::cos(2.0 * std::numbers::pi * u4));
|
|
r.I = static_cast<float>(((1.0 - twin_fraction) * true_intensity(hkl)
|
|
+ twin_fraction * true_intensity(twin.apply_to_hkl(hkl))) * jitter);
|
|
r.sigma = static_cast<float>(0.01 * r.I + 1.0);
|
|
out.push_back(r);
|
|
}
|
|
return out;
|
|
}
|
|
|
|
std::vector<MergedReflection> WilsonTetragonal(const char *true_group, double twin_fraction,
|
|
double anisotropy_b = 0.0, double scale_jitter = 0.0) {
|
|
return WilsonMerge(true_group, "P 4 2 2", gemmi::UnitCell(47, 47, 63, 90, 90, 90), twin_fraction,
|
|
anisotropy_b, scale_jitter);
|
|
}
|
|
}
|
|
|
|
// Reflections centric in the adopted group but acentric in a subgroup are their own twin mates under
|
|
// the operators the group adds. A perfect twin makes every operator statistic read "real", but these
|
|
// stay acentric under it and are centric only where the added operators are real.
|
|
TEST_CASE("Twin-immune zones tell a perfect twin from real symmetry", "[twinning]") {
|
|
const gemmi::SpaceGroup &p422 = gemmi::get_spacegroup_by_name("P 4 2 2");
|
|
const gemmi::UnitCell tetragonal(47, 47, 63, 90, 90, 90);
|
|
|
|
const auto real = AnalyzeTwinImmuneZones(WilsonTetragonal("P 4 2 2", 0.0), tetragonal, p422);
|
|
REQUIRE(real.zones.size() == 3); // over 4, over 222, over the diagonal 222
|
|
for (const auto &z : real.zones) {
|
|
CAPTURE(z.operators);
|
|
CHECK(z.mean_abs_e2_minus_1 == Catch::Approx(0.968).margin(0.08));
|
|
CHECK(z.evidence_nats > 0.0);
|
|
}
|
|
CHECK(real.control.mean_abs_e2_minus_1 == Catch::Approx(0.736).margin(0.05));
|
|
|
|
// A perfect P4 twin: whichever subgroup is taken as the parent, the operators added over it are
|
|
// not all real, and the zone reads acentric.
|
|
const auto twin = AnalyzeTwinImmuneZones(WilsonTetragonal("P 4", 0.5), tetragonal, p422);
|
|
REQUIRE(twin.zones.size() == 3);
|
|
for (const auto &z : twin.zones) {
|
|
CAPTURE(z.operators);
|
|
CHECK(z.mean_abs_e2_minus_1 == Catch::Approx(0.736).margin(0.08));
|
|
CHECK(z.evidence_nats < 0.0);
|
|
CHECK(z.calibrated_evidence_nats < 0.0);
|
|
}
|
|
// Wilson data normalise cleanly: the control reads the acentric expectation and the
|
|
// calibration takes nothing off.
|
|
CHECK(real.control_excess_per_reflection < 0.01);
|
|
CHECK(twin.control_excess_per_reflection < 0.01);
|
|
}
|
|
|
|
// The zone is read absolutely, so it is only as good as the normalisation. A strong anisotropy
|
|
// leaves every class of an isotropically normalised merge reading centric, the acentric control
|
|
// included, and a zone - a plane - more so than the control's sphere; it is fitted and taken out
|
|
// before anything is read. What no normalisation removes, the control certifies: its excess over the
|
|
// acentric expectation is taken off every zone, so a partial twin's zone reads acentric and a real
|
|
// operator's centric, as they do on clean data.
|
|
TEST_CASE("Twin-immune zones are normalised for anisotropy and calibrated by the acentric control",
|
|
"[twinning]") {
|
|
const gemmi::SpaceGroup &p422 = gemmi::get_spacegroup_by_name("P 4 2 2");
|
|
const gemmi::SpaceGroup &p4 = gemmi::get_spacegroup_by_name("P 4");
|
|
const gemmi::UnitCell tetragonal(47, 47, 63, 90, 90, 90);
|
|
|
|
SECTION("a 60 A^2 anisotropy along c* is taken out") {
|
|
const auto twin = AnalyzeTwinImmuneZone(WilsonTetragonal("P 4", 0.2, 60.0), tetragonal, p422, p4);
|
|
REQUIRE(twin.zones.size() == 1);
|
|
CHECK(twin.anisotropy_delta_b_A2 == Catch::Approx(60.0).margin(15.0));
|
|
CHECK(twin.control.mean_abs_e2_minus_1 < 0.76); // a partial twin's control: at or below acentric
|
|
CHECK(twin.zones[0].calibrated_evidence_nats < -20.0);
|
|
|
|
const auto real = AnalyzeTwinImmuneZone(WilsonTetragonal("P 4 2 2", 0.0, 60.0), tetragonal, p422, p4);
|
|
REQUIRE(real.zones.size() == 1);
|
|
CHECK(real.control.mean_abs_e2_minus_1 == Catch::Approx(0.736).margin(0.05));
|
|
CHECK(real.zones[0].calibrated_evidence_nats > 20.0);
|
|
}
|
|
|
|
SECTION("a scale jitter no normalisation removes is calibrated off") {
|
|
const auto real = AnalyzeTwinImmuneZone(WilsonTetragonal("P 4 2 2", 0.0, 0.0, 0.6), tetragonal, p422, p4);
|
|
REQUIRE(real.zones.size() == 1);
|
|
CHECK(real.control.mean_abs_e2_minus_1 > 0.80); // inflated: the normalisation is wrong
|
|
CHECK(real.control_excess_per_reflection > 0.02);
|
|
CHECK(real.zones[0].calibrated_evidence_nats > 20.0);
|
|
|
|
// The twin's control is deflated by the twinning as much as the jitter inflates it, so
|
|
// nothing is taken off - and its zone, acentric, still reads so.
|
|
const auto twin = AnalyzeTwinImmuneZone(WilsonTetragonal("P 4", 0.2, 0.0, 0.6), tetragonal, p422, p4);
|
|
REQUIRE(twin.zones.size() == 1);
|
|
CHECK(twin.zones[0].calibrated_evidence_nats < -20.0);
|
|
}
|
|
}
|
|
|
|
// A twin by a law OTHER than the operators a zone is read for averages the zone's reflections with
|
|
// unrelated ones like every other reflection, so a genuine 321 crystal twinned by a 622 law reads its
|
|
// 2-folds' zone acentric against the untwinned centric density. Read at the twin fraction, it reads
|
|
// centric again - while a crystal twinned by the zone's own operators, which cannot reach the zone,
|
|
// stays acentric.
|
|
TEST_CASE("Twin-immune zones are read at the twin fraction of the other twin laws", "[twinning]") {
|
|
const gemmi::SpaceGroup &p321 = gemmi::get_spacegroup_by_name("P 3 2 1");
|
|
const gemmi::SpaceGroup &p3 = gemmi::get_spacegroup_by_name("P 3");
|
|
const gemmi::UnitCell hexagonal(51, 51, 71, 90, 90, 120);
|
|
|
|
const auto genuine = WilsonMerge("P 3 2 1", "P 6 2 2", hexagonal, 0.1);
|
|
const auto untwinned_read = AnalyzeTwinImmuneZone(genuine, hexagonal, p321, p3);
|
|
REQUIRE(untwinned_read.zones.size() == 1);
|
|
CHECK(untwinned_read.zones[0].calibrated_evidence_nats < -20.0);
|
|
const auto twinned_read = AnalyzeTwinImmuneZone(genuine, hexagonal, p321, p3, nullptr, 1, 0.1);
|
|
CHECK(twinned_read.zones[0].calibrated_evidence_nats > 20.0);
|
|
|
|
// Twinned by the 2-folds themselves, with no other law to read the zone at.
|
|
const auto twin = WilsonMerge("P 3", "P 3 2 1", hexagonal, 0.1);
|
|
CHECK(AnalyzeTwinImmuneZone(twin, hexagonal, p321, p3).zones[0].calibrated_evidence_nats < -20.0);
|
|
}
|
|
|
|
// Averaging each intensity over its orbit under a group leaves the L-test where it was when the group
|
|
// is the crystal's, and narrows it towards the perfect-twin 0.375 when an operator of the group is
|
|
// not - read on the same pairs. A perfect twin of the subgroup already reads 0.375 unmerged, so there
|
|
// the averaging has nothing left to show, which is the limit no intensity statistic passes.
|
|
TEST_CASE("L-test under a merge tells a false operator from a real one", "[twinning]") {
|
|
std::vector<gemmi::Op> rotations;
|
|
for (const auto &op : gemmi::get_spacegroup_by_name("P 4 2 2").operations().sym_ops)
|
|
if (op.rot != gemmi::Op::identity().rot)
|
|
rotations.push_back(op);
|
|
|
|
const auto genuine = AnalyzeLTestUnderMerge(WilsonTetragonal("P 4 2 2", 0.0), rotations, 'P');
|
|
REQUIRE(genuine.pairs > 1000);
|
|
CHECK(genuine.mean_abs_l == Catch::Approx(0.5).margin(0.03));
|
|
CHECK(genuine.mean_abs_l_merged == Catch::Approx(genuine.mean_abs_l).margin(0.01));
|
|
|
|
const auto false_operator = AnalyzeLTestUnderMerge(WilsonTetragonal("P 4", 0.0), rotations, 'P');
|
|
CHECK(false_operator.mean_abs_l == Catch::Approx(0.5).margin(0.03));
|
|
CHECK(false_operator.mean_abs_l_merged == Catch::Approx(0.375).margin(0.02));
|
|
|
|
const auto perfect_twin = AnalyzeLTestUnderMerge(WilsonTetragonal("P 4", 0.5), rotations, 'P');
|
|
CHECK(perfect_twin.mean_abs_l == Catch::Approx(0.375).margin(0.02));
|
|
CHECK(perfect_twin.mean_abs_l_merged == Catch::Approx(perfect_twin.mean_abs_l).margin(0.01));
|
|
}
|
|
|
|
// The subgroups a twin-law merge can be made in: index 2, translations kept, and a two-fold of the
|
|
// group outside them to be the law. A four-fold over its two-fold offers no two-fold, so no law.
|
|
TEST_CASE("Index-2 subgroups with a two-fold twin law", "[twinning]") {
|
|
const auto names = [](const char *group) {
|
|
std::vector<std::string> out;
|
|
for (const auto &t : IndexTwoTwinSubgroups(gemmi::get_spacegroup_by_name(group))) {
|
|
CHECK_FALSE(t.laws.empty());
|
|
out.push_back(t.subgroup.xhm());
|
|
}
|
|
std::sort(out.begin(), out.end());
|
|
return out;
|
|
};
|
|
CHECK(names("P 2 2 2") == std::vector<std::string>{"P 1 1 2", "P 1 2 1", "P 2 1 1"});
|
|
CHECK(names("P 21 21 21") == std::vector<std::string>{"P 1 1 21", "P 1 21 1", "P 21 1 1"});
|
|
CHECK(names("P 4").empty());
|
|
const auto p422 = names("P 4 2 2");
|
|
CHECK(std::find(p422.begin(), p422.end(), "P 4") != p422.end());
|
|
CHECK(std::find(p422.begin(), p422.end(), "P 2 2 2") != p422.end());
|
|
CHECK(names("C 2 2 21").size() == 3);
|
|
}
|
|
|
|
// The perfect-twin hypotheses a refused promotion leaves beside the adopted group: the other index-2
|
|
// subgroups of the refused group, each with the two-folds of its coset as laws. A P3x21 adopted under a
|
|
// refused 622 leaves P6 (law k,h,-l among others) and P312; a group that is not an index-2 subgroup of
|
|
// the refused one leaves none.
|
|
TEST_CASE("Sibling twin hypotheses of a refused promotion", "[twinning]") {
|
|
const auto siblings = [](const char *adopted, const char *refused) {
|
|
return SiblingTwinSubgroups(gemmi::get_spacegroup_by_name(adopted), gemmi::get_spacegroup_by_name(refused));
|
|
};
|
|
const auto names = [](const std::vector<TwinSubgroup> &v) {
|
|
std::vector<std::string> out;
|
|
for (const auto &t : v) {
|
|
CHECK_FALSE(t.laws.empty());
|
|
out.push_back(t.subgroup.xhm());
|
|
}
|
|
std::sort(out.begin(), out.end());
|
|
return out;
|
|
};
|
|
const auto hexagonal = siblings("P 31 2 1", "P 6 2 2");
|
|
CHECK(names(hexagonal) == std::vector<std::string>{"P 3 1 2", "P 6"});
|
|
for (const auto &t : hexagonal) {
|
|
std::vector<std::string> laws;
|
|
for (const auto &l : t.laws)
|
|
laws.push_back(l.as_hkl().triplet('h'));
|
|
CAPTURE(t.subgroup.xhm(), laws);
|
|
if (t.subgroup.xhm() == "P 6")
|
|
CHECK(std::find(laws.begin(), laws.end(), "k,h,-l") != laws.end());
|
|
else
|
|
CHECK(std::find(laws.begin(), laws.end(), "-h,-k,l") != laws.end());
|
|
}
|
|
CHECK(names(siblings("P 1 21 1", "P 21 21 21")) == std::vector<std::string>{"P 1 1 21", "P 21 1 1"});
|
|
// A centred crystal's siblings are centred, whatever centring the refused group was given with.
|
|
const auto centred = names(siblings("I 41", "P 4 2 2"));
|
|
CHECK(std::find(centred.begin(), centred.end(), "I 2 2 2") != centred.end());
|
|
CHECK(siblings("P 4", "P 6 2 2").empty());
|
|
CHECK(siblings("P 6 2 2", "P 6 2 2").empty());
|
|
}
|
|
|
|
// Yeates' H on the pairs a law relates, in the group the law is not part of, reads the twin fraction.
|
|
TEST_CASE("Twin fraction of a law from its pairs' disagreement", "[twinning]") {
|
|
const gemmi::SpaceGroup &p4 = gemmi::get_spacegroup_by_name("P 4");
|
|
const gemmi::Op law = jfjoch_test::TwinLaw(p4, gemmi::get_spacegroup_by_name("P 4 2 2"));
|
|
for (double alpha : {0.0, 0.2, 0.4}) {
|
|
CAPTURE(alpha);
|
|
const auto f = EstimateTwinLawFraction(MergedIn(Tetragonal("P 4", alpha), "P 4"), p4, law);
|
|
CHECK(f.pairs > 1000);
|
|
CHECK(f.fraction == Catch::Approx(alpha).margin(0.03));
|
|
}
|
|
}
|