leonarski_fandClaude Opus 5
b2c0a5ea48
From two independent code reviews of the model-validation work. Nothing here changes a verdict: the acceptance set still reads +17.69 / +1.70 / -0.37 sigma and the rejected runs still write files byte-identical to a run with no model. THE RUN'S OUTPUT. Model validation runs BEFORE the reflection files are written, and two paths through it could throw: the null's replicates (rotated models fed to a scaling path that fails outright on data it cannot pair up - the adversarial input for it), and the mmCIF coordinate writer. Either would have taken the .mtz, .cif, .hkl and _unmerged.mtz with it, after the merge had already been paid for. A null that cannot be built is a question that could not be put, which is the NOT_TESTED state this design already has; a coordinate file that cannot be written is a lost convenience. Both now degrade instead of aborting. THE ORIGIN GAUGE. Translating the whole cell content along a free-origin direction - all three in P1, the unique axis in a polar group - leaves every |F| exactly unchanged. The code said the LM damping and the R-free gate made that harmless between them. Neither does: the gate is a function of |F| and is blind to exactly this, and the gauge column of the Jacobian is not zero but noise divided by the difference step. It is now projected out after every zone, against the group's own common fixed subspace. P2_1 alone is a large share of deposited structures, and the reported shift was partly fiction in every one of them. TWO CLAIMS THAT WERE FALSE. The report told the user R-work carries the decision "because nothing was refined against it", six lines from where six placement parameters are refined against it; the real argument is that the null is placed the same way, so the optimism is common-mode and cancels. And the constant's own comment quoted a +4 vs +1 sigma gap where the measurement is +17.7 vs +1.7. Also: the map file's phase columns are back to [0, 360), the convention they carried before sigma_A weighting; with no data space group the map file follows the reflections into P1 rather than taking the model's group; the sigma is floored against a near-zero null spread rather than only an exactly-zero one; the model is restored whether or not the solver reported a usable answer; and the zone list is the ladder walked rather than the ladder planned. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
762 lines
52 KiB
C++
762 lines
52 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <algorithm>
|
|
#include <fstream>
|
|
#include <sstream>
|
|
|
|
#include <spdlog/fmt/fmt.h>
|
|
|
|
#include "../common/JFJochMath.h" // PI - M_PI is not portable to MSVC
|
|
#include "../common/GitInfo.h"
|
|
#include "../common/time_utc.h"
|
|
#include "../image_analysis/scale_merge/AnisotropyAnalysis.h"
|
|
#include "../image_analysis/scale_merge/Merge.h"
|
|
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
|
|
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
|
|
|
|
#include "ResultReport.h"
|
|
|
|
namespace {
|
|
// The version of this file format. Bumped when a key is renamed or removed, a table column moves,
|
|
// or a reason code changes meaning - a consumer can gate on it.
|
|
constexpr int REPORT_VERSION = 6;
|
|
|
|
const char *BANNER = " ******************************************************************************";
|
|
|
|
void Section(std::ostream &os, const std::string &title) {
|
|
os << "\n" << BANNER << "\n " << title << "\n" << BANNER << "\n\n";
|
|
}
|
|
|
|
// Every number a consumer might want is written as one of these, so it is one grep away.
|
|
template <class T> void Key(std::ostream &os, const char *key, const T &value) {
|
|
os << key << "= " << value << "\n";
|
|
}
|
|
|
|
std::string CellString(const UnitCell &c) {
|
|
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,
|
|
const std::string &input_file,
|
|
const DiffractionExperiment &experiment,
|
|
const ProcessResult &result,
|
|
const RunProvenance &provenance) {
|
|
std::ostringstream os;
|
|
const bool rotation = experiment.IsRotationIndexing();
|
|
const bool merged = result.has_merge_statistics;
|
|
std::vector<std::string> warnings = result.warnings;
|
|
|
|
os << BANNER << "\n"
|
|
<< " RUGNUX PROCESSING REPORT\n"
|
|
<< BANNER << "\n\n"
|
|
<< " What this run determined, written next to its other output. The `KEY= value` lines and\n"
|
|
<< " the tables below are a stable interface - a script greps them, and REPORT_VERSION says\n"
|
|
<< " when that interface last changed. Rates and per-image progress are not here; they are\n"
|
|
<< " on stdout.\n\n";
|
|
|
|
Key(os, "REPORT_VERSION", REPORT_VERSION);
|
|
Key(os, "RUGNUX_VERSION", jfjoch_version());
|
|
if (!jfjoch_git_sha1().empty())
|
|
Key(os, "RUGNUX_GIT", jfjoch_git_sha1().substr(0, 6) + " " + jfjoch_git_date());
|
|
Key(os, "DATE", time_UTC(std::chrono::system_clock::now()));
|
|
Key(os, "INPUT_FILE", input_file);
|
|
Key(os, "OUTPUT_PREFIX", output_prefix);
|
|
// How the result was produced, what it cost and what it ran on, so the report stands on its own
|
|
// once the shell history it came from is gone. Absent rather than zero where the caller does not
|
|
// know them - the library and the viewer have no command line and no invocation to time.
|
|
if (!provenance.command_line.empty())
|
|
Key(os, "COMMAND_LINE", provenance.command_line);
|
|
if (provenance.wall_time_s > 0.0)
|
|
Key(os, "WALL_TIME", fmt::format("{:.2f}", provenance.wall_time_s));
|
|
if (provenance.gpu_count >= 0) {
|
|
Key(os, "GPU_COUNT", provenance.gpu_count);
|
|
// GPU_COUNT= 0 with no GPU= line is the CPU-only case, and saying so is the point: whether
|
|
// the GPUs were there is the first question about how long the run took.
|
|
if (!provenance.gpu_description.empty())
|
|
Key(os, "GPU", provenance.gpu_description);
|
|
}
|
|
|
|
// ---------------------------------------------------------------- 1. DATA SET
|
|
Section(os, "1. DATA SET");
|
|
Key(os, "EXPERIMENT_TYPE", rotation ? "ROTATION" : "STILLS");
|
|
Key(os, "IMAGES_PROCESSED", result.images_processed);
|
|
Key(os, "WAVELENGTH", fmt::format("{:.5f}", experiment.GetWavelength_A()));
|
|
if (const auto gonio = experiment.GetGoniometer()) {
|
|
Key(os, "OSCILLATION_RANGE", fmt::format("{:.4f}", gonio->GetIncrement_deg()));
|
|
Key(os, "STARTING_ANGLE", fmt::format("{:.3f}", gonio->GetStart_deg()));
|
|
const auto ax = gonio->GetAxis();
|
|
Key(os, "ROTATION_AXIS", fmt::format("{:.6f} {:.6f} {:.6f}", ax.x, ax.y, ax.z));
|
|
}
|
|
Key(os, "DETECTOR_DISTANCE", fmt::format("{:.3f}", result.used_distance_mm));
|
|
Key(os, "BEAM_CENTRE", fmt::format("{:.2f} {:.2f}", result.used_beam_x_pxl, result.used_beam_y_pxl));
|
|
// The tilt this run integrated at, and the two points it separates. BEAM_CENTRE is the PONI - the
|
|
// foot of the perpendicular from the sample - so on a tilted detector it is NOT where the direct
|
|
// beam lands, and the two were previously indistinguishable in this report because only one of
|
|
// them was printed. Degrees here; the JSON below carries radians, as the API spells it.
|
|
{
|
|
Key(os, "DETECTOR_TILT", fmt::format("{:.4f} {:.4f} {:.4f}", result.used_detector_tilt_deg[0],
|
|
result.used_detector_tilt_deg[1],
|
|
result.used_detector_tilt_deg[2]));
|
|
Key(os, "DIRECT_BEAM", fmt::format("{:.2f} {:.2f}", result.used_direct_beam_x_pxl,
|
|
result.used_direct_beam_y_pxl));
|
|
if (result.refined_detector_tilt_deg)
|
|
Key(os, "REFINED_DETECTOR_TILT",
|
|
fmt::format("{:.4f} {:.4f}", (*result.refined_detector_tilt_deg)[0],
|
|
(*result.refined_detector_tilt_deg)[1]));
|
|
}
|
|
os << "\n"
|
|
<< " The distance and beam centre above are the ones this result was integrated at, which on\n"
|
|
<< " a rotation run is the post-refined geometry rather than the values in the input file.\n"
|
|
<< " DETECTOR_TILT is rot1/rot2/rot3 in degrees, as the run INTEGRATED at them. Rotation\n"
|
|
<< " indexing does refine rot1/rot2, but the result is never written back onto the geometry,\n"
|
|
<< " so what it measured is reported separately as REFINED_DETECTOR_TILT.\n"
|
|
<< " BEAM_CENTRE is the PONI and DIRECT_BEAM is where the beam actually lands; they differ by\n"
|
|
<< " distance*tan(tilt)/pixel and are identical only when the tilt is zero. Quote whichever the\n"
|
|
<< " program you are feeding expects, and check which one it means.\n";
|
|
if (result.refined_detector_tilt_deg)
|
|
os << " REFINED_DETECTOR_TILT is rot1/rot2 in degrees as rotation indexing MEASURED them, and is\n"
|
|
<< " NOT what this run integrated at - nothing writes a refined tilt back onto the geometry.\n"
|
|
<< " One crystal does not measure a detector tilt. The tilt is aliased with the beam centre,\n"
|
|
<< " so what a single sweep determines is the DIRECT_BEAM above, not the split between the\n"
|
|
<< " two, and this fit stays close to the tilt it was started from - it moves a fraction of\n"
|
|
<< " the way to the true value, and the fraction varies from crystal to crystal. Comparing\n"
|
|
<< " one run against a powder calibration will therefore disagree with it by far more than\n"
|
|
<< " either number's uncertainty. What is worth comparing is the MEDIAN of this value over\n"
|
|
<< " several crystals collected on the same detector, which does track the calibration well\n"
|
|
<< " enough to show up a placeholder or a stale tilt in the file - it does not replace the\n"
|
|
<< " calibration. Do not feed a single run's value back into the instrument.\n"
|
|
<< " rot3 is omitted because a rotation about the beam is an exact null of this experiment\n"
|
|
<< " and the fit cannot move it.\n";
|
|
|
|
// The same geometry once more, as the object jfjoch_broker takes it in: the four required
|
|
// properties of dataset_settings in broker/jfjoch_api.yaml, spelled the way the API spells them.
|
|
// A run that refined the geometry is usually the best measurement of it anyone has, and without
|
|
// this the only way back into the instrument is to read two numbers off this report by eye and
|
|
// retype them. One line, valid JSON, so a script can lift it with a grep and POST it.
|
|
os << "\n";
|
|
{
|
|
// The rotations belong here whenever they are not zero: dataset_settings carries
|
|
// poni_rot1/2/3_rad, and a block that omits them describes a FLAT detector - a different
|
|
// geometry from the one this run used, silently, on every tilted setup. Omitted when zero
|
|
// because the API's own default is 0.0, so the shorter block means the same thing.
|
|
constexpr double RAD = PI / 180.0;
|
|
const auto &t = result.used_detector_tilt_deg;
|
|
std::string rot;
|
|
if (t[0] != 0.0 || t[1] != 0.0 || t[2] != 0.0)
|
|
rot = fmt::format(R"(, "poni_rot1_rad": {:.6f}, "poni_rot2_rad": {:.6f}, )"
|
|
R"("poni_rot3_rad": {:.6f})",
|
|
t[0] * RAD, t[1] * RAD, t[2] * RAD);
|
|
Key(os, "JFJOCH_DATASET_SETTINGS",
|
|
fmt::format(R"({{"beam_x_pxl": {:.2f}, "beam_y_pxl": {:.2f}, "detector_distance_mm": {:.3f}, )"
|
|
R"("incident_energy_keV": {:.4f}{}}})",
|
|
result.used_beam_x_pxl, result.used_beam_y_pxl, result.used_distance_mm,
|
|
experiment.GetDatasetSettings().GetPhotonEnergy_keV(), rot));
|
|
}
|
|
os << "\n"
|
|
<< " The geometry above as jfjoch_broker's dataset_settings, to carry a refined beam centre and\n"
|
|
<< " distance back to the instrument for the next collection. The PONI rotations ride with it\n"
|
|
<< " when they are non-zero, because without them the block describes a flat detector.\n";
|
|
|
|
if (result.spot_resolution_estimate_A.has_value()) {
|
|
os << "\n";
|
|
Key(os, "SPOT_RESOLUTION_ESTIMATE", fmt::format("{:.2f}", *result.spot_resolution_estimate_A));
|
|
os << "\n"
|
|
<< " How far the merged data are expected to reach, read off the found spots alone - no\n"
|
|
<< " lattice, no integration, no merge. It is a prediction, good to about 0.2 A on the\n"
|
|
<< " rotation data it was calibrated on, and it is not what the run achieved: compare it\n"
|
|
<< " with INCLUDE_RESOLUTION_RANGE in section 5. It is not limited to what this detector\n"
|
|
<< " records: where it reads finer than the high-resolution end of that range, the crystal\n"
|
|
<< " diffracts past the corner and the run is detector-limited.\n";
|
|
}
|
|
|
|
if (result.pass_count > 1) {
|
|
os << "\n";
|
|
Key(os, "PASS", fmt::format("{} of {}", result.pass_number, result.pass_count));
|
|
Key(os, "PASS_DECISION", result.pass_decision);
|
|
os << "\n"
|
|
<< " A rotation run integrates twice: once at the geometry in the input file, then again at\n"
|
|
<< " the post-refined geometry. Every number in this report describes the pass named above,\n"
|
|
<< " whose files are " << output_prefix << ".*; the header-geometry pass is not written.\n";
|
|
}
|
|
|
|
// ---------------------------------------------------------------- 2. INDEXING
|
|
Section(os, "2. INDEXING");
|
|
if (result.indexing_rate.has_value())
|
|
Key(os, "INDEXING_RATE", fmt::format("{:.4f}", result.indexing_rate.value()));
|
|
Key(os, "LATTICE_FOUND", (result.consensus_cell.has_value() ? "TRUE" : "FALSE"));
|
|
if (result.consensus_cell.has_value())
|
|
Key(os, "UNIT_CELL_CONSTANTS", CellString(*result.consensus_cell));
|
|
if (result.space_group.has_value()) {
|
|
Key(os, "SPACE_GROUP_NUMBER", result.space_group->number);
|
|
// The number alone does not name the setting, and the search can now return one that is not
|
|
// gemmi's reference setting - P 1 1 2_1 and P 1 2_1 1 are both number 4, on different axes.
|
|
Key(os, "SPACE_GROUP_NAME", result.space_group->xhm());
|
|
}
|
|
if (result.indexing_rate.value_or(0.0f) <= 0.0f)
|
|
warnings.emplace_back("No image indexed - no crystal lattice was determined from this dataset");
|
|
|
|
// ---------------------------------------------- 3. GEOMETRY POST-REFINEMENT
|
|
if (result.post_refine.has_value()) {
|
|
const auto &pr = *result.post_refine;
|
|
Section(os, "3. GEOMETRY POST-REFINEMENT");
|
|
os << " The rotation two-pass fits the detector distance and beam centre from the observed spot\n"
|
|
<< " positions, and the cell scale and rotation axis from the observed rocking angles. Each\n"
|
|
<< " step is committed only if it improves a held-out residual.\n\n";
|
|
Key(os, "POSTREFINE_EVENTS_USED", pr.events_used);
|
|
Key(os, "POSTREFINE_OBS_USED", pr.obs_used);
|
|
Key(os, "POSTREFINE_CELL_COMMITTED", pr.cell_refined ? "TRUE" : "FALSE");
|
|
Key(os, "POSTREFINE_DETECTOR_COMMITTED", pr.detector_refined ? "TRUE" : "FALSE");
|
|
Key(os, "POSTREFINE_DISTANCE", fmt::format("{:.3f} -> {:.3f}", pr.distance_before_mm,
|
|
pr.distance_after_mm));
|
|
Key(os, "POSTREFINE_BEAM_CENTRE", fmt::format("{:.2f} {:.2f} -> {:.2f} {:.2f}",
|
|
pr.beam_x_before_px, pr.beam_y_before_px,
|
|
pr.beam_x_after_px, pr.beam_y_after_px));
|
|
Key(os, "GONIOMETER_ROTATION_SCALE", fmt::format("{:.5f}", pr.rotation_scale));
|
|
Key(os, "GONIOMETER_ROTATION_SCALE_SUSPECT", pr.rotation_scale_suspect ? "TRUE" : "FALSE");
|
|
os << "\n GONIOMETER_ROTATION_SCALE is the factor by which the stage actually turned relative to\n"
|
|
<< " the angles stored in the file (which are the commanded ones). 1.0 = they agree. It drives\n"
|
|
<< " the second integration pass only when SUSPECT is TRUE - both cross-validated and outside\n"
|
|
<< " the tolerance - since a stage that is in fact well calibrated must be left alone. A\n"
|
|
<< " manual --rotation-scale replaces it and is applied to both passes.\n";
|
|
if (pr.rotation_scale_suspect)
|
|
warnings.emplace_back(fmt::format(
|
|
"The goniometer turned by a factor {:.5f} of the angles stored in the file - the "
|
|
"stage rotation looks mis-calibrated by {:+.2f}%. The correction was applied to this "
|
|
"run, but the fault is in the hardware and should be fixed there",
|
|
pr.rotation_scale, 100.0 * (pr.rotation_scale - 1.0)));
|
|
}
|
|
|
|
// ---------------------------------------------- 4. SPACE GROUP DETERMINATION
|
|
Section(os, "4. SPACE GROUP DETERMINATION");
|
|
|
|
// What the data could NOT decide, as keys rather than only as prose. SPACE_GROUP_NAME above is a
|
|
// scalar and reads like a determination; where several groups predict the same absences it is one
|
|
// of them, chosen by convention. A script that greps only the name records a coin-flip as an
|
|
// answer - measured, refining against the deposited model in the wrong enantiomorph gives
|
|
// R = 0.549 - so the ambiguity has to survive the same grep.
|
|
{
|
|
std::string alts;
|
|
if (result.space_group_search.has_value() && result.space_group.has_value())
|
|
for (const auto &alt : result.space_group_search->alternatives)
|
|
if (alt.number != result.space_group->number)
|
|
alts += (alts.empty() ? "" : " | ") + alt.xhm();
|
|
Key(os, "SPACE_GROUP_ALTERNATIVES", alts.empty() ? "NONE" : alts);
|
|
|
|
// The hand is a separate question from the alternatives list, and it is decidable from the
|
|
// 22 groups' numbers alone, so it is answered even where the search did not run. Merged
|
|
// intensities never decide it: an enantiomorphic pair has the same absences and the same
|
|
// Laue class, and only a model, a substructure or an anomalous signal names the hand.
|
|
const char *enantiomorph = "NOT_APPLICABLE";
|
|
if (result.space_group.has_value() && result.space_group->is_enantiomorphic()) {
|
|
// ASSUMED, not determined: nothing here measured the hand. A model that fits carries prior
|
|
// chemical knowledge these intensities do not - and cannot, since |Fcalc| is invariant
|
|
// under the change of hand - so taking its group is an assertion, not a measurement. It is
|
|
// only made where the model was shown to fit (section 10), and it names the hand either by
|
|
// the group actually written being the model's or by the model's having been adopted for it.
|
|
const auto &mv = result.model_validation;
|
|
if (mv.has_value() && mv->ok && mv->model_fits
|
|
&& (mv->adopted_model_enantiomorph
|
|
|| mv->model_space_group_number == result.space_group->number))
|
|
enantiomorph = "ASSUMED_FROM_MODEL";
|
|
else if (!result.space_group_search.has_value())
|
|
enantiomorph = "GIVEN"; // -S, or a reference MTZ: the user's assertion
|
|
else
|
|
enantiomorph = "UNDETERMINED";
|
|
}
|
|
Key(os, "SPACE_GROUP_ENANTIOMORPH", enantiomorph);
|
|
|
|
// A higher point group whose operators the intensities confirmed and whose promotion was
|
|
// refused. It exists only as prose in the search text below, so nothing can act on it - and
|
|
// a refusal is exactly the case where a user might want to try the higher group as well.
|
|
const bool refused = result.space_group_search.has_value()
|
|
&& !result.space_group_search->refused_point_group_hm.empty();
|
|
Key(os, "SPACE_GROUP_REFUSED_POINT_GROUP",
|
|
refused ? result.space_group_search->refused_point_group_hm : std::string("NONE"));
|
|
if (refused)
|
|
Key(os, "SPACE_GROUP_REFUSED_REASON", result.space_group_search->refused_reason);
|
|
os << "\n SPACE_GROUP_ALTERNATIVES names every group these data cannot separate from the one\n"
|
|
<< " adopted (enantiomorphic partners, origin-ambiguous pairs, groups a gap in the data\n"
|
|
<< " leaves untested); NONE means the absences single the answer out. SPACE_GROUP_ENANTIOMORPH\n"
|
|
<< " is UNDETERMINED whenever the group is one of the 22 that come in enantiomorphic pairs and\n"
|
|
<< " nothing outside the merged intensities named the hand - the intensities cannot, so this\n"
|
|
<< " is the normal outcome, not a failure. SPACE_GROUP_REFUSED_POINT_GROUP is a higher point\n"
|
|
<< " group the operator correlations supported and the consistency tests would not take.\n";
|
|
}
|
|
|
|
if (result.space_group_search.has_value()) {
|
|
Key(os, "SPACE_GROUP_SEARCH", "DE_NOVO");
|
|
os << "\n" << SearchSpaceGroupResultToText(*result.space_group_search) << "\n";
|
|
// A centering the data could not test must not read like one they confirmed. The group may
|
|
// still be right - the lattice metric says so - but nothing in these intensities backs it,
|
|
// and that belongs beside the warnings rather than in a table column alone.
|
|
const auto &search = *result.space_group_search;
|
|
if (search.best_space_group.has_value())
|
|
for (const auto &c : search.candidates)
|
|
if (c.space_group.number == search.best_space_group->number && c.centering_untested)
|
|
warnings.emplace_back(fmt::format(
|
|
"The {} centering of {} was NOT confirmed from these data: the crystal was "
|
|
"indexed and integrated on the primitive sub-cell, so the reflections a "
|
|
"{}-centred lattice extinguishes are not in this merge at all. It comes "
|
|
"from the lattice metric. The point group is confirmed from the "
|
|
"intensities; the centering is not",
|
|
search.best_space_group->centring_type(),
|
|
search.best_space_group->short_name(),
|
|
search.best_space_group->centring_type()));
|
|
} else if (result.space_group.has_value()) {
|
|
Key(os, "SPACE_GROUP_SEARCH", "FIXED");
|
|
os << "\n The space group was given, not determined here.\n";
|
|
} else {
|
|
Key(os, "SPACE_GROUP_SEARCH", "NONE");
|
|
os << "\n No space group was determined.\n";
|
|
}
|
|
|
|
// ---------------------------------------------------- 5. SCALING AND MERGING
|
|
Section(os, "5. SCALING AND MERGING");
|
|
if (!merged) {
|
|
Key(os, "MERGE", "NOT_PERFORMED");
|
|
os << "\n No scaling or merging was performed on this run, so there are no merging statistics, no\n"
|
|
<< " error model, and no sweep-quality diagnosis below. The integrated reflections are in\n"
|
|
<< " " << output_prefix << "_process.h5.\n";
|
|
} else {
|
|
const auto &o = result.merge_statistics.overall;
|
|
Key(os, "MERGE", "PERFORMED");
|
|
Key(os, "INCLUDE_RESOLUTION_RANGE", fmt::format("{:.3f} {:.3f}", o.d_max, o.d_min));
|
|
if (result.resolution_fit_A)
|
|
Key(os, "FITTED_RESOLUTION", fmt::format("{:.2f}", *result.resolution_fit_A));
|
|
Key(os, "FRIEDELS_LAW", experiment.GetScalingSettings().GetMergeFriedel() ? "TRUE" : "FALSE");
|
|
Key(os, "UNIQUE_REFLECTIONS", o.unique_reflections);
|
|
Key(os, "TOTAL_OBSERVATIONS", o.total_observations);
|
|
Key(os, "COMPLETENESS", o.possible_unique_reflections > 0
|
|
? fmt::format("{:.1f}", 100.0 * o.unique_reflections / o.possible_unique_reflections)
|
|
: std::string("nan"));
|
|
Key(os, "MULTIPLICITY", o.unique_reflections > 0
|
|
? fmt::format("{:.2f}", static_cast<double>(o.total_observations) / o.unique_reflections)
|
|
: std::string("nan"));
|
|
Key(os, "I_OVER_SIGMA", fmt::format("{:.2f}", o.mean_i_over_sigma));
|
|
Key(os, "R_MEAS", fmt::format("{:.4f}", o.r_meas));
|
|
Key(os, "CC_HALF", fmt::format("{:.4f}", o.cc_half));
|
|
Key(os, "SIGANO", fmt::format("{:.3f}", o.abs_diff_over_sigma_anomalous));
|
|
Key(os, "WILSON_B", fmt::format("{:.2f}", result.merge_statistics.wilson_b));
|
|
// The error model in XDS's convention, so the numbers are directly comparable with a CORRECT.LP.
|
|
Key(os, "ERROR_MODEL_A", fmt::format("{:.4f}", result.error_model_a));
|
|
Key(os, "ERROR_MODEL_B", fmt::format("{:.4e}", result.error_model_b));
|
|
Key(os, "ISA", fmt::format("{:.2f}", result.error_model_isa));
|
|
if (result.error_model_isa_asymptotic > 0.0)
|
|
Key(os, "ISA_ASYMPTOTIC", fmt::format("{:.2f}", result.error_model_isa_asymptotic));
|
|
Key(os, "REFERENCE_DATA_USED", result.has_reference ? "TRUE" : "FALSE");
|
|
// The shell table straight off the statistics rather than result.merge_statistics_text: that
|
|
// string also carries the twinning analysis and the advisories, which have sections of their own.
|
|
os << "\n";
|
|
if (result.resolution_fit_A)
|
|
os << fmt::format(
|
|
" INCLUDE_RESOLUTION_RANGE is the range the reflections were WRITTEN to;"
|
|
" FITTED_RESOLUTION\n is where the CC1/2 fall-off crosses {:.2f}, and is the number to"
|
|
" quote. The data are kept one\n shell past it on purpose: a shell that is included can"
|
|
" still be downweighted or dropped by\n refinement, while one that was truncated cannot"
|
|
" be put back.\n\n",
|
|
experiment.GetScalingSettings().GetResolutionCCTarget());
|
|
os << " ERROR_MODEL_A / ERROR_MODEL_B are in XDS's convention, sigma^2 = a*(sigma0^2 + b*I^2),\n"
|
|
<< " so ISA = 1/sqrt(a*b) means what CORRECT.LP's ISa means. ISA_ASYMPTOTIC, where present,\n"
|
|
<< " is the strong-reflection tier only.\n\n"
|
|
<< result.merge_statistics;
|
|
}
|
|
|
|
// --------------------------------------------------------------- 6. TWINNING
|
|
if (merged && result.twinning.l_test_pairs > 0) {
|
|
Section(os, "6. TWINNING");
|
|
Key(os, "TWINNING_SUSPECTED", result.twinning.twinning_suspected ? "TRUE" : "FALSE");
|
|
Key(os, "L_TEST_MEAN_ABS_L", fmt::format("{:.4f}", result.twinning.mean_abs_l));
|
|
Key(os, "L_TEST_MEAN_L_SQUARED", fmt::format("{:.4f}", result.twinning.mean_l_squared));
|
|
Key(os, "SECOND_MOMENT_I", fmt::format("{:.4f}", result.twinning.second_moment));
|
|
// The same two statistics measured before the group was chosen. Printed next to the ones above
|
|
// rather than instead of them, because they answer different questions and can disagree: these
|
|
// are the only ones a promotion cannot have contaminated.
|
|
if (result.pre_promotion_twinning) {
|
|
Key(os, "L_TEST_MEAN_ABS_L_BEFORE_SEARCH",
|
|
fmt::format("{:.4f}", result.pre_promotion_twinning->mean_abs_l));
|
|
Key(os, "SECOND_MOMENT_I_BEFORE_SEARCH",
|
|
fmt::format("{:.4f}", result.pre_promotion_twinning->second_moment));
|
|
}
|
|
Key(os, "ESTIMATED_TWIN_FRACTION", fmt::format("{:.3f}", result.twinning.estimated_twin_fraction));
|
|
os << "\n" << TwinningAnalysisToText(result.twinning) << "\n";
|
|
if (result.pre_promotion_twinning)
|
|
os << " The _BEFORE_SEARCH pair was measured on the merge the space-group search was given,\n"
|
|
<< " before any point group was adopted. That ordering is the whole point: the statistics\n"
|
|
<< " above are computed in the Laue class this run ADOPTED, so if the search promoted the\n"
|
|
<< " point group they can only report that no twin law exists inside the class it chose -\n"
|
|
<< " and a twin is precisely what would have caused that promotion. Where the two pairs\n"
|
|
<< " disagree, believe the _BEFORE_SEARCH one about whether the crystal is twinned.\n";
|
|
if (result.twinning.twinning_suspected)
|
|
warnings.emplace_back(fmt::format(
|
|
"Twinning is indicated (<|L|> = {:.3f}, <I^2>/<I>^2 = {:.3f}, estimated twin "
|
|
"fraction {:.2f}) - refine against the merged data with care",
|
|
result.twinning.mean_abs_l, result.twinning.second_moment,
|
|
result.twinning.estimated_twin_fraction));
|
|
}
|
|
|
|
// ------------------------------------------------------- 7. RADIATION DAMAGE
|
|
if (!result.radiation_damage_text.empty()) {
|
|
Section(os, "7. RADIATION DAMAGE");
|
|
// A number, or a word saying why there is none: NOT_A_TREND where the per-batch curve was measured
|
|
// but no straight line describes it (damage is progressive, so that curve is not dose), NOT_MEASURED
|
|
// where the monitor could not run at all.
|
|
const double db = result.merge_statistics.radiation_damage_delta_b;
|
|
Key(os, "RADIATION_DAMAGE_RELATIVE_B",
|
|
std::isfinite(db) ? fmt::format("{:.2f}", db)
|
|
: result.merge_statistics.radiation_damage_b_batch.empty() ? std::string("NOT_MEASURED")
|
|
: std::string("NOT_A_TREND"));
|
|
os << "\n" << result.radiation_damage_text << "\n";
|
|
}
|
|
|
|
// ------------------------------------------------------ 8. SWEEP QUALITY
|
|
const auto &sq = result.merge_statistics.sweep_quality;
|
|
Section(os, "8. SWEEP QUALITY");
|
|
os << " Stretches of the sweep over which the crystal delivered much less than the rest of the run.\n"
|
|
<< " REASON comes from a closed vocabulary, listed below so a consumer can tell an unknown code\n"
|
|
<< " from a missing one. SEVERITY is the fraction of the run's typical diffracting power missing\n"
|
|
<< " over the range (0 = as good as the run, 1 = nothing at all); SCALE and CC are the range's\n"
|
|
<< " mean per-image scale and CC-to-merge relative to the run median; INDEXED is the fraction of\n"
|
|
<< " the range's frames that were scaled at all. Nothing is excluded on the strength of this.\n\n";
|
|
Key(os, "SWEEP_QUALITY_STATUS", sq.measured ? "COMPUTED" : "NOT_COMPUTED");
|
|
Key(os, "SWEEP_QUALITY_COUNT", sq.ranges.size());
|
|
{
|
|
std::string codes;
|
|
for (int r = 0; r <= static_cast<int>(SweepQualityReason::RadiationDamage); ++r)
|
|
codes += (codes.empty() ? "" : " ")
|
|
+ std::string(SweepQualityReasonCode(static_cast<SweepQualityReason>(r)));
|
|
Key(os, "SWEEP_QUALITY_REASONS", codes);
|
|
}
|
|
if (sq.measured) {
|
|
Key(os, "SWEEP_ROTATION", fmt::format("{:.1f}", sq.sweep_deg));
|
|
Key(os, "FLUX_PEAK_TO_TROUGH", fmt::format("{:.2f}", sq.flux_peak_to_trough));
|
|
Key(os, "SCALE_MODULATION_PEAK_TO_TROUGH", fmt::format("{:.2f}", sq.modulation_peak_to_trough));
|
|
}
|
|
os << "\n"
|
|
<< " FIRST_IMAGE LAST_IMAGE N_IMAGES ROTATION REASON SEVERITY SCALE CC INDEXED\n"
|
|
<< " ----------- ----------- --------- -------- -------------------- -------- ------ ------ --------\n";
|
|
for (const auto &r : sq.ranges) {
|
|
os << fmt::format(" {:11d} {:11d} {:9d} {:8.1f} {:<20} {:8.2f} {:6.2f} {:6.2f} {:8.2f}\n",
|
|
r.first_image, r.last_image, r.last_image - r.first_image + 1, r.rotation_deg,
|
|
SweepQualityReasonCode(r.reason), r.severity, r.mean_relative_scale,
|
|
r.mean_relative_cc, r.indexed_fraction);
|
|
warnings.push_back(fmt::format(
|
|
"Frames {}-{} {} ({:.1f} deg, scale {:.2f} and CC {:.2f} of the run, {:.0f}% scaled)",
|
|
r.first_image, r.last_image, SweepQualityReasonText(r.reason), r.rotation_deg,
|
|
r.mean_relative_scale, r.mean_relative_cc, 100.0 * r.indexed_fraction));
|
|
}
|
|
os << " ----------- ----------- --------- -------- -------------------- -------- ------ ------ --------\n";
|
|
|
|
// ---------------------------------------------------------- 9. DIFFRACTION ANISOTROPY
|
|
const auto &an = result.merge_statistics.anisotropy;
|
|
if (merged && an.n_reflections > 0) {
|
|
Section(os, "9. DIFFRACTION ANISOTROPY");
|
|
os << " How much the fall-off depends on direction, and whether that is established above this\n"
|
|
<< " data set's own systematic error. Nothing here corrects an intensity or removes a\n"
|
|
<< " reflection: the merged data and the written files do not depend on direction at all.\n"
|
|
<< " ANISOTROPY_DELTA_B is the range of the principal components of the anisotropy tensor,\n"
|
|
<< " on the ordinary crystallographic B scale (the same scale as phenix.xtriage's B_cart and\n"
|
|
<< " ctruncate's anisotropic B), fitted on intensities with nothing dropped;\n"
|
|
<< " ANISOTROPY_SIGNIFICANCE gates ANISOTROPY_DELTA_B_LINEAR, the part of it that follows\n"
|
|
<< " exp(-1/2 s^T B s), which is not the same number. A 1 in ANISOTROPY_D_MIN_CENSORED marks\n"
|
|
<< " a direction whose limit is the edge of the measured data rather than the crystal's own.\n\n";
|
|
Key(os, "ANISOTROPY_VERDICT", AnisotropyVerdictCode(an.verdict));
|
|
Key(os, "ANISOTROPY_FREE_DIRECTIONS", an.n_free_parameters);
|
|
Key(os, "ANISOTROPY_DELTA_B", fmt::format("{:.2f}", an.delta_b));
|
|
Key(os, "ANISOTROPY_DELTA_B_LINEAR", fmt::format("{:.2f}", an.delta_b_linear));
|
|
Key(os, "ANISOTROPY_PRINCIPAL_B", fmt::format("{:.2f} {:.2f} {:.2f}",
|
|
an.eigenvalue[0] - an.eigenvalue[2],
|
|
an.eigenvalue[1] - an.eigenvalue[2], 0.0));
|
|
Key(os, "ANISOTROPY_FOLD_WEAKENING", fmt::format("{:.1f}", an.fold_weakening));
|
|
Key(os, "ANISOTROPY_D_MIN_PRINCIPAL", fmt::format("{:.2f} {:.2f} {:.2f}", an.d_min_axis[0],
|
|
an.d_min_axis[1], an.d_min_axis[2]));
|
|
Key(os, "ANISOTROPY_D_MIN_CENSORED", fmt::format("{} {} {}", an.d_min_censored[0] ? 1 : 0,
|
|
an.d_min_censored[1] ? 1 : 0,
|
|
an.d_min_censored[2] ? 1 : 0));
|
|
Key(os, "ANISOTROPY_D_MIN_SPREAD", fmt::format("{:.2f}", an.d_min_spread));
|
|
// The finest of the three principal limits, on its own line so it can be grepped. The cut this
|
|
// run applied is isotropic, so on an anisotropic crystal the outer shells are complete in count
|
|
// and empty in signal along the weak directions; this says how far the crystal actually reaches
|
|
// where it reaches furthest.
|
|
Key(os, "ANISOTROPY_D_MIN_BEST",
|
|
fmt::format("{:.2f}", *std::min_element(an.d_min_axis, an.d_min_axis + 3)));
|
|
Key(os, "ANISOTROPY_SHAPE", AnisotropyShapeCode(an.shape));
|
|
Key(os, "ANISOTROPY_SHAPE_INTERCEPT", fmt::format("{:.3f}", an.shape_intercept));
|
|
Key(os, "ANISOTROPY_SHAPE_INTERCEPT_Z", fmt::format("{:.1f}", an.shape_intercept_z));
|
|
Key(os, "ANISOTROPY_SHAPE_SLOPE", fmt::format("{:.2f}", an.shape_slope));
|
|
Key(os, "ANISOTROPY_SHAPE_RESIDUAL", fmt::format("{:.1f}", an.shape_residual));
|
|
Key(os, "ANISOTROPY_N_OBSERVATIONS", an.n_observations);
|
|
Key(os, "ANISOTROPY_SIGMA_SYSTEMATIC", fmt::format("{:.3f}", an.sigma_systematic));
|
|
Key(os, "ANISOTROPY_FORBIDDEN_Z", fmt::format("{:.1f}", an.forbidden_z));
|
|
Key(os, "ANISOTROPY_FLOOR", fmt::format("{:.3f}", an.floor));
|
|
Key(os, "ANISOTROPY_SIGNIFICANCE", fmt::format("{:.2f}", an.significance));
|
|
Key(os, "ANISOTROPY_DETECTION_LIMIT", fmt::format("{:.2f}", an.detection_limit));
|
|
os << "\n" << AnisotropyToText(an) << "\n";
|
|
// Three different resolutions appear in this report and they answer three different questions.
|
|
// Say so here rather than let a reader assume one of them is "the" resolution.
|
|
if (result.resolution_fit_A)
|
|
os << fmt::format(
|
|
" Three resolutions, three questions. INCLUDE_RESOLUTION_RANGE ({:.2f} A) is what was\n"
|
|
" WRITTEN, deliberately one shell past the fit. FITTED_RESOLUTION ({:.2f} A) is where\n"
|
|
" the isotropic CC1/2 fall-off crosses its target - the overall measurement.\n"
|
|
" ANISOTROPY_D_MIN_BEST ({:.2f} A) is how far the crystal reaches along its strongest\n"
|
|
" direction. It is read off the fitted tensor, not off the cut, so it can land either\n"
|
|
" side of the other two - and where ANISOTROPY_D_MIN_CENSORED flags that direction it\n"
|
|
" is the edge of the measured data rather than the crystal's limit. None of the three\n"
|
|
" is wrong; quoting one without saying which it is, is.\n\n",
|
|
result.merge_statistics.overall.d_min, *result.resolution_fit_A,
|
|
*std::min_element(an.d_min_axis, an.d_min_axis + 3));
|
|
// CC1/2 that falls and then climbs again is not a fall-off, so the single number the cut is read
|
|
// off does not describe these data. Reported as an observation, with no advice to re-cut: the
|
|
// generous cut is deliberate and refinement can downweight what it does not want.
|
|
{
|
|
// A climb only counts while the shell still carries signal: past CC1/2 ~ 0.3 the data are
|
|
// beyond any cut this run would take and the number is oscillating in noise, which says
|
|
// nothing about the fall-off. 0.05 is well above that surviving noise.
|
|
constexpr double CC_HALF_RISE = 0.05, CC_HALF_ALIVE = 0.30;
|
|
const auto &sh = result.merge_statistics.shells;
|
|
size_t fell = 0, rose_after = 0;
|
|
for (size_t i = 1; i < sh.size(); ++i) {
|
|
if (sh[i].cc_half < sh[i - 1].cc_half) fell = i;
|
|
else if (fell > 0 && sh[i].cc_half > sh[i - 1].cc_half + CC_HALF_RISE
|
|
&& sh[i].cc_half > CC_HALF_ALIVE) rose_after = i;
|
|
}
|
|
if (rose_after > 0)
|
|
warnings.emplace_back(fmt::format(
|
|
"CC1/2 is not monotone with resolution (it climbs again at {:.2f}-{:.2f} A) - the "
|
|
"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) {
|
|
// 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<int>(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 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
|
|
if (result.model_validation.has_value()) {
|
|
const auto &mv = *result.model_validation;
|
|
Section(os, "10. MODEL VALIDATION");
|
|
Key(os, "MODEL_FILE", mv.model_path);
|
|
if (!mv.ok) {
|
|
Key(os, "MODEL_VALIDATION", "NOT_PERFORMED");
|
|
Key(os, "MODEL_VALIDATION_REASON", mv.failure_reason);
|
|
os << "\n A model was given but could not be used, so there are no R-factors and no maps.\n"
|
|
<< " Everything else in this report is unaffected: the merge does not depend on the model.\n";
|
|
} else {
|
|
// The counterpart of NOT_PERFORMED above, so a consumer can grep one key for either answer
|
|
// instead of having to infer success from the absence of a failure line.
|
|
Key(os, "MODEL_VALIDATION", "PERFORMED");
|
|
Key(os, "MODEL_SPACE_GROUP_NUMBER", mv.model_space_group_number);
|
|
Key(os, "R_WORK", fmt::format("{:.4f}", mv.r_work));
|
|
Key(os, "R_FREE", fmt::format("{:.4f}", mv.r_free));
|
|
Key(os, "R_WORK_REFLECTIONS", mv.n_work);
|
|
Key(os, "R_FREE_REFLECTIONS", mv.n_free);
|
|
Key(os, "BULK_SOLVENT_K_SOL", fmt::format("{:.3f}", mv.k_sol));
|
|
Key(os, "BULK_SOLVENT_B_SOL", fmt::format("{:.1f}", mv.b_sol));
|
|
Key(os, "SCALE_OVERALL", fmt::format("{:.4f}", mv.k_overall));
|
|
// Does this model describe these data? Not a threshold on R - the R a model that explains
|
|
// nothing reaches depends on the model as much as on the data - but the same model refitted
|
|
// from random orientations, which is the only null that fits both. NOT_TESTED is a third
|
|
// answer and not a missing one: the model claimed nothing that needed arbitrating, so the
|
|
// null was not built, and a consumer must not read that as a model the data refused.
|
|
Key(os, "MODEL_FIT",
|
|
!mv.fit_tested ? "NOT_TESTED" : (mv.model_fits ? "ACCEPTED" : "REJECTED"));
|
|
if (mv.fit_tested) {
|
|
Key(os, "MODEL_FIT_STATISTIC", "R_WORK");
|
|
Key(os, "MODEL_FIT_VALUE", fmt::format("{:.4f}", mv.r_work));
|
|
Key(os, "MODEL_FIT_NULL_MEAN", fmt::format("{:.4f}", mv.null_r_work_mean));
|
|
Key(os, "MODEL_FIT_NULL_SD", fmt::format("{:.4f}", mv.null_r_work_sd));
|
|
Key(os, "MODEL_FIT_NULL_REPLICATES", mv.null_replicates);
|
|
Key(os, "MODEL_FIT_SIGMA", fmt::format("{:+.2f}", mv.r_work_sigma));
|
|
}
|
|
// The model is placed against the data as one rigid body, and only where the free
|
|
// reflections say it helped - so R_FREE_BEFORE_RIGID_BODY says what the placement bought.
|
|
Key(os, "RIGID_BODY", mv.rigid_body_applied ? "APPLIED" : "NOT_APPLIED");
|
|
if (mv.rigid_body_applied) {
|
|
Key(os, "RIGID_BODY_ROTATION_DEG", fmt::format("{:.3f}", mv.rigid_body_angle_deg));
|
|
Key(os, "RIGID_BODY_SHIFT_A", fmt::format("{:.3f}", mv.rigid_body_shift_A));
|
|
Key(os, "R_FREE_BEFORE_RIGID_BODY", fmt::format("{:.4f}", mv.r_free_before_rigid_body));
|
|
}
|
|
Key(os, "MAP_COEFFICIENTS", "2mFo-DFc / mFo-DFc");
|
|
Key(os, "MAP_SIGMA_A_SHELLS", mv.sigma_a_shells);
|
|
Key(os, "MAP_MEAN_FOM", fmt::format("{:.3f}", mv.mean_fom));
|
|
Key(os, "MEAN_ATOM_DENSITY_SIGMA", fmt::format("{:.2f}", mv.mean_atom_density_sigma));
|
|
// The anomalous scatterers the data themselves found, named by the model's atoms.
|
|
if (!mv.anomalous_sites.empty()) {
|
|
Key(os, "ANOMALOUS_BIJVOET_PAIRS", mv.anomalous_pairs);
|
|
for (size_t i = 0; i < mv.anomalous_sites.size(); i++)
|
|
// Two digits so the ten keys are the same width and the values line up.
|
|
Key(os, fmt::format("ANOMALOUS_SITE_{:02}", i + 1).c_str(),
|
|
fmt::format("{:<18} {:6.2f} sigma", mv.anomalous_sites[i].label,
|
|
mv.anomalous_sites[i].sigma));
|
|
}
|
|
// What was applied to the written reflections, so a reader can tell whether the file is
|
|
// in the indexing it was merged in or in the model's. The enantiomorph is a label only.
|
|
const bool took_indexing = !(mv.indexing_op == gemmi::Op::identity());
|
|
Key(os, "MODEL_DECISIONS_TAKEN",
|
|
mv.adopted_model_enantiomorph
|
|
? (took_indexing ? "ENANTIOMORPH+INDEXING" : "ENANTIOMORPH")
|
|
: (took_indexing ? "INDEXING" : "NONE"));
|
|
Key(os, "MODEL_ENANTIOMORPH_ADOPTED", mv.adopted_model_enantiomorph ? "TRUE" : "FALSE");
|
|
Key(os, "MODEL_INDEXING_OPERATOR", mv.indexing_op.triplet());
|
|
// The merohedral choice beside what the same choice looks like when it is made by a model
|
|
// in a random orientation - which also picks a winner, and by a comparable lead.
|
|
if (mv.indexing_probed) {
|
|
Key(os, "MODEL_INDEXING_MARGIN", fmt::format("{:.4f}", mv.indexing_margin));
|
|
Key(os, "MODEL_INDEXING_MARGIN_NULL",
|
|
fmt::format("{:.4f} +- {:.4f}", mv.indexing_margin_null_mean,
|
|
mv.indexing_margin_null_sd));
|
|
Key(os, "MODEL_INDEXING_MARGIN_SIGMA",
|
|
fmt::format("{:+.2f}", mv.indexing_margin_sigma));
|
|
}
|
|
if (!mv.maps_prefix.empty())
|
|
Key(os, "MAPS_PREFIX", mv.maps_prefix);
|
|
os << "\n R-free here measures the merged intensities against an external structure, which is\n"
|
|
<< " what CC1/2 and R_meas cannot do - they only measure the data against themselves. The\n"
|
|
<< " model is not refined: it is scaled to the data with a flat bulk solvent and an overall\n"
|
|
<< " anisotropic B, so these R-factors are higher than a refined structure's and are a\n"
|
|
<< " data-quality reading, not a refinement result.\n";
|
|
if (mv.fit_tested)
|
|
os << "\n The model is a hypothesis, and MODEL_FIT says whether these data accept it. There is no\n"
|
|
<< " value of R that settles that on its own - what a model which explains nothing reaches\n"
|
|
<< " depends on its atom count and B-factors as much as on the data - so the same model was\n"
|
|
<< " refitted, and re-placed, from " << mv.null_replicates << " random orientations about its own centroid, and\n"
|
|
<< " MODEL_FIT_SIGMA is how far the real fit sits above that null. R-work carries the\n"
|
|
<< " decision because the null was placed the same way the real fit was, so what those six\n"
|
|
<< " placement parameters buy is bought on both sides and cancels - and it has far more\n"
|
|
<< " reflections than R-free. MODEL_DECISIONS_TAKEN names what the model was allowed to change about the\n"
|
|
<< " written reflections; R-factors, maps and the rigid-body placement are reported either\n"
|
|
<< " way, because they describe the model, not the data.\n";
|
|
else
|
|
os << "\n MODEL_FIT= NOT_TESTED, which is not a failed test. A model can change only two things\n"
|
|
<< " about the written reflections - the space-group label, where it asserts the other\n"
|
|
<< " enantiomorph, and the indexing - and this one asserted neither: it names no other\n"
|
|
<< " hand for these data, and prefers the indexing they were merged in. With nothing to\n"
|
|
<< " arbitrate there was nothing to arbitrate it against, so the null the other two answers\n"
|
|
<< " are measured against was not built and the run did not pay for it. The R-factors, the\n"
|
|
<< " maps and the rigid-body placement above are exactly what they would have been.\n";
|
|
if (mv.fit_tested && !mv.model_fits)
|
|
os << "\n The model was tried and REJECTED. It was scaled and placed against these data exactly\n"
|
|
<< " as an accepted one would have been, and reached R-work " << fmt::format("{:.4f}", mv.r_work) << " where the same\n"
|
|
<< " model in " << mv.null_replicates << " random orientations reached "
|
|
<< fmt::format("{:.4f} +- {:.4f}", mv.null_r_work_mean, mv.null_r_work_sd) << " - "
|
|
<< fmt::format("{:+.2f}", mv.r_work_sigma) << " sigma, which is no\n"
|
|
<< " better than chance. Nothing downstream moved: the space group is the one the data\n"
|
|
<< " were merged in, the indexing is the one they were merged in, and the reflection\n"
|
|
<< " files are byte for byte what a run with no model would have written. The R-factors\n"
|
|
<< " and the maps above still describe this model against these data - they are the\n"
|
|
<< " negative result, not a failure of the run.\n";
|
|
if (!mv.anomalous_sites.empty())
|
|
os << "\n The anomalous sites are the highest peaks of the anomalous difference map -\n"
|
|
<< " F(+)-F(-) on the model phase turned back by 90 degrees - read at the model's own\n"
|
|
<< " atom centres, so each one is named rather than left as a coordinate. A dataset with\n"
|
|
<< " no anomalous signal still lists ten sites: it is their height, a few sigma at most,\n"
|
|
<< " that says so. The map itself is written as _anom.ccp4, where a scatterer the model\n"
|
|
<< " does not contain would show up as a peak on nothing.\n";
|
|
if (!(mv.indexing_op == gemmi::Op::identity()))
|
|
os << "\n The written reflections were reindexed into the model's frame - the operator above\n"
|
|
<< " says how - so the reflection files, the R-factors and the maps all describe one\n"
|
|
<< " indexing.\n";
|
|
if (mv.adopted_model_enantiomorph)
|
|
os << "\n The written reflections carry the model's enantiomorph as their space group. That is a\n"
|
|
<< " change of label and nothing else: the two groups of an enantiomorphic pair have the same\n"
|
|
<< " rotation operations, so no reflection moved. Reindexing by the change-of-hand operator\n"
|
|
<< " would have swapped I(+) with I(-) - flipping the anomalous differences, not correcting\n"
|
|
<< " them - on the strength of a label the space-group search reports as undetermined.\n";
|
|
else if (mv.model_enantiomorph_candidate)
|
|
os << "\n The model asserts the other enantiomorph of the written space group, and that assertion\n"
|
|
<< " was NOT taken up. Adopting it would rewrite the group of every reflection written here on\n"
|
|
<< " the strength of a model these data did not accept, and merged intensities cannot check it:\n"
|
|
<< " |Fcalc| is invariant under the change of hand, so R-free reads the same for both. The\n"
|
|
<< " anomalous difference map is the only measurement here that is sensitive to the hand, and\n"
|
|
<< " where it says the two disagree it vetoes the adoption outright.\n";
|
|
if (mv.anomalous_hands_disagree)
|
|
warnings.emplace_back(fmt::format(
|
|
"The anomalous density at the model's atoms is inverted ({} reads {:.1f} sigma, deeper "
|
|
"than the highest peak) - the data and the model are in opposite hands. Either the model "
|
|
"is the wrong enantiomorph for this crystal or the data were indexed in the wrong hand; "
|
|
"the reflections have not been reindexed, which would have hidden which of the two it is",
|
|
mv.anomalous_deepest_site, mv.anomalous_deepest_sigma));
|
|
}
|
|
}
|
|
|
|
// --------------------------------------------------------------- 11. WARNINGS
|
|
if (result.cancelled)
|
|
warnings.emplace_back(fmt::format("Processing was cancelled after {} images - this report "
|
|
"describes an incomplete run", result.images_processed));
|
|
Section(os, "11. WARNINGS");
|
|
os << " Everything that needs a person's attention, one line each, marked so a script can find\n"
|
|
<< " them with a single grep for \"WARNING:\".\n\n";
|
|
Key(os, "WARNING_COUNT", warnings.size());
|
|
os << "\n";
|
|
for (const auto &w : warnings)
|
|
os << "WARNING: " << w << "\n";
|
|
if (warnings.empty())
|
|
os << " (none)\n";
|
|
|
|
os << "\n" << BANNER << "\n END OF REPORT\n" << BANNER << "\n";
|
|
return os.str();
|
|
}
|
|
|
|
void WriteResultReport(const std::string &output_prefix,
|
|
const std::string &input_file,
|
|
const DiffractionExperiment &experiment,
|
|
const ProcessResult &result,
|
|
Logger &logger,
|
|
const RunProvenance &provenance) {
|
|
if (output_prefix.empty())
|
|
return; // "compute the statistics, persist nothing"
|
|
|
|
const std::string filename = output_prefix + "_report.txt";
|
|
// The report is unconditional, so it must never be the reason a run fails: a run that produced a
|
|
// good .mtz must survive an unwritable path or a full disk. Report the failure and carry on.
|
|
try {
|
|
std::ofstream file(filename);
|
|
file.exceptions(std::ios::failbit | std::ios::badbit);
|
|
file << RenderResultReport(output_prefix, input_file, experiment, result, provenance);
|
|
} catch (const std::exception &e) {
|
|
logger.Warning("Could not write the results report {}: {}", filename, e.what());
|
|
}
|
|
}
|