diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 705aa8663..398cca775 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -17,6 +17,7 @@ * `rugnux` fits the direction of the goniometer axis and not its length, so the cell chosen by the first pass is the one its own refinement scored. * `rugnux` reports the detector geometry it measured - the direct beam, the tilt and the beam centre - and what a single sweep can and cannot determine. * `rugnux` reports the resolution the CC1/2 fit reached and the strong-direction diffraction limit. +* The rugnux anisotropy warning says which reciprocal axis each of the two diffraction limits it quotes belongs to. * `rugnux` reports twinning measured before and after the space group was decided, and no longer reports it when the L-test contradicts it. * `rugnux --mode calibration` writes `.json` beside the `.poni`, holding the geometry as a `jfjoch_broker` `dataset_settings` body. * `rugnux --mode calibration` says when a fit is not a measurement: it writes no `.poni`, exits non-zero, and records `converged` in the `.json`. diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index ccb34fcc9..66dd9e8b3 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -37,6 +37,30 @@ namespace { return fmt::format("{:.3f} {:.3f} {:.3f} {:.3f} {:.3f} {:.3f}", c.a, c.b, c.c, c.alpha, c.beta, c.gamma); } + + // Which reciprocal axis a principal anisotropy direction lies along. The tensor is fitted on + // s = frac.mat * (h,k,l), so in that Cartesian frame a*, b*, c* are the rows of frac.mat and the + // eigenvector is named by whichever it makes the smallest angle with. The label is exact in every + // Laue class the tensor has a free direction in except triclinic; the cosine says how well it fits. + std::string ReciprocalAxisLabel(const UnitCell &cell, const double v[3]) { + static const char *NAME[3] = {"a*", "b*", "c*"}; + const gemmi::UnitCell gc = cell; + int best = 0; + double best_cos = -1.0; + for (int i = 0; i < 3; ++i) { + const gemmi::Vec3 axis = gc.frac.mat.left_multiply( + gemmi::Vec3(i == 0 ? 1.0 : 0.0, i == 1 ? 1.0 : 0.0, i == 2 ? 1.0 : 0.0)); + const double len = axis.length(); + if (!(len > 0.0)) + continue; + const double c = std::fabs((axis.x * v[0] + axis.y * v[1] + axis.z * v[2]) / len); + if (c > best_cos) { + best_cos = c; + best = i; + } + } + return fmt::format("{} (cos {:.2f})", NAME[best], best_cos); + } } std::string RenderResultReport(const std::string &output_prefix, @@ -521,13 +545,24 @@ std::string RenderResultReport(const std::string &output_prefix, "fall-off the resolution cut is read off does not describe these data", sh[rose_after].d_max, sh[rose_after].d_min)); } - if (an.verdict == AnisotropyVerdict::Detected && an.d_min_spread > 0.5) + if (an.verdict == AnisotropyVerdict::Detected && an.d_min_spread > 0.5) { + // Say WHICH direction each limit belongs to. Without it the warning states that the crystal + // is anisotropic and leaves the reader no way to act on it; the eigenvectors are measured + // here and reach the mmCIF, so the name costs nothing. + const double *worst = std::max_element(an.d_min_axis, an.d_min_axis + 3); + const double *best = std::min_element(an.d_min_axis, an.d_min_axis + 3); + const auto along = [&](const double *it) { + const int n = static_cast(it - an.d_min_axis); + return result.consensus_cell + ? ReciprocalAxisLabel(*result.consensus_cell, an.eigenvector[n]) + : fmt::format("principal direction {}", n + 1); + }; warnings.emplace_back(fmt::format( - "Diffraction is anisotropic (deltaB {:.1f} A^2; the diffraction limit runs from " - "{:.2f} to {:.2f} A depending on direction) - refinement and map interpretation " - "should allow for it; no intensity has been corrected for it here", - an.delta_b, *std::max_element(an.d_min_axis, an.d_min_axis + 3), - *std::min_element(an.d_min_axis, an.d_min_axis + 3))); + "Diffraction is anisotropic (deltaB {:.1f} A^2; the diffraction limit is {:.2f} A " + "along {} and {:.2f} A along {}) - refinement and map interpretation should allow " + "for it; no intensity has been corrected for it here", + an.delta_b, *worst, along(worst), *best, along(best))); + } } // ------------------------------------------------------- 10. MODEL VALIDATION