diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 826d5ba7..15348a57 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -3,6 +3,7 @@ #include "WriteReflections.h" #include "scale_merge/Merge.h" +#include "scale_merge/TwinningAnalysis.h" #include #include @@ -62,6 +63,7 @@ void WriteMmcifReflections(const std::vector &reflections, const DiffractionExperiment &experiment, const MergeStatistics &statistics, const std::string &isa, + const TwinningAnalysisResult &twinning, const std::string &filename) { std::ofstream out(filename); @@ -136,6 +138,16 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_reflns.pdbx_Rrim_I_all " << Fmt(ov.r_meas, 4) << "\n"; out << "_reflns.pdbx_CC_half " << Fmt(ov.cc_half, 4) << "\n"; out << "_reflns.pdbx_diffrn_ISa " << CifStr(isa) << " # asymptotic I/sigma (Diederichs)\n"; + // Twinning indicators (no standard CIF item, written as free-text pdbx values like ISa above). + if (twinning.l_test_pairs > 0) { + out << "_reflns.pdbx_L_test_mean_abs_L " << Fmt(twinning.mean_abs_l, 3) + << " # Padilla-Yeates <|L|> (untwinned 0.500, perfect twin 0.375)\n"; + out << "_reflns.pdbx_L_test_mean_L_squared " << Fmt(twinning.mean_l_squared, 3) + << " # (untwinned 0.333, perfect twin 0.200)\n"; + } + if (twinning.moment_reflections > 0) + out << "_reflns.pdbx_second_moment_I " << Fmt(twinning.second_moment, 3) + << " # /^2 (untwinned 2.00, perfect twin 1.50)\n"; out << "#\n"; out << "loop_\n"; @@ -249,13 +261,14 @@ void WriteReflections(const std::vector &reflections, const DiffractionExperiment &experiment, const MergeStatistics &statistics, const std::string &isa, + const TwinningAnalysisResult &twinning, const std::string &filename) { switch (experiment.GetScalingSettings().GetFileFormat()) { case IntensityFormat::Text: WriteHKLReflections(reflections, filename + ".hkl"); break; case IntensityFormat::mmCIF: - WriteMmcifReflections(reflections, unitCell, experiment, statistics, isa, filename + ".cif"); + WriteMmcifReflections(reflections, unitCell, experiment, statistics, isa, twinning, filename + ".cif"); break; case IntensityFormat::MTZ: WriteMtzReflections(reflections, unitCell, experiment, filename + ".mtz"); diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 8215dd93..4c03d092 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -11,12 +11,14 @@ #include "../common/DiffractionExperiment.h" struct MergeStatistics; +struct TwinningAnalysisResult; void WriteMmcifReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, const MergeStatistics &statistics, const std::string &isa, + const TwinningAnalysisResult &twinning, const std::string &filename); void WriteMtzReflections(const std::vector &reflections, @@ -32,4 +34,5 @@ void WriteReflections(const std::vector &reflections, const DiffractionExperiment &experiment, const MergeStatistics &statistics, const std::string &isa, + const TwinningAnalysisResult &twinning, const std::string &filename); \ No newline at end of file diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp index 97b70c30..8dd392d7 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -49,7 +49,14 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector& merg if (UsableIntensity(r) && r.I > 0.0 && acentric(r)) intensity.emplace(PackHKL(r.h, r.k, r.l), r.I); - const std::array, 3> steps{{{2, 0, 0}, {0, 2, 0}, {0, 0, 2}}}; + // 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, 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; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 80b5dd9c..fae65aa4 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -609,7 +609,7 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { phase("Writing reflections"); WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, result.error_model_isa > 0 ? fmt::format("{:.2f}", result.error_model_isa) : "?", - config_.output_prefix); + result.twinning, config_.output_prefix); // Per-image scaling table (G, B-factor, mosaicity, wedge, CC) for inspection / XDS // comparison. The offline self-scaling result is otherwise not exposed (process.h5's // per-image arrays are only filled on the online per-image path). Sourced from the diff --git a/tools/rugnux_cli.cpp b/tools/rugnux_cli.cpp index bb3563ad..49c12685 100644 --- a/tools/rugnux_cli.cpp +++ b/tools/rugnux_cli.cpp @@ -1014,12 +1014,13 @@ int main(int argc, char **argv) { const auto twin_sg_number = experiment.GetSpaceGroupNumber(); const gemmi::SpaceGroup *twin_sg = twin_sg_number ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr; - std::cout << std::endl << TwinningAnalysisToText(AnalyzeTwinning(merged_reflections, twin_sg)) << std::endl; + const auto twinning = AnalyzeTwinning(merged_reflections, twin_sg); + std::cout << std::endl << TwinningAnalysisToText(twinning) << std::endl; if (!output_prefix.empty()) WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics, error_model_isa > 0 ? fmt::format("{:.2f}", error_model_isa) : "?", - output_prefix); + twinning, output_prefix); return 0; }