diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 19cb14139..b5a4e5b3f 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -10,6 +10,7 @@ * `rugnux` writes a batch header in the unmerged MTZ for every image the observations span, not only for the images that produced one, so a scaling program reads one run per sweep. * `rugnux` leaves an event out of the unmerged MTZ when less of its rocking curve was captured than `--min-captured-fraction`, as the merge does, since the exported rows are declared full. * `rugnux` writes `FreeR_flag` with 0 for the test set and 1 for the working set, the CCP4 convention REFMAC5 defaults to; it was the other way round. +* `rugnux -S` refuses or re-seats a fixed space group whose symmetry axes the indexed cell does not carry, not only one whose centring differs; `--mode scale` refuses it too. * `rugnux` handles symmetry better: the lattice, the point group, the setting and the systematic absences. * `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. diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 5370496a1..5a6fe21fd 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -64,6 +64,34 @@ #include #include +double MetricViolation(const UnitCell &uc, const gemmi::SpaceGroup &sg) { + const double a = uc.a, b = uc.b, c = uc.c; + const double d2r = PI / 180.0; + const double ab = a * b * std::cos(uc.gamma * d2r); + const double ac = a * c * std::cos(uc.beta * d2r); + const double bc = b * c * std::cos(uc.alpha * d2r); + const double g[3][3] = {{a * a, ab, ac}, {ab, b * b, bc}, {ac, bc, c * c}}; + const double scale = std::max({a * a, b * b, c * c}); + if (!(scale > 0.0)) + return 0.0; + double worst = 0.0; + for (const gemmi::Op &op : sg.operations()) { + double r[3][3]; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + r[i][j] = static_cast(op.rot[i][j]) / gemmi::Op::DEN; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) { + double t = 0.0; + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + t += r[k][i] * g[k][l] * r[l][j]; + worst = std::max(worst, std::fabs(t - g[i][j]) / scale); + } + } + return worst; +} + namespace { // A CUDA out-of-memory, or a host allocation that failed, says nothing about the image that // happened to be in flight and everything about the machine - the next image will hit it too. @@ -3332,15 +3360,21 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // observations thrown away (2.14 M -> 1.08 M) and ISa 28.3 -> 24.6. Ask the character table for the // group's own class instead and re-seat the reflections into it. // - // Keyed on the CENTRING differing, which is the whole of the harm: a fixed group of lower symmetry - // than the lattice but with the same centring (-S P21 on an orthorhombic-P lattice, say) merges - // correctly where it stands, and re-seating it would only permute axes that are already right. - // A group the RUN determined is left alone - the arm above and the centring check below are its - // handling, and a de-novo run must behave exactly as it did. + // Two ways an indexed lattice can fail to carry the fixed group's Bravais lattice, and one + // re-seating for both. The CENTRING can differ, as above. Or the centring can agree while the + // METRIC does not - the axes permuted, a=37.909 b=78.031 c=77.594 declared P 41 21 2, which + // needs a=b with the 4-fold along c. That one is invisible to a centring test and merges + // anyway, in a frame the reflections are not in: measured, --mode scale -S P41212 on such a + // run reports COMPLETENESS= 195.3. + // A fixed group of lower symmetry than the lattice but with the same centring (-S P21 on an + // orthorhombic-P lattice, say) satisfies both tests where it stands, merges correctly, and is + // left alone. So is a group the RUN determined - it was determined FROM this lattice, and the + // arm above and the check below are its handling; a de-novo run must behave exactly as it did. if (user_fixed_sg_ && end_msg.rotation_lattice.has_value() && end_msg.rotation_lattice_type.has_value()) { const gemmi::SpaceGroup &fixed_sg = *user_fixed_sg_; const char indexed_centering = end_msg.rotation_lattice_type->centering; - if (indexed_centering != fixed_sg.centring_type()) { + const double violation = end_msg.unit_cell ? MetricViolation(*end_msg.unit_cell, fixed_sg) : 0.0; + if (indexed_centering != fixed_sg.centring_type() || violation > MAX_METRIC_VIOLATION) { // A trigonal-P group sits on a hexagonal-P lattice; every other system names its own. const auto want_system = fixed_sg.crystal_system() == gemmi::CrystalSystem::Trigonal && fixed_sg.centring_type() == 'P' @@ -3359,9 +3393,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b fixed_sg.centring_type()); if (cand && reindex_into(*cand)) { const auto &uc = *result.consensus_cell; - logger.Info("Reindexed the {}-centred indexed lattice into the {}-centred setting the " - "fixed space group {} needs: a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} " - "beta={:.2f} gamma={:.2f}", indexed_centering, cand->centering, + logger.Info("Reindexed the indexed {}-centred lattice (metric violation {:.3f} against " + "the group's own operators) into the {}-centred setting the fixed space " + "group {} needs: a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} beta={:.2f} " + "gamma={:.2f}", indexed_centering, violation, cand->centering, fixed_sg.xhm(), uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma); } } @@ -3382,25 +3417,34 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // question that was not asked. Name the cell that WAS indexed so the user can act on it. if (user_fixed_sg_ && end_msg.rotation_lattice_type.has_value() && end_msg.unit_cell.has_value()) { const gemmi::SpaceGroup &fixed_sg = *user_fixed_sg_; - if (end_msg.rotation_lattice_type->centering != fixed_sg.centring_type()) { + const bool centring_wrong = + end_msg.rotation_lattice_type->centering != fixed_sg.centring_type(); + // The same two tests as the re-seating above, asked again now that it has had its chance. + const double violation = MetricViolation(*end_msg.unit_cell, fixed_sg); + if (centring_wrong || violation > MAX_METRIC_VIOLATION) { const auto &uc = *end_msg.unit_cell; + const std::string why = centring_wrong + ? fmt::format("its lattice is {}-centred and this crystal indexes as {}-centred", + fixed_sg.centring_type(), end_msg.rotation_lattice_type->centering) + : fmt::format("this crystal's cell is not one the group's own operators leave " + "invariant (metric violation {:.3f}, tolerance {:.3f}) - the symmetry " + "axes do not run along the axes the group puts them on", + violation, MAX_METRIC_VIOLATION); if (prepass_result_) { - logger.Warning("Two-pass: the second pass indexed a {}-centred lattice, but the fixed " - "space group {} has a {}-centred lattice - the group cannot describe " - "this lattice, so this pass will not be adopted", - end_msg.rotation_lattice_type->centering, fixed_sg.xhm(), - fixed_sg.centring_type()); + logger.Warning("Two-pass: the second pass indexed a lattice the fixed space group {} " + "cannot describe - {} - so this pass will not be adopted", + fixed_sg.xhm(), why); result.lattice_conflicts_with_fixed_sg = true; } else { throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, fmt::format( - "The space group {} was fixed for this run, and its lattice is {}-centred - but " - "this crystal indexes as a {}-centred {} lattice (a={:.3f} b={:.3f} c={:.3f} " - "alpha={:.2f} beta={:.2f} gamma={:.2f}), and no setting of that lattice carries the " - "group's. Merging in {} would apply its absence rule to reflections that are not in " - "its frame, so nothing it reported would describe this crystal. Re-run without a fixed " - "space group to have it determined from the data, or fix one whose lattice this " - "crystal has.", - fixed_sg.xhm(), fixed_sg.centring_type(), end_msg.rotation_lattice_type->centering, + "The space group {} was fixed for this run, but {}: the crystal indexes as a " + "{}-centred {} lattice (a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} beta={:.2f} " + "gamma={:.2f}), and no setting of that lattice carries the group's. Merging in {} " + "would apply its symmetry to reflections that are not in its frame, so nothing it " + "reported would describe this crystal - a completeness above 100% is the usual " + "symptom. Re-run without a fixed space group to have it determined from the data, " + "or fix one whose lattice this crystal has.", + fixed_sg.xhm(), why, end_msg.rotation_lattice_type->centering, gemmi::crystal_system_str(end_msg.rotation_lattice_type->crystal_system), uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, fixed_sg.xhm())); } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 5e6865b55..721aec11a 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -319,6 +319,25 @@ struct ProcessResult { std::vector warnings; }; +// How far a cell stands from the metric its space group requires. A group's own rotations must leave +// the metric tensor G invariant - R^T G R = G - so a cell that fails that cannot be described by that +// group whatever its centring is: -S P41212 on a=37.909 b=78.031 c=77.594 has the right centring and +// the wrong axis, the 4-fold running along a where the group puts it along c, and the group's +// operators do not act on the reflections' own indices. Returned relative to the largest diagonal +// element of G, so it is dimensionless and comparable across cells. +double MetricViolation(const UnitCell &uc, const gemmi::SpaceGroup &sg); + +// The metric never fits exactly - the cell is refined against the data, not constrained to the group - +// so the test needs a tolerance, and it is used to REFUSE, so the tolerance has to sit above what a +// correct answer reaches rather than below what a wrong one does. Measured over 113 corpus runs, every +// group determined from its own cell scores under 0.032 (the worst is a C222 call on a cell whose +// alpha is 87.75); the known axis-permuted case scores 0.764. 0.1 leaves 3x headroom over the first +// and 7.6x under the second. +// SearchSpaceGroup's CellHostsRotations asks the same question at 2e-3, normalised per element pair: +// that is the right bound for ADMITTING an extra candidate setting, and it would refuse 6 of those +// 113 runs if it were used here. +constexpr double MAX_METRIC_VIOLATION = 0.1; + // Callbacks for progress and live results. Methods may be called from worker threads, so an // implementation must be thread-safe. The default no-ops suit the CLIs. class RugnuxObserver { diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index df3031a6e..789a8cbf6 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1557,6 +1557,30 @@ static int RunRugnux(int argc, char **argv) { } } + // The cell has to carry the group's own symmetry axes, not merely its centring. --mode scale + // takes both from the file unless -S / -C override them, and an override that puts a 4-fold on + // an axis the crystal does not have merges the data through operators that do not act on its + // indices: measured, -S P41212 on a run indexed a=37.909 b=78.031 c=77.594 reported + // COMPLETENESS= 195.3, since the possible-reflection count is then computed in one frame and + // the observations in another. Nothing downstream can notice, so the check belongs here, with + // both the cell and the group in hand. + if (const auto merge_sg = experiment.GetGemmiSpaceGroup()) { + const auto &uc = experiment.GetUnitCell().value(); + const double violation = MetricViolation(uc, *merge_sg); + if (violation > MAX_METRIC_VIOLATION) { + logger.Error("The space group {} cannot describe this cell: a={:.3f} b={:.3f} c={:.3f} " + "alpha={:.2f} beta={:.2f} gamma={:.2f} is not left invariant by the group's " + "own operators (metric violation {:.3f}, tolerance {:.3f}), so its symmetry " + "axes do not run along the axes the group puts them on. Merging would count " + "the possible reflections in one frame and the observed ones in another - a " + "completeness above 100% is the usual symptom. Merge in a group this cell " + "has, or reprocess the images with --mode mx.", + merge_sg->xhm(), uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, + violation, MAX_METRIC_VIOLATION); + exit(EXIT_FAILURE); + } + } + auto refl_stats = UpdateReflectionResolution(experiment.GetUnitCell().value(), reflections); logger.Info("Read {} reflections from {} images", refl_stats.n_reflections, refl_stats.n_images); experiment.ImagesPerTrigger(refl_stats.n_images);