-S: check the metric the fixed group needs, not only its centring
Both guards on a user-fixed space group keyed on the centring differing, so an axis permutation within the same centring passed: a run that indexed a=37.909 b=78.031 c=77.594 and was told -S P41212 - which needs a=b with the 4-fold along c - merged through operators that do not act on its own indices, and --mode scale on it reported COMPLETENESS= 195.3, an arithmetically impossible number, without complaint. MetricViolation asks the setting-independent question instead: a group's rotations must leave the cell's metric tensor invariant. The re-seating arm now runs on either failure, so a permuted cell can be reindexed into the setting the group needs rather than merely refused; the refusal arm catches what re-seating could not fix. --mode scale never reaches either arm, so it gets the same test where it fixes its cell and its group, which is the only place an impossible completeness could still be produced. The tolerance is a refusal bound, so it sits above what a correct answer reaches: over 113 corpus runs every group determined from its own cell scores under 0.032 and the permuted case scores 0.764. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
This commit is contained in:
Claude Opus 5
@@ -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.
|
||||
|
||||
+67
-23
@@ -64,6 +64,34 @@
|
||||
#include <array>
|
||||
#include <map>
|
||||
|
||||
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<double>(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()));
|
||||
}
|
||||
|
||||
@@ -319,6 +319,25 @@ struct ProcessResult {
|
||||
std::vector<std::string> 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 {
|
||||
|
||||
@@ -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);
|
||||
|
||||
Reference in New Issue
Block a user