diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 71a9a68b..e9ecfafd 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -291,6 +291,60 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { return RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); } +namespace { + // Angle between the crystal's symmetry axes and the spindle. A rotation sweep never records the + // reflections whose reciprocal vector lies close to the spindle - the blind cusp - and normally + // symmetry fills it in from an equivalent elsewhere in reciprocal space. It cannot when a symmetry + // axis IS the spindle: the cusp then maps onto itself, every reflection in it is equivalent only to + // other reflections in it, and the hole stays empty however long the sweep runs. Returns the + // smallest angle in degrees between any proper rotation axis of the space group and the goniometer + // axis, with the order of the axis that achieves it. + std::optional> ClosestSymmetryAxisToSpindle(const gemmi::SpaceGroup &sg, + const CrystalLattice &lattice, + const Coord &spindle) { + const double spindle_len = spindle.Length(); + if (spindle_len < 1e-9) + return std::nullopt; + std::optional> best; + for (const auto &op : sg.operations().derive_symmorphic().sym_ops) { + if (op.rot == gemmi::Op::identity().rot) + continue; + // Order of the rotation, then project onto its invariant direction by summing its powers: + // (1/n) sum_k W^k annihilates everything except the axis. + gemmi::Op acc = gemmi::Op::identity(), w{op.rot, {0, 0, 0}, op.notation}; + int order = 0; + gemmi::Op cur = gemmi::Op::identity(); + std::array, 3> sum{}; + for (int k = 0; k < 6; ++k) { + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + sum[i][j] += static_cast(cur.rot[i][j]) / gemmi::Op::DEN; + cur = cur.combine(w).wrap(); + ++order; + if (cur.rot == gemmi::Op::identity().rot) + break; + } + if (order < 2) + continue; + // Any non-degenerate column of the projector is the axis in fractional direct coordinates. + for (int c = 0; c < 3; ++c) { + const Coord axis = lattice.Vec0() * static_cast(sum[0][c]) + + lattice.Vec1() * static_cast(sum[1][c]) + + lattice.Vec2() * static_cast(sum[2][c]); + const double len = axis.Length(); + if (len < 1e-6 * lattice.Vec0().Length()) + continue; + const double cosang = std::abs((axis * spindle) / (len * spindle_len)); + const double ang = std::acos(std::min(1.0, cosang)) * 180.0 / M_PI; + if (!best.has_value() || ang < best->first) + best = std::make_pair(ang, order); + break; + } + } + return best; + } +} + ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, bool geometry_prepass) { Logger logger("Rugnux"); ProcessResult result; @@ -1370,6 +1424,39 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.twinning.laue_class_was_chosen_by_promotion = promoted_point_group; stats_text << TwinningAnalysisToText(result.twinning) << "\n"; + // Symmetry axis vs spindle. Reported always on rotation data (it is a property of how the + // crystal was mounted, which the user can change), warned about when the two nearly coincide. + if (experiment_.IsRotationIndexing() && twin_sg && end_msg.rotation_lattice.has_value()) { + if (const auto gonio = experiment_.GetGoniometer()) { + const auto closest = ClosestSymmetryAxisToSpindle(*twin_sg, *end_msg.rotation_lattice, + gonio->GetAxis()); + if (closest.has_value()) { + // Practical bound, not a derived one: the blind cone's half-angle is the maximum + // Bragg angle (~15 deg for 2 A data at 1 A), and measured on this battery a 13.6 deg + // case shows the loss while a 16.2 deg one is 99.7% complete. + constexpr double WARN_DEG = 15.0; + const auto [angle, order] = *closest; + if (angle < WARN_DEG) { + const std::string msg = fmt::format( + "The crystal's {}-fold axis is only {:.1f} deg from the spindle. A rotation sweep " + "never records the reflections whose reciprocal vector lies within the Bragg angle " + "of the spindle, and symmetry normally supplies them from an equivalent elsewhere; " + "it cannot here, because that axis maps the blind region onto itself. Those " + "reflections stay missing however long the sweep runs. The loss is confined to " + "that cone rather than spread over the data, so overall completeness may still " + "look reasonable - check it near the spindle direction. A second sweep on a " + "different axis, or re-mounting, recovers them.", + order, angle); + logger.Warning("{}", msg); + stats_text << " !! " << msg << "\n\n"; + } else { + stats_text << "Closest symmetry axis to the spindle: " << order << "-fold at " + << std::fixed << std::setprecision(1) << angle << " deg\n\n"; + } + } + } + } + // Indexing-ambiguity (alternative-indexing) advisory. When the lattice metric symmetry exceeds // the Laue symmetry the crystal can be validly indexed in several hands related by twin-law // operators. For an obvious merohedral case (P3/P4/P6...) users expect this; but a PSEUDO-merohedral