Merge branch 'rc186-twin-x' into rc186-merge
This commit is contained in:
leonarski_f committed
commit
2e5a8812ee
7 files changed
+118
-23
No files matched your search
@@ -62,7 +62,7 @@ angle < 0.44$) and the averaging closed a third or more of the gap from the unme
|
||||
|
||||
### 13.2 Twinning check, and translational pseudo-symmetry
|
||||
|
||||
A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$ — 0.500 and 0.333 untwinned, 0.375 and 0.200 for a perfect twin) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ (2.0 for untwinned acentric data, 1.5 for a perfect twin) are written to the merged mmCIF as a twinning diagnostic. Both are taken on intensities divided by their resolution-shell mean, over the shells whose $\langle I/\sigma\rangle$ reaches 1, with Wilson outliers rejected; the $L$-test pairs are therefore compared on one scale even where two index steps span a steep fall-off, as in phenix.xtriage and ctruncate, and the selection is by shell, never by the individual reflection's $I/\sigma$, which would cut the weak tail and bias $\langle|L|\rangle$ down. The twin fraction is quoted from the statistic that carries the verdict — the $L$-test unless the call rests on the second moment alone — and a second moment is not turned into a fraction under a detected pseudo-translation, which inflates it. A twin law exists only where the **lattice** has a rotation the group lacks, so whether one is possible is read off the cell, not off the name of the Laue class: the lattice's two-folds are found as for the free set (Le Page, within 3° obliquity — GEMMI's `find_twin_laws`), and where none lies outside the group no twin is called. That covers the holohedral classes ($4/mmm$ on a tetragonal lattice, $6/mmm$, $m\bar{3}m$, $\bar{3}m$ on a rhombohedral one) and equally $mmm$ or $2/m$ on a lattice with no special metric, while a pseudo-merohedral $2/m$ or $mmm$ lattice — a monoclinic $\beta$ near 90°, a C-centred cell with $b \approx a\sqrt{3}$ — admits a law and stays eligible. The $L$-test is still read where no law is possible, for a different question: merging $I(h)$ with $I(Th)$ under a false operator $T$ gives $(I(h)+I(Th))/2$ whatever the twin fraction, which has the perfect-twin distribution, so $\langle|L|\rangle$ below 0.42 there is reported as **an adopted operator averaging unequal intensities** — the space group is too high, or a twin law was absorbed into the point group — a warning, never a change to the space group. A genuine operator leaves the untwinned 0.5. Reflections that overlap along a very long axis narrow the distribution the same way, which is one reason it stays a warning. The same numbers measured on the P1 merge the space-group search was given, before any point group was adopted, are reported beside them (`_BEFORE_SEARCH`), with that merge's own pseudo-translation declared and the reflections its lattice centring extinguishes left out. Between that bar and the twin band's 0.44 the verdict is `NO_TWIN_LAW`: the intensities read twin-like, but the lattice admits no law for the group, so they are not a twin in it — the group may be too high (a subgroup's twin law absorbed, or non-crystallographic symmetry taken for real), or something other than twinning narrowed the distribution. It is reported as such rather than as a twin that "needs a twin law" no refinement program could be given.
|
||||
A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$ — 0.500 and 0.333 untwinned, 0.375 and 0.200 for a perfect twin) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ (2.0 for untwinned acentric data, 1.5 for a perfect twin) are written to the merged mmCIF as a twinning diagnostic. Both are taken on intensities divided by their resolution-shell mean, over the shells whose intensities are measured (at least 85% of a shell's reflections at $I/\sigma > 3$, phenix.xtriage's cut for its twin tests; noise pulls $\langle|L|\rangle$ below 0.5 like a twin does), with Wilson outliers rejected; the $L$-test pairs are therefore compared on one scale even where two index steps span a steep fall-off, as in phenix.xtriage and ctruncate, and the selection is by shell, never by the individual reflection's $I/\sigma$, which would cut the weak tail and bias $\langle|L|\rangle$ down. The twin fraction is quoted from the statistic that carries the verdict — the $L$-test unless the call rests on the second moment alone — and a second moment is not turned into a fraction under a detected pseudo-translation, which inflates it. A twin law exists only where the **lattice** has a rotation the group lacks, so whether one is possible is read off the cell, not off the name of the Laue class: the lattice's two-folds are found as for the free set (Le Page, within 3° obliquity — GEMMI's `find_twin_laws`), and where none lies outside the group no twin is called. That covers the holohedral classes ($4/mmm$ on a tetragonal lattice, $6/mmm$, $m\bar{3}m$, $\bar{3}m$ on a rhombohedral one) and equally $mmm$ or $2/m$ on a lattice with no special metric, while a pseudo-merohedral $2/m$ or $mmm$ lattice — a monoclinic $\beta$ near 90°, a C-centred cell with $b \approx a\sqrt{3}$ — admits a law and stays eligible. The $L$-test is still read where no law is possible, for a different question: merging $I(h)$ with $I(Th)$ under a false operator $T$ gives $(I(h)+I(Th))/2$ whatever the twin fraction, which has the perfect-twin distribution, so $\langle|L|\rangle$ below 0.42 there is reported as **an adopted operator averaging unequal intensities** — the space group is too high, or a twin law was absorbed into the point group — a warning, never a change to the space group. A genuine operator leaves the untwinned 0.5. Reflections that overlap along a very long axis narrow the distribution the same way, which is one reason it stays a warning. The same numbers measured on the P1 merge the space-group search was given, before any point group was adopted, are reported beside them (`_BEFORE_SEARCH`), with that merge's own pseudo-translation declared and the reflections its lattice centring extinguishes left out. Between that bar and the twin band's 0.44 the verdict is `NO_TWIN_LAW`: the intensities read twin-like, but the lattice admits no law for the group, so they are not a twin in it — the group may be too high (a subgroup's twin law absorbed, or non-crystallographic symmetry taken for real), or something other than twinning narrowed the distribution. It is reported as such rather than as a twin that "needs a twin law" no refinement program could be given.
|
||||
|
||||
**The twin hypothesis.** Where the adopted group may have absorbed a twin law — a twin gate of the search refused it and was overruled (§13.1), or a group above $P1$ was adopted on a merge that read twin-like before the search — the same observations are also merged in a subgroup and written beside the answer as `<prefix>_twin_<SG>.mtz` / `.cif`, in the subgroup's standard setting, to the answer's resolution and with the answer's free set (each reflection takes the flag of the reflection of the answer's merge its orbit falls in, so twin mates share it). The subgroup is the one a decision named — the group a refused promotion was adopted over, or the indexed lattice's group a metric re-ask promoted from — and otherwise, among the index-2 subgroups of the adopted group that the search confirmed as point-group finalists, the one whose lost two-fold agrees worst with the intensities (lowest operator CC): the most likely twin law. One file is written; the other subgroups that qualified are named in the report. The law is in the mmCIF's `_pdbx_reflns_twin` with its fraction, from Yeates' $H$ on the pairs it relates in the subgroup merge — $H = |I(h)-I(Th)|/(I(h)+I(Th))$ is uniform on $[0, 1-2\alpha]$, so $\alpha = 1/2 - \mathrm{median}\,H$, a lower bound since measurement error only adds to $H$. Nothing the answer is made of depends on it: it is a second merge from the same scaled observations, and the answer's files are byte for byte what they are without it. Refinement programs want twinned data merged in the true lower group, never detwinned, which is why the file is the merge and not a correction of the answer; a near-perfect twin of a subgroup and the higher group predict the same intensities, so which of the two is the crystal is a question for a model (§14.9).
|
||||
|
||||
|
||||
@@ -1048,6 +1048,20 @@ than 50 reflections for the second moment (a few dozen unique reflections, as a
|
||||
gives): ⟨|L|⟩ over N pairs is known to 0.29/√N, so below 50 neither reading separates the
|
||||
untwinned 0.5 from a perfect twin's 0.375.
|
||||
|
||||
Both statistics are read on the resolution shells whose intensities are measured - at least 85% of a
|
||||
shell's reflections at I/σ above 3, the cut phenix.xtriage applies to its twin tests - because noise
|
||||
moves ⟨|L|⟩ the way a twin does: a pair of noisy intensities kept at I > 0 reads below 0.5, and shells
|
||||
admitted on a mean I/σ of 1 pulled untwinned crystals into the twin band.
|
||||
|
||||
**`TWIN_LAW_n=`** lists every twin law the lattice admits for the group (an hkl operator in the data's
|
||||
basis, the one a twin refinement takes), each with the number of reflection pairs it relates, their
|
||||
median disagreement H = |I₁−I₂|/(I₁+I₂), the twin fraction 0.5 − H that implies if only twinning makes
|
||||
the pairs agree (the pairs are the stronger reflections, and untwinned crystals read up to about 0.25
|
||||
on it), their intensity-weighted R and the correlation of their normalised intensities. Where the
|
||||
space-group search refused a higher group on its twin-immune zone, **`TWIN_LAW=`** headlines the law
|
||||
of that list whose pairs agree best, at its own fraction, and the mmCIF carries it as
|
||||
`_pdbx_reflns_twin`.
|
||||
|
||||
### The twin hypothesis (`TWIN_SUBGROUP_*`)
|
||||
|
||||
Where the adopted space group may have absorbed a twin law - a twin gate of the space-group search
|
||||
|
||||
@@ -949,7 +949,13 @@ RowModulation FindRowModulation(const gemmi::UnitCell& cell, const std::array<in
|
||||
for (size_t i = 0; i < m.n; ++i) {
|
||||
if (!m.pass_cc[i] || visited[i] == epoch)
|
||||
continue;
|
||||
const auto m2 = op.apply_to_hkl(gemmi::Op::Miller{{m.H[i], m.K[i], m.L[i]}});
|
||||
// A pseudo-merohedral law of a centred lattice has half-integer entries in the lattice's own
|
||||
// basis (OperatorScorer), so it maps an index outside the centring onto non-integers; such an
|
||||
// image is no reflection, and apply_to_hkl would silently truncate it onto one.
|
||||
const auto scaled = op.apply_to_hkl_without_division(gemmi::Op::Miller{{m.H[i], m.K[i], m.L[i]}});
|
||||
if (scaled[0] % gemmi::Op::DEN != 0 || scaled[1] % gemmi::Op::DEN != 0 || scaled[2] % gemmi::Op::DEN != 0)
|
||||
continue;
|
||||
const auto m2 = gemmi::Op::divide_hkl_by_DEN(scaled);
|
||||
const HKLKey k2 = Canonicalize(m2[0], m2[1], m2[2], opt.merge_friedel);
|
||||
if (k2 == m.key[i])
|
||||
continue; // reflection lies on this rotation axis
|
||||
@@ -3220,12 +3226,16 @@ std::optional<SpaceGroupOperatorScore> OperatorScorer::operator()(const gemmi::M
|
||||
// gemmi applies an Op to Miller indices through the TRANSPOSE of its rot (Op::apply_to_hkl), and
|
||||
// writes an hkl triplet from the transpose as well, so the operator whose Miller-index matrix is M
|
||||
// is the Op holding M^T - and its triplet then reads as the rows of M, as everywhere else here.
|
||||
// The entries are kept in gemmi's 1/DEN units, not rounded to integers first: a pseudo-merohedral
|
||||
// law of a centred lattice (a C-centred cell with a = sqrt(3) b, say) has half-integer entries in
|
||||
// that basis, and rounding them scored a different operator - its pairs read unrelated (H 0.61, CC
|
||||
// 0.02) on a crystal twinned by it at 0.3.
|
||||
gemmi::Op op{};
|
||||
op.notation = 'h';
|
||||
op.tran = {0, 0, 0};
|
||||
for (int i = 0; i < 3; ++i)
|
||||
for (int j = 0; j < 3; ++j)
|
||||
op.rot[i][j] = static_cast<int>(std::lround(hkl_matrix[j][i])) * gemmi::Op::DEN;
|
||||
op.rot[i][j] = static_cast<int>(std::lround(hkl_matrix[j][i] * gemmi::Op::DEN));
|
||||
|
||||
const NormalisedMerge &nm = *impl->nm;
|
||||
std::vector<uint32_t> visited(nm.n, 0);
|
||||
|
||||
@@ -229,16 +229,24 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector<MergedReflection>& merg
|
||||
// on the same scale. 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 these twinning indicators rather than data-quality artefacts -
|
||||
// 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.
|
||||
// a whole dataset's value). To keep these twinning indicators rather than data-quality artefacts
|
||||
// we 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, and read only the shells whose
|
||||
// intensities are MEASURED: at least 85% of a shell's reflections at I/sigma above 3, the cut
|
||||
// phenix.xtriage applies to its twin tests. Both statistics need that. Noise moves <|L|> the way
|
||||
// a twin does: a pair of noisy intensities kept at I > 0 reads |I1-I2|/(I1+I2) below the untwinned
|
||||
// 0.5 (simulated on Wilson intensities, <|L|> = 0.474 at sigma/<I> = 0.1, 0.449 at 0.3, 0.437 at
|
||||
// 0.5, 0.426 at 1), so shells admitted at a mean I/sigma of 1 - the rule before - pulled an
|
||||
// untwinned crystal's <|L|> into the twin band: five untwinned crystals read 0.413-0.440 over
|
||||
// their whole range against 0.448-0.492 over the shells xtriage reads, and were called twins,
|
||||
// while the twinned ones read within 0.01 either way.
|
||||
//
|
||||
// The shells are equal in width in 1/d^2. Equal-COUNT shells were measured as the alternative
|
||||
// (the first equal-width shell spans a wide range of d, and its fall-off inflates the moment):
|
||||
// they move the moment the wrong way on the same data, widening its gap to phenix.xtriage by
|
||||
// 0.01-0.04 rather than closing it, so the shells were kept as they are.
|
||||
constexpr double min_shell_isig = 1.0; // shells below this are noise, not signal
|
||||
constexpr double measured_i_over_sigma = 3.0; // a reflection counts as measured above this
|
||||
constexpr double min_measured_fraction = 0.85; // and a shell is read when this many of its are
|
||||
constexpr double wilson_outlier_e2 = 8.0; // reject improbably strong reflections (P ~ e^-8)
|
||||
const int n_shells = std::max(1, resolution_shells);
|
||||
double min_s = std::numeric_limits<double>::infinity();
|
||||
@@ -256,26 +264,36 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector<MergedReflection>& 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<int>(t * n_shells)));
|
||||
};
|
||||
// Group intensities by shell, and accumulate <I/sigma> to gauge each shell's signal.
|
||||
// Group intensities by shell, and count how many of each shell's are measured.
|
||||
std::vector<std::vector<double>> shell_I(n_shells);
|
||||
std::vector<int> shell_with_sigma(n_shells, 0), shell_measured(n_shells, 0);
|
||||
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 (!used(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_with_sigma[b] += 1;
|
||||
shell_isig_sum[b] += r.I / r.sigma;
|
||||
shell_isig_n[b] += 1;
|
||||
if (r.I / r.sigma > measured_i_over_sigma)
|
||||
shell_measured[b] += 1;
|
||||
}
|
||||
}
|
||||
// The mean intensity of each shell that carries signal; 0 marks a shell that does not.
|
||||
std::vector<char> shell_read(n_shells, 0);
|
||||
for (int b = 0; b < n_shells; ++b)
|
||||
shell_read[b] = shell_with_sigma[b] > 0 && shell_measured[b] >= min_measured_fraction * shell_with_sigma[b];
|
||||
// Weak data can have no measured shell at all. Then the statistics are read as they were before,
|
||||
// on the shells whose mean I/sigma reaches 1, rather than not at all: the reading is the one such
|
||||
// data always gave, and a verdict is still better than a report without the statistics.
|
||||
if (std::none_of(shell_read.begin(), shell_read.end(), [](char c) { return c != 0; }))
|
||||
for (int b = 0; b < n_shells; ++b)
|
||||
shell_read[b] = shell_with_sigma[b] > 0 && shell_isig_sum[b] >= 1.0 * shell_with_sigma[b];
|
||||
// The mean intensity of each shell that is read; 0 marks a shell that is not.
|
||||
std::vector<double> shell_mean(n_shells, 0.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)
|
||||
if (intensities.empty() || !shell_read[b])
|
||||
continue;
|
||||
double sum = 0.0;
|
||||
for (double I : intensities)
|
||||
|
||||
@@ -2608,8 +2608,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
|
||||
if (!result.twinning.twin_law_operator.empty()) {
|
||||
Add(s, KeyText("TWIN_LAW", fmt::format("{} fraction {:.2f}", result.twinning.twin_law_operator,
|
||||
result.twinning.twin_law_fraction)));
|
||||
Add(s, Prose(fmt::format(" Twin law identified: operator {} at a fraction of about {:.2f} (a lower bound from\n"
|
||||
" the operator's own disagreement). The mmCIF carries it as _pdbx_reflns_twin; refine\n"
|
||||
Add(s, Prose(fmt::format(" Twin law identified: operator {} at a fraction of about {:.2f} (0.5 - H of the\n"
|
||||
" pairs it relates, as TWIN_LAW_n reads it). The mmCIF carries it as _pdbx_reflns_twin; refine\n"
|
||||
" against these data in this space group with that law.",
|
||||
result.twinning.twin_law_operator, result.twinning.twin_law_fraction)));
|
||||
}
|
||||
|
||||
@@ -4505,7 +4505,11 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
|
||||
// would not do: the call replaces the whole struct.)
|
||||
result.twinning.laue_class_was_chosen_by_promotion = promoted_point_group;
|
||||
// A refusal that stood on the twin-immune zone names the law: the first operator the
|
||||
// refused group adds, at the fraction its own disagreement implies (a lower bound).
|
||||
// refused group adds, at the fraction its own disagreement implies (a lower bound). Both are
|
||||
// replaced below by the law whose own pairs agree best, once every law of the lattice has been
|
||||
// read: the fraction here is the mean H over ALL the operators the refused group added, or the
|
||||
// L-test's, and neither is the headlined law's own - a 3-fold crystal twinned at 0.4 by one of
|
||||
// its three laws was headlined at 0.20 beside its law's own H reading 0.40.
|
||||
if (twin_sg && result.space_group_search
|
||||
&& result.space_group_search->refused_point_group_representative
|
||||
&& result.twin_zone_refusal_stood) {
|
||||
@@ -4628,11 +4632,14 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
|
||||
};
|
||||
std::optional<OperatorScorer> law_correlation;
|
||||
std::string laws;
|
||||
std::string best_law; // the law whose pairs agree best, and its H
|
||||
double best_law_h = 1.0;
|
||||
for (const auto &op : twin_ops) {
|
||||
std::string t = op.triplet();
|
||||
std::replace(t.begin(), t.end(), 'x', 'h');
|
||||
std::replace(t.begin(), t.end(), 'y', 'k');
|
||||
std::replace(t.begin(), t.end(), 'z', 'l');
|
||||
// find_twin_laws returns real-space operations; the Miller-index operator is the
|
||||
// transpose (an hkl triplet, as TWIN_ZONE and TWIN_LAW print theirs), and the two differ
|
||||
// for every law that is not symmetric - a pseudo-merohedral one of a centred lattice,
|
||||
// say, whose real-space triplet relates unrelated reflections when read as hkl.
|
||||
const std::string t = op.as_hkl().triplet();
|
||||
laws += (laws.empty() ? "" : " ; ") + t;
|
||||
// What each law reads on the merge: the disagreement H of the reflections it
|
||||
// relates (Yeates (1997) Methods Enzymol. 276, 344-358: a twin of fraction a reads
|
||||
@@ -4649,16 +4656,27 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
|
||||
law_opts.nthreads = static_cast<size_t>(std::max(1, config_.nthreads));
|
||||
if (!law_correlation)
|
||||
law_correlation.emplace(expanded_merge(), law_opts);
|
||||
if (const auto s = (*law_correlation)(m); s.has_value())
|
||||
if (const auto s = (*law_correlation)(m); s.has_value()) {
|
||||
result.twin_law_evidence.push_back(fmt::format(
|
||||
"{}: {} pairs, H {:.3f} (twin fraction {:.2f} if only twinning makes the pairs "
|
||||
"agree; untwinned crystals read up to about 0.25), R {:.3f}, CC {:.3f}",
|
||||
t, s->n_pairs, s->h_stat, std::max(0.0, 0.5 - s->h_stat), s->r_stat, s->cc));
|
||||
else
|
||||
if (s->h_stat < best_law_h) {
|
||||
best_law_h = s->h_stat;
|
||||
best_law = t;
|
||||
}
|
||||
} else
|
||||
result.twin_law_evidence.push_back(fmt::format("{}: too few pairs to read", t));
|
||||
}
|
||||
for (const auto &e : result.twin_law_evidence)
|
||||
logger.Info("Twin law {}", e);
|
||||
// The twin law a refusal identified (above) is the one of these laws whose pairs agree
|
||||
// best - a twin law makes its own pairs agree, the lattice's other laws relate unrelated
|
||||
// intensities - at the fraction its own H implies, the number TWIN_LAW_n prints for it.
|
||||
if (!result.twinning.twin_law_operator.empty() && !best_law.empty()) {
|
||||
result.twinning.twin_law_operator = best_law;
|
||||
result.twinning.twin_law_fraction = std::max(0.0, 0.5 - best_law_h);
|
||||
}
|
||||
// A model can resolve it too, but whether it does is only known once the model has
|
||||
// been fitted (validation step below), so the warning stands here exactly as in a run
|
||||
// without a model, and is withdrawn there if the model's indexing probe decided it.
|
||||
|
||||
@@ -1025,3 +1025,38 @@ TEST_CASE("AnalyzeCentreOfSymmetry separates acentric from centric Wilson data")
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// A pseudo-merohedral twin law of a centred lattice has half-integer entries in the lattice's own
|
||||
// basis (a C-centred cell with a = sqrt(3) b is metrically hexagonal, and the hexagonal two-folds read
|
||||
// as h/2 +- 3k/2 ...): it maps the centring-present indices onto integers and the absent ones onto
|
||||
// halves. Scored as given it relates its own pairs; rounded to integers it is a different operator.
|
||||
TEST_CASE("Operator correlation of a half-integer pseudo-merohedral law") {
|
||||
const gemmi::Op law = gemmi::parse_triplet("h/2-k/2,-3/2*h-k/2,-l", 'h');
|
||||
gemmi::Mat33 m; // the Miller-index matrix of the law: h -> M h, rows as the triplet reads
|
||||
for (int i = 0; i < 3; ++i)
|
||||
for (int j = 0; j < 3; ++j)
|
||||
m[i][j] = static_cast<double>(law.rot[j][i]) / gemmi::Op::DEN;
|
||||
// Friedel-merged C-centred reflections whose intensity is a function of the orbit under the law.
|
||||
std::vector<MergedReflection> merged;
|
||||
for (int h = -8; h <= 8; ++h)
|
||||
for (int k = -8; k <= 8; ++k)
|
||||
for (int l = -8; l <= 8; ++l) {
|
||||
if ((h + k) % 2 != 0 || (h == 0 && k == 0 && l == 0))
|
||||
continue;
|
||||
const gemmi::Op::Miller hkl{{h, k, l}};
|
||||
const gemmi::Op::Miller mate = law.apply_to_hkl(hkl);
|
||||
const gemmi::Op::Miller neg{{-h, -k, -l}}, mate_neg{{-mate[0], -mate[1], -mate[2]}};
|
||||
if (hkl < neg)
|
||||
continue; // one Friedel mate
|
||||
const gemmi::Op::Miller rep = std::max({hkl, neg, mate, mate_neg});
|
||||
merged.push_back(MergedReflection{.h = h, .k = k, .l = l,
|
||||
.I = static_cast<float>(SyntheticIntensityFromAsu(rep)),
|
||||
.sigma = 1.0, .d = static_cast<float>(CalcSyntheticD(h, k, l))});
|
||||
}
|
||||
const auto score = OperatorCorrelation(merged, m);
|
||||
REQUIRE(score.has_value());
|
||||
CHECK(score->n_pairs > 300);
|
||||
CHECK(score->cc > 0.99); // on shell-normalised E^2, and the synthetic d does not follow the law exactly
|
||||
CHECK(score->h_stat == Catch::Approx(0.0).margin(1e-6));
|
||||
CHECK(score->op_triplet_hkl == "h/2-k/2,-3/2*h-k/2,-l");
|
||||
}
|
||||
Reference in new issue
Block a user