leonarski_fandClaude Opus 5.5
3908956540
SHELXL refused the .ins rugnux writes beside the .hkl: its REM line was
over 80 characters ("** INPUT INSTRUCTION REM Rugnux 1.0.0-rc.... IS LONGER
THAN 80 CHARACTERS **"), and wherever the run had no cell e.s.d.s the ZERR
line was left out, which SHELXL stops on ("** ZERR MISSING **"). SHELXT read
both versions, so the problem showed only once a SHELXT solution, or the
user's own model, went on to SHELXL.
- The remarks are split into lines within 80 characters, TITL is cut to 80.
- ZERR is always written, after CELL; where the e.s.d.s were not measured
(stills, --mode scale) they are 0 and a REM says so.
- The joint merge of several sweeps measured no e.s.d.s (its fulls carry no
detector positions), so its .ins had no ZERR on every multi-sweep run; it
now takes the reference sweep's, which also reach its mmCIF and report.
- TEMP (Celsius) is written from the sample temperature the data file
records, the same value as the mmCIF's _diffrn.ambient_temp; SHELXT carries
it into its .res and SHELXL into its CIF (_diffrn_ambient_temperature
100(2)). Without one a REM says SHELXL will assume 20 C.
- The NXmx reader takes /entry/sample/temperature in its declared unit: K
(or none) as is, Celsius converted, an unknown unit warned about and not
read. It was read as K whatever the file said.
The .hkl was already right (3I4,2F8.2, |I| and sigma scaled to at most
9999, 0 0 0 terminator, unmerged at the measured index on rotation data so
Friedel mates stay apart) and is unchanged. Checked with SHELXT (solutions
matching the COD structures atom for atom), SHELXL from the SHELXT solution
and from the COD model, SHELXC/SHELXD on two anomalous protein sets, and
SHELXL reading a 1.2 million-record protein .hkl.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UqXYdL4qBGYRiSXSGXfCLu
1483 lines
83 KiB
C++
1483 lines
83 KiB
C++
// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "WriteReflections.h"
|
|
#include "geom_refinement/CellEsd.h"
|
|
#include "scale_merge/Merge.h"
|
|
#include "scale_merge/HKLKey.h"
|
|
#include "scale_merge/TwinningAnalysis.h"
|
|
#include "bragg_integration/SystematicAbsence.h"
|
|
#include "../common/ParallelFor.h"
|
|
#include "../writer/WriteRugnuxFiles.h"
|
|
|
|
#include <algorithm>
|
|
#include <charconv>
|
|
#include <cmath>
|
|
#include <cctype>
|
|
#include <cstring>
|
|
#include <numeric>
|
|
#include <cstdio>
|
|
#include <map>
|
|
#include <memory>
|
|
#include <tuple>
|
|
#include <filesystem>
|
|
#include <fstream>
|
|
#include <iomanip>
|
|
#include <sstream>
|
|
#include <stdexcept>
|
|
#include <ctime>
|
|
#include <chrono>
|
|
|
|
#include <gemmi/mtz.hpp>
|
|
|
|
#include "../common/GitInfo.h"
|
|
|
|
namespace {
|
|
|
|
/// Current date in ISO-8601 (YYYY-MM-DD) for the _audit block.
|
|
std::string CurrentDateISO() {
|
|
auto now = std::chrono::system_clock::now();
|
|
auto t = std::chrono::system_clock::to_time_t(now);
|
|
std::tm tm{};
|
|
#ifdef _WIN32
|
|
gmtime_s(&tm, &t);
|
|
#else
|
|
gmtime_r(&t, &tm);
|
|
#endif
|
|
char buf[32];
|
|
std::strftime(buf, sizeof(buf), "%Y-%m-%d", &tm);
|
|
return buf;
|
|
}
|
|
|
|
/// The C locale's "%.*f" of a finite double. std::to_chars where the standard library has it for
|
|
/// floating point: the same digits - both round the exact binary value, halfway cases to even - at a
|
|
/// quarter of snprintf's time, and with no locale to consult. snprintf elsewhere (libc++ before macOS
|
|
/// 13.3). 512 characters hold any finite double.
|
|
std::string FixedDigits(double val, int decimals) {
|
|
char buf[512];
|
|
#ifdef __cpp_lib_to_chars
|
|
const auto r = std::to_chars(buf, buf + sizeof(buf), val, std::chars_format::fixed, decimals);
|
|
return std::string(buf, r.ptr);
|
|
#else
|
|
const int n = std::snprintf(buf, sizeof(buf), "%.*f", decimals, val);
|
|
return std::string(buf, n);
|
|
#endif
|
|
}
|
|
|
|
/// Format a double with given decimal places; returns "?" for non-finite.
|
|
/// Not an ostringstream: the reflection loop below calls this twelve times per reflection, and
|
|
/// building a stream (and its locale) per call dominated the time spent writing a merged file.
|
|
std::string Fmt(double val, int decimals = 4) {
|
|
if (!std::isfinite(val))
|
|
return "?";
|
|
return FixedDigits(val, decimals);
|
|
}
|
|
|
|
/// Quote a CIF string value; returns "?" for empty.
|
|
std::string CifStr(const std::string& s) {
|
|
if (s.empty())
|
|
return "?";
|
|
// If it contains spaces or special chars, single-quote it
|
|
if (s.find(' ') != std::string::npos ||
|
|
s.find('\'') != std::string::npos ||
|
|
s.find('#') != std::string::npos)
|
|
return "'" + s + "'";
|
|
return s;
|
|
}
|
|
|
|
std::string Upper(std::string s) {
|
|
for (char &c : s)
|
|
c = static_cast<char>(std::toupper(static_cast<unsigned char>(c)));
|
|
return s;
|
|
}
|
|
|
|
// The general class of the detector, _diffrn_detector.detector in the PDBx enumeration, where its make
|
|
// and model name it; empty where they do not.
|
|
std::string DetectorClass(const std::string &description) {
|
|
const std::string d = Upper(description);
|
|
if (d.find("CCD") != std::string::npos)
|
|
return "CCD";
|
|
if (d.find("IMAGE PLATE") != std::string::npos)
|
|
return "IMAGE PLATE";
|
|
for (const char *pixel : {"EIGER", "PILATUS", "JUNGFRAU"})
|
|
if (d.find(pixel) != std::string::npos)
|
|
return "PIXEL";
|
|
return "";
|
|
}
|
|
|
|
// The general class of the source, _diffrn_source.source in the PDBx enumeration, from the NXmx
|
|
// /entry/source/type the file states; empty where it states none or one PDBx has no class for.
|
|
std::string SourceClass(const std::string &nxmx_type) {
|
|
const std::string t = Upper(nxmx_type);
|
|
if (t.find("SYNCHROTRON") != std::string::npos) return "SYNCHROTRON";
|
|
if (t.find("FREE-ELECTRON") != std::string::npos || t.find("FREE ELECTRON") != std::string::npos)
|
|
return "FREE ELECTRON LASER";
|
|
if (t.find("ROTATING ANODE") != std::string::npos) return "ROTATING ANODE";
|
|
if (t.find("FIXED TUBE") != std::string::npos) return "SEALED TUBE";
|
|
if (t.find("METAL JET") != std::string::npos) return "LIQUID ANODE";
|
|
return "";
|
|
}
|
|
|
|
// The date part of an ISO-8601 time stamp (YYYY-MM-DD...), which is all the PDBx collection date
|
|
// takes; empty for anything else.
|
|
std::string IsoDate(const std::string ×tamp) {
|
|
if (timestamp.size() < 10 || timestamp[4] != '-' || timestamp[7] != '-')
|
|
return "";
|
|
return timestamp.substr(0, 10);
|
|
}
|
|
|
|
// The commit this build was made from, with its date where the build knows it ("unknown" for a source
|
|
// tree without git history).
|
|
std::string GitDescription() {
|
|
return jfjoch_git_date().empty() ? jfjoch_git_sha1() : jfjoch_git_sha1() + " " + jfjoch_git_date();
|
|
}
|
|
|
|
// The MTZ HISTORY records saying where a file came from, after the one naming the program: the
|
|
// commit, the command line and the input file. MTZ holds at most 30 records of 80 characters
|
|
// (mtzformat.html), so a long value runs on over as many records as it needs, up to that limit, each
|
|
// continuation starting with '+'. Readers strip the blanks at both ends of a record, so a record is
|
|
// never cut just after a space: the space starts the next one, behind its '+'.
|
|
void AddProvenanceHistory(std::vector<std::string> &history, const OutputProvenance &provenance) {
|
|
const auto add = [&history](const std::string &label, const std::string &value) {
|
|
if (value.empty())
|
|
return;
|
|
std::string rest = label + value;
|
|
while (!rest.empty() && history.size() < 30) {
|
|
size_t n = std::min<size_t>(80, rest.size());
|
|
while (n > 1 && n < rest.size() && rest[n - 1] == ' ')
|
|
--n;
|
|
history.push_back(rest.substr(0, n));
|
|
rest = n < rest.size() ? "+" + rest.substr(n) : std::string();
|
|
}
|
|
};
|
|
add("Rugnux git ", GitDescription());
|
|
add("Command: ", provenance.command_line);
|
|
add("Input: ", provenance.input_file);
|
|
add("Input file: ", provenance.input_identity);
|
|
add("Collected: ", provenance.collection_date);
|
|
}
|
|
|
|
// An e.s.d. as written: "." for a parameter the symmetry fixes (its e.s.d. is zero), as SHELX's ZERR
|
|
// writes 0 for it.
|
|
std::string EsdStr(double esd, int decimals) {
|
|
return esd < CELL_ESD_FIXED ? std::string(".") : Fmt(esd, decimals);
|
|
}
|
|
|
|
// _entry, _cell (with the e.s.d.s where there are any, and the _cell_measurement they came from) and
|
|
// _symmetry.
|
|
void WriteBlockHeader(std::ostream &out, const std::string &entry_id, const UnitCell &unitCell,
|
|
const gemmi::SpaceGroup &sg, const CellEsd *esd = nullptr) {
|
|
out << "_entry.id " << entry_id << "\n";
|
|
out << "#\n";
|
|
out << "_cell.entry_id " << entry_id << "\n";
|
|
out << "_cell.length_a " << Fmt(unitCell.a, 4) << "\n";
|
|
out << "_cell.length_b " << Fmt(unitCell.b, 4) << "\n";
|
|
out << "_cell.length_c " << Fmt(unitCell.c, 4) << "\n";
|
|
out << "_cell.angle_alpha " << Fmt(unitCell.alpha, 3) << "\n";
|
|
out << "_cell.angle_beta " << Fmt(unitCell.beta, 3) << "\n";
|
|
out << "_cell.angle_gamma " << Fmt(unitCell.gamma, 3) << "\n";
|
|
if (esd) {
|
|
out << "_cell.length_a_esd " << EsdStr(esd->esd[0], 5) << "\n";
|
|
out << "_cell.length_b_esd " << EsdStr(esd->esd[1], 5) << "\n";
|
|
out << "_cell.length_c_esd " << EsdStr(esd->esd[2], 5) << "\n";
|
|
out << "_cell.angle_alpha_esd " << EsdStr(esd->esd[3], 4) << "\n";
|
|
out << "_cell.angle_beta_esd " << EsdStr(esd->esd[4], 4) << "\n";
|
|
out << "_cell.angle_gamma_esd " << EsdStr(esd->esd[5], 4) << "\n";
|
|
out << "#\n";
|
|
out << "_cell_measurement.entry_id " << entry_id << "\n";
|
|
out << "_cell_measurement.reflns_used " << esd->reflections << "\n";
|
|
out << "_cell_measurement.theta_min " << Fmt(esd->theta_min_deg, 2) << "\n";
|
|
out << "_cell_measurement.theta_max " << Fmt(esd->theta_max_deg, 2) << "\n";
|
|
}
|
|
out << "#\n";
|
|
out << "_symmetry.entry_id " << entry_id << "\n";
|
|
out << "_symmetry.space_group_name_H-M " << CifStr(sg.hm) << "\n";
|
|
out << "_symmetry.Int_Tables_number " << sg.number << "\n";
|
|
out << "#\n";
|
|
}
|
|
|
|
// One output row per reflection in the standard CCP4 anomalous layout: the merged mean (IMEAN / F)
|
|
// plus the two Bijvoet mates (I(+)/I(-), F(+)/F(-)). The merge keeps the two mates as separate rows
|
|
// (I+ under the ASU representative hkl, I- under -hkl); this collapses them into one row so the MTZ
|
|
// and SHELX writers share one row list. has_anom is set false when NO reflection carries an anomalous
|
|
// split (e.g. the stills path) - callers then omit the +/- columns.
|
|
struct MergedOutRow {
|
|
int h = 0, k = 0, l = 0;
|
|
float Imean = NAN, sImean = NAN, Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN;
|
|
float Fmean = NAN, sFmean = NAN, Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN;
|
|
// The value written into FreeR_flag: 0 for the TEST set, 1 for the working set. That is the CCP4
|
|
// convention the column NAME belongs to - CCP4's own freerflag writes 0..19 with 0 as the ~5%
|
|
// test bin, and REFMAC's default FREE 0 reads it that way, so the opposite (phenix/CNS) numbering
|
|
// under this name stops REFMAC dead: "more than half of reflections are in free R set", then
|
|
// "Cannot switch free R flag", exit 1. phenix auto-detects either numbering with equal confidence
|
|
// (measured on both), and rugnux's own reference-MTZ reader already takes 0 as the test set, so
|
|
// this is the numbering that works everywhere and the one the rest of the code assumes.
|
|
int rfree = 1;
|
|
};
|
|
|
|
std::vector<MergedOutRow> BuildMergedRows(const std::vector<MergedReflection> &reflections,
|
|
const DiffractionExperiment &experiment,
|
|
bool &has_anom) {
|
|
std::vector<MergedOutRow> out_rows;
|
|
has_anom = true;
|
|
|
|
if (experiment.GetScalingSettings().GetMergeFriedel()) {
|
|
// Friedel-merged: IMEAN is the already-merged intensity (r.I). I(+)/I(-) are carried verbatim
|
|
// from the Bijvoet split the merge kept (rotation always does; scaled non-anomalously), so a weak
|
|
// anomalous signal is preserved without reprocessing. A reflection with only one mate, or a
|
|
// centric, gets a missing value (NaN) for the absent hand. When NO reflection has an anomalous
|
|
// split (e.g. the stills path, which does not compute one) the anomalous columns are omitted.
|
|
has_anom = std::any_of(reflections.begin(), reflections.end(),
|
|
[](const MergedReflection& r){ return std::isfinite(r.I_plus) || std::isfinite(r.I_minus); });
|
|
out_rows.reserve(reflections.size());
|
|
for (const auto& r : reflections)
|
|
out_rows.push_back({r.h, r.k, r.l, r.I, r.sigma, r.I_plus, r.sigma_plus, r.I_minus,
|
|
r.sigma_minus, r.F, r.sigmaF, r.F_plus, r.sigmaF_plus, r.F_minus,
|
|
r.sigmaF_minus, r.rfree_flag ? 0 : 1});
|
|
} else {
|
|
// Anomalous: group the two mates by their (shared) Friedel-merged ASU representative, then form
|
|
// IMEAN / F as their inverse-variance Friedel mean. A single generator gives both the group key
|
|
// (its hkl, identical for +hkl and -hkl) and which mate this row is (.plus).
|
|
const HKLKeyGenerator key_gen(false, experiment.GetSpaceGroupOrP1());
|
|
struct AnomRow {
|
|
int h = 0, k = 0, l = 0;
|
|
float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN;
|
|
float Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN;
|
|
int rfree = 1;
|
|
};
|
|
std::map<std::tuple<int, int, int>, AnomRow> rows;
|
|
for (const auto& r : reflections) {
|
|
const HKLKey key = key_gen(r);
|
|
AnomRow& row = rows[{key.h, key.k, key.l}];
|
|
row.h = key.h; row.k = key.k; row.l = key.l;
|
|
row.rfree = r.rfree_flag ? 0 : 1;
|
|
if (key.plus) { row.Ip = r.I; row.sIp = r.sigma; row.Fp = r.F; row.sFp = r.sigmaF; }
|
|
else { row.Im = r.I; row.sIm = r.sigma; row.Fm = r.F; row.sFm = r.sigmaF; }
|
|
}
|
|
|
|
// Friedel-mean of the two mates by inverse variance (the single mate, if only one was measured).
|
|
const auto combine = [](float a, float sa, float b, float sb, float& val, float& sig) {
|
|
const bool ok_a = std::isfinite(a) && sa > 0.0f;
|
|
const bool ok_b = std::isfinite(b) && sb > 0.0f;
|
|
if (ok_a && ok_b) {
|
|
const double wa = 1.0 / (static_cast<double>(sa) * sa);
|
|
const double wb = 1.0 / (static_cast<double>(sb) * sb);
|
|
val = static_cast<float>((wa * a + wb * b) / (wa + wb));
|
|
sig = static_cast<float>(1.0 / std::sqrt(wa + wb));
|
|
} else if (ok_a) { val = a; sig = sa; }
|
|
else if (ok_b) { val = b; sig = sb; }
|
|
else { val = NAN; sig = NAN; }
|
|
};
|
|
|
|
out_rows.reserve(rows.size());
|
|
for (const auto& [hkl, row] : rows) {
|
|
float i_mean, sig_i_mean, f_mean, sig_f_mean;
|
|
combine(row.Ip, row.sIp, row.Im, row.sIm, i_mean, sig_i_mean);
|
|
combine(row.Fp, row.sFp, row.Fm, row.sFm, f_mean, sig_f_mean);
|
|
out_rows.push_back({row.h, row.k, row.l, i_mean, sig_i_mean, row.Ip, row.sIp, row.Im,
|
|
row.sIm, f_mean, sig_f_mean, row.Fp, row.sFp, row.Fm, row.sFm, row.rfree});
|
|
}
|
|
}
|
|
return out_rows;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
|
|
void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
|
|
const UnitCell &unitCell,
|
|
const DiffractionExperiment &experiment,
|
|
const MergeStatistics &statistics,
|
|
const ErrorModelReport &error_model,
|
|
const TwinningAnalysisResult &twinning,
|
|
const OutputProvenance &provenance,
|
|
const std::string &filename,
|
|
size_t nthreads,
|
|
const CellEsd *cell_esd) {
|
|
|
|
std::ofstream out(filename);
|
|
if (!out)
|
|
throw std::runtime_error("WriteMmcifReflections: cannot open " + filename);
|
|
|
|
out << std::fixed;
|
|
|
|
// The keys PDBx asks for are written with the one value this file has for each: one entry, one
|
|
// crystal, one diffraction experiment (diffrn_id 1, which _pdbx_reflns_twin below also points at)
|
|
// and one wavelength.
|
|
const std::string entry_id = "sample";
|
|
const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1();
|
|
|
|
// ---------- data block ----------
|
|
out << "data_sample" << "\n";
|
|
out << "#\n";
|
|
WriteBlockHeader(out, entry_id, unitCell, sg, cell_esd);
|
|
|
|
// ---------- _audit ----------
|
|
out << "_audit.revision_id 1\n";
|
|
out << "_audit.creation_date " << CurrentDateISO() << "\n";
|
|
out << "_audit.update_record 'Initial release'\n";
|
|
out << "#\n";
|
|
|
|
// ---------- _software ----------
|
|
// Rugnux integrated the images (data reduction) and scaled and merged them (data scaling), the
|
|
// two rows the wwPDB asks a data-processing program to be named under; the description is the
|
|
// commit the build was made from.
|
|
const std::string software_description = CifStr("git " + GitDescription());
|
|
out << "loop_\n";
|
|
out << "_software.pdbx_ordinal\n";
|
|
out << "_software.name\n";
|
|
out << "_software.version\n";
|
|
out << "_software.classification\n";
|
|
out << "_software.description\n";
|
|
out << "1 Rugnux " << CifStr(jfjoch_version()) << " 'data reduction' " << software_description << "\n";
|
|
out << "2 Rugnux " << CifStr(jfjoch_version()) << " 'data scaling' " << software_description << "\n";
|
|
out << "#\n";
|
|
|
|
// ---------- provenance ----------
|
|
// How this file was made and from what: the command line and the input file, which no PDBx item
|
|
// holds, under the jfjoch local prefix as the other private items here. The two that can carry
|
|
// any character are text fields.
|
|
if (!provenance.command_line.empty())
|
|
out << "_jfjoch_provenance.command_line\n;" << provenance.command_line << "\n;\n";
|
|
if (!provenance.input_file.empty())
|
|
out << "_jfjoch_provenance.input_file\n;" << provenance.input_file << "\n;\n";
|
|
if (!provenance.input_identity.empty())
|
|
out << "_jfjoch_provenance.input_identity " << CifStr(provenance.input_identity) << "\n";
|
|
out << "#\n";
|
|
|
|
// ---------- _exptl_crystal / _diffrn ----------
|
|
out << "_exptl_crystal.id 1\n";
|
|
out << "#\n";
|
|
out << "_diffrn.id 1\n";
|
|
out << "_diffrn.crystal_id 1\n";
|
|
if (const auto t = experiment.GetSampleTemperature_K())
|
|
out << "_diffrn.ambient_temp " << Fmt(*t, 1) << "\n";
|
|
out << "#\n";
|
|
|
|
// ---------- _diffrn_radiation / _diffrn_source / _diffrn_detector ----------
|
|
out << "_diffrn_radiation.diffrn_id 1\n";
|
|
out << "_diffrn_radiation.wavelength_id 1\n";
|
|
out << "_diffrn_radiation.pdbx_diffrn_protocol 'SINGLE WAVELENGTH'\n";
|
|
out << "#\n";
|
|
out << "_diffrn_radiation_wavelength.id 1\n";
|
|
out << "_diffrn_radiation_wavelength.wavelength " << Fmt(experiment.GetWavelength_A(), 5) << "\n";
|
|
out << "#\n";
|
|
|
|
// The source and the detector as the input file names them: the facility and the beamline (NXmx
|
|
// /entry/source/name and /entry/instrument/name) and the detector's make and model. The general
|
|
// classes PDBx enumerates are written only where those names settle them.
|
|
out << "_diffrn_source.diffrn_id 1\n";
|
|
if (const std::string source = SourceClass(experiment.GetSourceType()); !source.empty())
|
|
out << "_diffrn_source.source " << CifStr(source) << "\n";
|
|
if (!experiment.GetSourceName().empty())
|
|
out << "_diffrn_source.pdbx_synchrotron_site " << CifStr(experiment.GetSourceName()) << "\n";
|
|
if (!experiment.GetInstrumentName().empty())
|
|
out << "_diffrn_source.pdbx_synchrotron_beamline " << CifStr(experiment.GetInstrumentName()) << "\n";
|
|
out << "_diffrn_source.pdbx_wavelength_list " << Fmt(experiment.GetWavelength_A(), 5) << "\n";
|
|
out << "#\n";
|
|
|
|
out << "_diffrn_detector.diffrn_id 1\n";
|
|
if (const std::string cls = DetectorClass(experiment.GetDetectorDescription()); !cls.empty())
|
|
out << "_diffrn_detector.detector " << CifStr(cls) << "\n";
|
|
out << "_diffrn_detector.type " << CifStr(experiment.GetDetectorDescription()) << "\n";
|
|
if (const std::string date = IsoDate(provenance.collection_date); !date.empty())
|
|
out << "_diffrn_detector.pdbx_collection_date " << date << "\n";
|
|
// Detector geometry actually used for integration (refined, when geometry refinement ran - the rotation
|
|
// two-pass or the stills global refinement update it on experiment_ before the written pass). jfjoch_
|
|
// local-data-name items: feedback of the distance / beam centre the data was reduced with.
|
|
out << "_diffrn_detector.jfjoch_distance_mm " << Fmt(experiment.GetDetectorDistance_mm(), 4) << "\n";
|
|
out << "_diffrn_detector.jfjoch_beam_center_x_pxl " << Fmt(experiment.GetBeamX_pxl(), 2) << "\n";
|
|
out << "_diffrn_detector.jfjoch_beam_center_y_pxl " << Fmt(experiment.GetBeamY_pxl(), 2) << "\n";
|
|
out << "#\n";
|
|
|
|
// ---------- merging statistics (_reflns overall + _reflns_shell loop) ----------
|
|
// cc_half and r_meas are stored as fractions (0-1), which is the mmCIF convention. ISa (the
|
|
// Diederichs asymptotic I/sigma, 1/b of the a*sigma^2 + (b*I)^2 error model) and the twinning
|
|
// indicators below have no standard mmCIF item. They are written under the "jfjoch" reserved
|
|
// prefix (_reflns.jfjoch_*), the IUCr-sanctioned local-data-name extension for private items -
|
|
// NOT the "pdbx_" prefix, which is owned by the wwPDB PDBx/mmCIF dictionary and must not label
|
|
// items that dictionary does not define. (The other pdbx_ items here are genuine PDBx items.)
|
|
const auto mult = [](const MergeStatisticsShell &s) {
|
|
return s.unique_reflections > 0 ? static_cast<double>(s.total_observations) / s.unique_reflections : 0.0; };
|
|
const auto compl_pct = [](const MergeStatisticsShell &s) {
|
|
return s.possible_unique_reflections > 0
|
|
? 100.0 * static_cast<double>(s.unique_reflections) / s.possible_unique_reflections : 0.0; };
|
|
if (!statistics.shells.empty()) {
|
|
const auto &ov = statistics.overall;
|
|
// Anomalous signal-to-noise (SigAno) is written only when an anomalous split was made, so a
|
|
// non-anomalous merge keeps its previous stats block / shell-loop columns unchanged.
|
|
const bool has_anom = std::isfinite(ov.abs_diff_over_sigma_anomalous);
|
|
out << "_reflns.entry_id " << entry_id << "\n";
|
|
out << "_reflns.pdbx_ordinal 1\n";
|
|
out << "_reflns.pdbx_diffrn_id 1\n";
|
|
out << "_reflns.d_resolution_high " << Fmt(ov.d_min, 2) << "\n";
|
|
out << "_reflns.d_resolution_low " << Fmt(ov.d_max, 2) << "\n";
|
|
out << "_reflns.number_obs " << ov.unique_reflections << "\n";
|
|
out << "_reflns.pdbx_number_measured_all " << ov.total_observations << "\n";
|
|
out << "_reflns.pdbx_redundancy " << Fmt(mult(ov), 2) << "\n";
|
|
out << "_reflns.percent_possible_obs " << Fmt(compl_pct(ov), 1) << "\n";
|
|
out << "_reflns.pdbx_netI_over_sigmaI " << Fmt(ov.mean_i_over_sigma, 2) << "\n";
|
|
out << "_reflns.pdbx_Rmerge_I_obs " << Fmt(ov.r_merge, 4) << "\n";
|
|
out << "_reflns.pdbx_Rrim_I_all " << Fmt(ov.r_meas, 4) << "\n";
|
|
out << "_reflns.pdbx_Rpim_I_all " << Fmt(ov.r_pim, 4) << "\n";
|
|
out << "_reflns.pdbx_CC_half " << Fmt(ov.cc_half, 4) << "\n";
|
|
out << "_reflns.pdbx_CC_star " << Fmt(CCStar(ov.cc_half), 4) << "\n";
|
|
if (has_anom)
|
|
out << "_reflns.pdbx_absDiff_over_sigma_anomalous " << Fmt(ov.abs_diff_over_sigma_anomalous, 3)
|
|
<< " # SigAno = <|dano|>/<sigma(dano)>\n";
|
|
// ISa in XDS's sense: the whole-range 1/sqrt(a*b) of the error model below, so this item can
|
|
// be read straight against a CORRECT.LP. The strong-reflection asymptote - a tier XDS does not
|
|
// have, and always the more optimistic of the two - is written separately rather than here.
|
|
if (!statistics.stills)
|
|
out << "_reflns.jfjoch_diffrn_ISa " << CifStr(error_model.isa)
|
|
<< " # 1/sqrt(a*b), the XDS convention\n";
|
|
if (!error_model.isa_asymptotic.empty())
|
|
out << "_reflns.jfjoch_diffrn_ISa_asymptotic " << CifStr(error_model.isa_asymptotic)
|
|
<< " # strong-reflection asymptote (Diederichs); rotation path only\n";
|
|
if (!error_model.a.empty())
|
|
out << "_reflns.jfjoch_error_model_a " << CifStr(error_model.a)
|
|
<< " # sigma^2 = a*(sigma0^2 + b*I^2), XDS convention\n";
|
|
if (!error_model.b.empty())
|
|
out << "_reflns.jfjoch_error_model_b " << CifStr(error_model.b) << "\n";
|
|
if (!error_model.c.empty())
|
|
out << "_reflns.jfjoch_error_model_c " << CifStr(error_model.c)
|
|
<< " # stills partiality term: + c*((1-p)/p)^2*I^2\n";
|
|
// Dataset-wide isotropic Wilson B-factor estimate (standard PDBx item), analogous to XDS's
|
|
// "WILSON LINE ... B=". Emitted only when the log-linear fit succeeded.
|
|
if (std::isfinite(statistics.wilson_b) && statistics.wilson_b > 0.0)
|
|
out << "_reflns.B_iso_Wilson_estimate " << Fmt(statistics.wilson_b, 2) << "\n";
|
|
// Twinning indicators (no standard mmCIF item; same jfjoch local prefix as ISa above).
|
|
if (twinning.l_test_pairs > 0) {
|
|
out << "_reflns.jfjoch_L_test_mean_abs_L " << Fmt(twinning.mean_abs_l, 3)
|
|
<< " # Padilla-Yeates <|L|> (untwinned 0.500, perfect twin 0.375)\n";
|
|
out << "_reflns.jfjoch_L_test_mean_L_squared " << Fmt(twinning.mean_l_squared, 3)
|
|
<< " # <L^2> (untwinned 0.333, perfect twin 0.200)\n";
|
|
}
|
|
if (twinning.moment_reflections > 0)
|
|
out << "_reflns.jfjoch_second_moment_I " << Fmt(twinning.second_moment, 3)
|
|
<< " # <I^2>/<I>^2 (untwinned 2.00, perfect twin 1.50)\n";
|
|
// A twin law the run identified (TwinningAnalysisResult::twin_law_operator): the standard
|
|
// PDBx loop, one row per domain, the second domain related by the operator. Only meaningful
|
|
// in the true (sub)group's file, which is the one this is - in the higher group the operator
|
|
// is a symmetry element. Merohedral: the law is a lattice operation of the crystal system.
|
|
if (!twinning.twin_law_operator.empty()) {
|
|
out << "#\n";
|
|
out << "loop_\n";
|
|
out << "_pdbx_reflns_twin.domain_id\n";
|
|
out << "_pdbx_reflns_twin.crystal_id\n";
|
|
out << "_pdbx_reflns_twin.diffrn_id\n";
|
|
out << "_pdbx_reflns_twin.operator\n";
|
|
out << "_pdbx_reflns_twin.fraction\n";
|
|
out << "_pdbx_reflns_twin.type\n";
|
|
out << "1 1 1 'h,k,l' " << Fmt(1.0 - twinning.twin_law_fraction, 3) << " merohedral\n";
|
|
out << "2 1 1 " << CifStr(twinning.twin_law_operator) << " " << Fmt(twinning.twin_law_fraction, 3)
|
|
<< " merohedral\n";
|
|
}
|
|
// Radiation-damage monitor (rotation): the relative Debye-Waller B change from the first to the last
|
|
// frame (A^2). A large magnitude flags a dose-dependent resolution-scale change = radiation damage;
|
|
// positive is the typical direction (high-resolution intensity fades with dose). No standard mmCIF item.
|
|
if (std::isfinite(statistics.radiation_damage_delta_b))
|
|
out << "_reflns.jfjoch_radiation_damage_relative_B " << Fmt(statistics.radiation_damage_delta_b, 2)
|
|
<< " # relative-B first->last over the run (A^2); + = high-res fades with dose\n";
|
|
// Diffraction anisotropy. The eigen-decomposition of the anisotropy tensor has standard PDBx
|
|
// items; the eigenvalues there must be non-negative, and only the deviatoric part of the tensor
|
|
// is determined at all (its isotropic part is degenerate with the overall scale), so they are
|
|
// written relative to the weakest direction - eigenvalue_3 is 0 by construction and
|
|
// eigenvalue_1 is the anisotropic deltaB. The eigenvectors are in the PDB orthogonalisation
|
|
// convention, which is the one gemmi (and hence rugnux) uses throughout.
|
|
//
|
|
// The directional diffraction LIMITS are deliberately NOT written as
|
|
// _reflns.pdbx_aniso_diffraction_limit_*: the dictionary defines those as the semi-axes of an
|
|
// ellipsoid fitted to a diffraction cut-off surface, which is a different construction from the
|
|
// one below and one rugnux does not perform - it cuts nothing on a directional criterion. They
|
|
// go under the jfjoch local prefix with their own definition instead.
|
|
const auto &an = statistics.anisotropy;
|
|
if (an.n_reflections > 0 && an.n_cells > 0 && std::isfinite(an.delta_b)) {
|
|
out << "_reflns.pdbx_orthogonalization_convention pdb\n";
|
|
for (int i = 0; i < 3; ++i) {
|
|
out << "_reflns.pdbx_aniso_B_tensor_eigenvalue_" << (i + 1) << " "
|
|
<< Fmt(an.eigenvalue[i] - an.eigenvalue[2], 2) << "\n";
|
|
for (int j = 0; j < 3; ++j)
|
|
out << "_reflns.pdbx_aniso_B_tensor_eigenvector_" << (i + 1) << "_ortho[" << (j + 1)
|
|
<< "] " << Fmt(an.eigenvector[i][j], 4) << "\n";
|
|
}
|
|
out << "_reflns.jfjoch_aniso_delta_B " << Fmt(an.delta_b, 2)
|
|
<< " # range of the principal components (A^2), fitted on intensities\n";
|
|
if (std::isfinite(an.delta_b_linear))
|
|
out << "_reflns.jfjoch_aniso_delta_B_linear " << Fmt(an.delta_b_linear, 2)
|
|
<< " # deltaB implied by the s^2 slope alone: what a Debye-Waller B accounts for\n";
|
|
out << "_reflns.jfjoch_aniso_shape " << AnisotropyShapeCode(an.shape)
|
|
<< " # resolution signature of the directional deficit\n";
|
|
if (std::isfinite(an.floor))
|
|
out << "_reflns.jfjoch_aniso_floor " << Fmt(an.floor, 3)
|
|
<< " # deltaB this data set's own systematic error could manufacture (A^2)\n";
|
|
if (std::isfinite(an.significance))
|
|
out << "_reflns.jfjoch_aniso_significance " << Fmt(an.significance, 2)
|
|
<< " # deltaB(linear) / floor\n";
|
|
out << "_reflns.jfjoch_aniso_verdict " << AnisotropyVerdictCode(an.verdict)
|
|
<< " # DETECTED / NOT_DETECTED / CANNOT_DETERMINE: deltaB(linear)/floor > 3.5, or the directional limits\n";
|
|
for (int i = 0; i < 3; ++i)
|
|
if (std::isfinite(an.d_min_axis[i]))
|
|
out << "_reflns.jfjoch_aniso_d_min_" << (i + 1) << " "
|
|
<< Fmt(an.d_min_axis[i], 2)
|
|
<< " # <I/sigma(I)> = 2 in a 20 deg cone about eigenvector " << (i + 1)
|
|
<< (an.d_min_censored[i] ? "; at the edge of the measured data, so a bound" : "")
|
|
<< "\n";
|
|
}
|
|
out << "#\n";
|
|
|
|
// Per-batch relative-B curve (the radiation-damage monitor, rotation): one relative Debye-Waller B
|
|
// per rotation-range batch, measured before any correction. rotation_start_deg = id * batch_deg.
|
|
if (!statistics.radiation_damage_b_batch.empty()) {
|
|
out << "loop_\n";
|
|
out << "_jfjoch_radiation_damage_batch.id\n";
|
|
out << "_jfjoch_radiation_damage_batch.rotation_start_deg\n";
|
|
out << "_jfjoch_radiation_damage_batch.relative_B\n";
|
|
for (size_t i = 0; i < statistics.radiation_damage_b_batch.size(); ++i)
|
|
out << " " << (i + 1) << " "
|
|
<< Fmt(static_cast<double>(i) * statistics.radiation_damage_batch_deg, 1) << " "
|
|
<< Fmt(statistics.radiation_damage_b_batch[i], 2) << "\n";
|
|
out << "#\n";
|
|
}
|
|
|
|
out << "loop_\n";
|
|
out << "_reflns_shell.pdbx_ordinal\n";
|
|
out << "_reflns_shell.pdbx_diffrn_id\n";
|
|
out << "_reflns_shell.d_res_high\n";
|
|
out << "_reflns_shell.d_res_low\n";
|
|
out << "_reflns_shell.number_measured_obs\n";
|
|
out << "_reflns_shell.number_unique_obs\n";
|
|
out << "_reflns_shell.pdbx_redundancy\n";
|
|
out << "_reflns_shell.percent_possible_obs\n";
|
|
out << "_reflns_shell.meanI_over_sigI_obs\n";
|
|
out << "_reflns_shell.Rmerge_I_obs\n";
|
|
out << "_reflns_shell.pdbx_Rrim_I_all\n";
|
|
out << "_reflns_shell.pdbx_Rpim_I_all\n";
|
|
out << "_reflns_shell.pdbx_CC_half\n";
|
|
out << "_reflns_shell.pdbx_CC_star\n";
|
|
if (has_anom)
|
|
out << "_reflns_shell.pdbx_absDiff_over_sigma_anomalous\n";
|
|
int ordinal = 0;
|
|
for (const auto &s : statistics.shells) {
|
|
if (s.unique_reflections == 0)
|
|
continue;
|
|
out << ++ordinal << " 1 "
|
|
<< Fmt(s.d_min, 2) << " " << Fmt(s.d_max, 2) << " "
|
|
<< s.total_observations << " " << s.unique_reflections << " "
|
|
<< Fmt(mult(s), 2) << " " << Fmt(compl_pct(s), 1) << " "
|
|
<< Fmt(s.mean_i_over_sigma, 2) << " " << Fmt(s.r_merge, 4) << " " << Fmt(s.r_meas, 4) << " "
|
|
<< Fmt(s.r_pim, 4) << " " << Fmt(s.cc_half, 4) << " " << Fmt(CCStar(s.cc_half), 4);
|
|
if (has_anom)
|
|
out << " " << Fmt(s.abs_diff_over_sigma_anomalous, 3);
|
|
out << "\n";
|
|
}
|
|
out << "#\n";
|
|
}
|
|
|
|
// ---------- _refln loop ----------
|
|
out << "loop_\n";
|
|
out << "_refln.index_h\n";
|
|
out << "_refln.index_k\n";
|
|
out << "_refln.index_l\n";
|
|
out << "_refln.intensity_meas\n";
|
|
out << "_refln.intensity_sigma\n";
|
|
out << "_refln.pdbx_I_plus\n";
|
|
out << "_refln.pdbx_I_plus_sigma\n";
|
|
out << "_refln.pdbx_I_minus\n";
|
|
out << "_refln.pdbx_I_minus_sigma\n";
|
|
out << "_refln.F_meas_au\n";
|
|
out << "_refln.F_meas_sigma_au\n";
|
|
out << "_refln.pdbx_F_plus\n";
|
|
out << "_refln.pdbx_F_plus_sigma\n";
|
|
out << "_refln.pdbx_F_minus\n";
|
|
out << "_refln.pdbx_F_minus_sigma\n";
|
|
out << "_refln.status\n";
|
|
|
|
// One row per unique reflection, twelve formatted floats each - tens of megabytes on a crowded
|
|
// crystal, and the largest single-threaded stretch left in a run. Nothing about a row depends on
|
|
// any other, so each worker formats its own block into its own string and the blocks go to the
|
|
// file in order. The columns are written exactly as the stream wrote them: the same "%.4f" (or
|
|
// "?" where the value is not finite), right-aligned in the same width.
|
|
{
|
|
const size_t nrow = reflections.size();
|
|
const size_t nw = std::max<size_t>(nthreads, 1);
|
|
const int nch = static_cast<int>(ThreadsForWork(nrow, nw, 4096));
|
|
std::vector<std::string> block(nch);
|
|
ParallelChunks(nch, nw, [&](int tlo, int thi) {
|
|
for (int t = tlo; t < thi; ++t) {
|
|
const size_t lo = nrow * t / nch, hi = nrow * (t + 1) / nch;
|
|
std::string &s = block[t];
|
|
s.reserve((hi - lo) * 208);
|
|
const auto column = [&s](const std::string &v, size_t width) {
|
|
if (v.size() < width) s.append(width - v.size(), ' ');
|
|
s.append(v);
|
|
s.push_back(' ');
|
|
};
|
|
for (size_t i = lo; i < hi; ++i) {
|
|
const auto &r = reflections[i];
|
|
column(std::to_string(r.h), 5);
|
|
column(std::to_string(r.k), 5);
|
|
column(std::to_string(r.l), 5);
|
|
column(Fmt(r.I, 4), 14);
|
|
column(Fmt(r.sigma, 4), 14);
|
|
column(Fmt(r.I_plus, 4), 14);
|
|
column(Fmt(r.sigma_plus, 4), 14);
|
|
column(Fmt(r.I_minus, 4), 14);
|
|
column(Fmt(r.sigma_minus, 4), 14);
|
|
column(Fmt(r.F, 4), 14);
|
|
column(Fmt(r.sigmaF, 4), 14);
|
|
column(Fmt(r.F_plus, 4), 14);
|
|
column(Fmt(r.sigmaF_plus, 4), 14);
|
|
column(Fmt(r.F_minus, 4), 14);
|
|
column(Fmt(r.sigmaF_minus, 4), 14);
|
|
// PDBx/mmCIF _refln.status: 'f' marks the free (test) set, 'o' the working
|
|
// set. This is the only place the free set is carried in the file - cif2mtz and
|
|
// gemmi both read it, and neither knows any other reflection-level free-set item.
|
|
s.append(r.rfree_flag ? "f\n" : "o\n");
|
|
}
|
|
}
|
|
});
|
|
for (const std::string &s : block)
|
|
out.write(s.data(), static_cast<std::streamsize>(s.size()));
|
|
}
|
|
|
|
out << "#\n";
|
|
out << "# End of reflections\n";
|
|
|
|
out.close();
|
|
}
|
|
|
|
MtzTable MergedMtzTable(const std::vector<MergedReflection> &reflections,
|
|
const UnitCell &unitCell,
|
|
const DiffractionExperiment &experiment,
|
|
const OutputProvenance &provenance) {
|
|
MtzTable m;
|
|
m.unit_cell = unitCell;
|
|
m.space_group_name = experiment.GetSpaceGroupOrP1().xhm();
|
|
// Producing-software provenance in the MTZ header (title + HISTORY, the CCP4 convention).
|
|
m.title = "Rugnux merged reflections";
|
|
m.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction");
|
|
AddProvenanceHistory(m.history, provenance);
|
|
// HKL_base first, so the data dataset is id 1: a dataset written as id 0 is the reserved base dataset
|
|
// as far as CCP4's mtzlib is concerned, and the wavelength on it is dropped - mtzinfo then reports the
|
|
// Cu K-alpha default 1.54187 A, which is what a program that takes f'/f'' from the file would phase on.
|
|
m.datasets = {{"HKL_base", "HKL_base", "HKL_base", 0.0f},
|
|
{"native", experiment.GetSampleName(), "native", experiment.GetWavelength_A()}};
|
|
|
|
// One row per reflection in the CCP4 anomalous layout (IMEAN + I(+)/I(-), and the same split for
|
|
// the French-Wilson amplitude), which aimless / ctruncate / mtz2sca / ANODE read directly.
|
|
bool has_anom = true;
|
|
const std::vector<MergedOutRow> out_rows = BuildMergedRows(reflections, experiment, has_anom);
|
|
std::vector<std::pair<const char *, char>> columns{{"H", 'H'}, {"K", 'H'}, {"L", 'H'},
|
|
{"IMEAN", 'J'}, {"SIGIMEAN", 'Q'}};
|
|
if (has_anom)
|
|
columns.insert(columns.end(), {{"I(+)", 'K'}, {"SIGI(+)", 'M'}, {"I(-)", 'K'}, {"SIGI(-)", 'M'}});
|
|
columns.insert(columns.end(), {{"F", 'F'}, {"SIGF", 'Q'}}); // F: the French-Wilson amplitude
|
|
if (has_anom)
|
|
columns.insert(columns.end(), {{"F(+)", 'G'}, {"SIGF(+)", 'L'}, {"F(-)", 'G'}, {"SIGF(-)", 'L'}});
|
|
columns.push_back({"FreeR_flag", 'I'});
|
|
for (const auto &[label, type] : columns)
|
|
m.columns.push_back({label, type, type == 'H' ? 0 : 1, {}});
|
|
for (const auto &row : out_rows) {
|
|
std::vector<float> v{static_cast<float>(row.h), static_cast<float>(row.k), static_cast<float>(row.l),
|
|
row.Imean, row.sImean};
|
|
if (has_anom)
|
|
v.insert(v.end(), {row.Ip, row.sIp, row.Im, row.sIm});
|
|
v.insert(v.end(), {row.Fmean, row.sFmean});
|
|
if (has_anom)
|
|
v.insert(v.end(), {row.Fp, row.sFp, row.Fm, row.sFm});
|
|
v.push_back(static_cast<float>(row.rfree));
|
|
for (size_t c = 0; c < v.size(); c++)
|
|
m.columns[c].values.push_back(v[c]);
|
|
}
|
|
return m;
|
|
}
|
|
|
|
void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
|
|
const UnitCell &unitCell,
|
|
const DiffractionExperiment &experiment,
|
|
const OutputProvenance &provenance,
|
|
const std::string &filename) {
|
|
WriteMtz(MergedMtzTable(reflections, unitCell, experiment, provenance), filename);
|
|
}
|
|
|
|
namespace {
|
|
// The SHELX HKLF 4 record, FORMAT(3I4,2F8.2): h k l I sigma(I). "%.2f" right-aligned in the fixed
|
|
// field is what `fixed` + `setprecision(2)` + `setw` made the stream write.
|
|
void AppendHklf4(std::string &s, int h, int k, int l, double I, double sigma) {
|
|
const auto column = [&s](const std::string &v, size_t width) {
|
|
if (v.size() < width) s.append(width - v.size(), ' ');
|
|
s.append(v);
|
|
};
|
|
const auto num2 = [](double v) { return FixedDigits(v, 2); };
|
|
column(std::to_string(h), 4);
|
|
column(std::to_string(k), 4);
|
|
column(std::to_string(l), 4);
|
|
column(num2(I), 8);
|
|
column(num2(sigma), 8);
|
|
s.push_back('\n');
|
|
}
|
|
|
|
// An HKLF 4 file of nrow source rows, each appending its records with append_row(i, s). Built in
|
|
// parallel into per-worker blocks and handed to the file in order, exactly as the mmCIF rows are,
|
|
// and closed with the 0 0 0 end-of-data record.
|
|
template <class AppendRow>
|
|
void WriteHklf4File(const std::string &filename, size_t nrow, size_t nthreads, AppendRow append_row) {
|
|
std::ofstream out(filename);
|
|
if (!out)
|
|
throw std::runtime_error("WriteShelxHklReflections: cannot open " + filename);
|
|
const size_t nw = std::max<size_t>(nthreads, 1);
|
|
const int nch = static_cast<int>(ThreadsForWork(nrow, nw, 4096));
|
|
std::vector<std::string> block(nch);
|
|
ParallelChunks(nch, nw, [&](int tlo, int thi) {
|
|
for (int t = tlo; t < thi; ++t) {
|
|
const size_t lo = nrow * t / nch, hi = nrow * (t + 1) / nch;
|
|
std::string &s = block[t];
|
|
s.reserve((hi - lo) * 2 * 29);
|
|
for (size_t i = lo; i < hi; ++i)
|
|
append_row(i, s);
|
|
}
|
|
});
|
|
for (const std::string &s : block)
|
|
out.write(s.data(), static_cast<std::streamsize>(s.size()));
|
|
std::string tail;
|
|
AppendHklf4(tail, 0, 0, 0, 0.0, 0.0);
|
|
out.write(tail.data(), static_cast<std::streamsize>(tail.size()));
|
|
}
|
|
|
|
// Intensities are put on a common scale so the largest value fits the F8.2 field (the absolute scale
|
|
// is irrelevant to every SHELX program, which refines or ignores it); I and sigma share the scale.
|
|
double Hklf4Scale(double max_abs) {
|
|
return (std::isfinite(max_abs) && max_abs > 0.0) ? 9999.0 / max_abs : 1.0;
|
|
}
|
|
|
|
bool UsableForHkl(float v, float s) { return std::isfinite(v) && std::isfinite(s) && s > 0.0f; }
|
|
}
|
|
|
|
void WriteShelxHklReflections(const std::vector<MergedReflection> &reflections,
|
|
const DiffractionExperiment &experiment,
|
|
const std::string &filename,
|
|
size_t nthreads) {
|
|
bool has_anom = true;
|
|
const std::vector<MergedOutRow> rows = BuildMergedRows(reflections, experiment, has_anom);
|
|
|
|
// One record per merged reflection. The Bijvoet mates are written separately - I(+) at +hkl, I(-)
|
|
// at -hkl - so the anomalous differences survive; a reflection with no anomalous split is written
|
|
// once as its mean.
|
|
double max_abs = 0.0;
|
|
for (const auto& r : rows) {
|
|
if (UsableForHkl(r.Ip, r.sIp)) max_abs = std::max({max_abs, std::fabs(double(r.Ip)), double(r.sIp)});
|
|
if (UsableForHkl(r.Im, r.sIm)) max_abs = std::max({max_abs, std::fabs(double(r.Im)), double(r.sIm)});
|
|
if (!UsableForHkl(r.Ip, r.sIp) && !UsableForHkl(r.Im, r.sIm) && UsableForHkl(r.Imean, r.sImean))
|
|
max_abs = std::max({max_abs, std::fabs(double(r.Imean)), double(r.sImean)});
|
|
}
|
|
const double scale = Hklf4Scale(max_abs);
|
|
WriteHklf4File(filename, rows.size(), nthreads, [&](size_t i, std::string &s) {
|
|
const auto &r = rows[i];
|
|
const bool plus = UsableForHkl(r.Ip, r.sIp);
|
|
const bool minus = UsableForHkl(r.Im, r.sIm);
|
|
if (plus) AppendHklf4(s, r.h, r.k, r.l, scale * r.Ip, scale * r.sIp);
|
|
if (minus) AppendHklf4(s, -r.h, -r.k, -r.l, scale * r.Im, scale * r.sIm);
|
|
if (!plus && !minus && UsableForHkl(r.Imean, r.sImean))
|
|
AppendHklf4(s, r.h, r.k, r.l, scale * r.Imean, scale * r.sImean);
|
|
});
|
|
}
|
|
|
|
void WriteShelxHklReflections(const std::vector<ScaledFull> &fulls,
|
|
const std::string &filename,
|
|
size_t nthreads) {
|
|
// One record per full, at the index it was measured at: SHELXL / SHELXC average the equivalents
|
|
// themselves and report Rint and Rsigma from them, which a merged file hides. No batch number - in
|
|
// HKLF 4 that column selects a BASF scale factor in SHELXL, which these fulls do not want.
|
|
double max_abs = 0.0;
|
|
for (const auto &f : fulls)
|
|
if (UsableForHkl(f.I, f.sigma))
|
|
max_abs = std::max({max_abs, std::fabs(double(f.I)), double(f.sigma)});
|
|
const double scale = Hklf4Scale(max_abs);
|
|
WriteHklf4File(filename, fulls.size(), nthreads, [&](size_t i, std::string &s) {
|
|
const auto &f = fulls[i];
|
|
if (UsableForHkl(f.I, f.sigma))
|
|
AppendHklf4(s, f.h, f.k, f.l, scale * f.I, scale * f.sigma);
|
|
});
|
|
}
|
|
|
|
void WriteShelxIns(const UnitCell &cell, const CellEsd *esd, const DiffractionExperiment &experiment,
|
|
const std::string &title, const std::string &filename) {
|
|
const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1();
|
|
const gemmi::GroupOps ops = sg.operations();
|
|
// LATT is the centring, positive where the group holds an inversion centre at the origin, which
|
|
// SHELX then generates itself - so only one of each -1 pair is listed (the proper rotations). A
|
|
// centring SHELX has no code for gets LATT -1 with every operator spelled out.
|
|
const std::string centrings = "PIRFABC";
|
|
const char centring = ops.find_centering();
|
|
const size_t latt_code = centring ? centrings.find(centring) : std::string::npos;
|
|
bool origin_inversion = false;
|
|
if (const gemmi::Op *inv = ops.find_by_rotation(gemmi::Op::inversion_rot()))
|
|
for (const auto &t : ops.cen_ops)
|
|
origin_inversion |= inv->tran == t;
|
|
const bool centro = latt_code != std::string::npos && origin_inversion;
|
|
int latt = latt_code != std::string::npos ? static_cast<int>(latt_code) + 1 : 1;
|
|
std::vector<gemmi::Op> symm;
|
|
for (const auto &op : latt_code != std::string::npos ? ops.sym_ops : ops.all_ops_sorted())
|
|
if (!(op == gemmi::Op::identity()) && !(centro && op.det_rot() < 0))
|
|
symm.push_back(op);
|
|
|
|
// SFAC and UNIT have to be the formula, which the diffraction does not tell: an organic one is put
|
|
// in its place at 18 A^3 per non-hydrogen atom, for the user to replace.
|
|
const int nonh = std::max(1, static_cast<int>(std::lround(gemmi::UnitCell(cell).volume / 18.0)));
|
|
const int n_c = std::max(1, static_cast<int>(std::lround(0.7 * nonh)));
|
|
const int n_n = std::max(1, static_cast<int>(std::lround(0.15 * nonh)));
|
|
const int n_o = std::max(1, nonh - n_c - n_n);
|
|
|
|
std::ofstream out(filename);
|
|
if (!out)
|
|
throw std::runtime_error("WriteShelxIns: cannot open " + filename);
|
|
// SHELXL stops on an instruction longer than 80 characters, so every line is kept within that.
|
|
char line[256];
|
|
out << ("TITL " + title + " in " + sg.xhm()).substr(0, 80) << "\n";
|
|
out << "REM Rugnux " << jfjoch_version() << ": CELL, ZERR, LATT and SYMM from the data.\n";
|
|
out << "REM SFAC/UNIT are a placeholder, an organic formula at 18 A^3 per non-H atom:\n";
|
|
out << "REM replace them with the real formula before refining with SHELXL.\n";
|
|
if (!esd)
|
|
out << "REM The cell e.s.d.s were not measured in this run: ZERR gives them as 0.\n";
|
|
std::snprintf(line, sizeof(line), "CELL %.5f %.5f %.5f %.5f %.4f %.4f %.4f\n", experiment.GetWavelength_A(),
|
|
cell.a, cell.b, cell.c, cell.alpha, cell.beta, cell.gamma);
|
|
out << line;
|
|
// ZERR is required by SHELXL (it stops with "ZERR MISSING"). Its first number is Z, the formula units in
|
|
// the cell: one per general position is assumed. Where the run did not measure the e.s.d.s they are 0.
|
|
const std::array<double, 6> e = esd ? esd->esd : std::array<double, 6>{};
|
|
std::snprintf(line, sizeof(line), "ZERR %d %.5f %.5f %.5f %.4f %.4f %.4f\n", ops.order(), e[0], e[1], e[2],
|
|
e[3], e[4], e[5]);
|
|
out << line;
|
|
out << "LATT " << (centro ? latt : -latt) << "\n";
|
|
for (const auto &op : symm) {
|
|
std::string t = op.triplet();
|
|
for (auto &ch : t)
|
|
ch = static_cast<char>(std::toupper(static_cast<unsigned char>(ch)));
|
|
out << "SYMM " << t << "\n";
|
|
}
|
|
out << "SFAC C H N O\n";
|
|
out << "UNIT " << n_c << " " << n_c << " " << n_n << " " << n_o << "\n";
|
|
// TEMP, in Celsius, from the sample temperature the data file records (the mmCIF's _diffrn.ambient_temp,
|
|
// in K): SHELXL takes its starting U values and the CIF's temperature from it, and assumes 20 C without it.
|
|
if (const auto t = experiment.GetSampleTemperature_K()) {
|
|
std::snprintf(line, sizeof(line), "TEMP %.1f\n", *t - 273.15);
|
|
out << line;
|
|
} else
|
|
out << "REM No TEMP: the data file records no sample temperature (SHELXL assumes 20 C).\n";
|
|
out << "HKLF 4\n";
|
|
out << "END\n";
|
|
}
|
|
|
|
namespace {
|
|
|
|
// The detector position an observation is written at: a reflection too weak to have a measured
|
|
// centroid still has a predicted one.
|
|
float DetectorX(const Reflection &r) { return std::isfinite(r.observed_x) ? r.observed_x : r.predicted_x; }
|
|
float DetectorY(const Reflection &r) { return std::isfinite(r.observed_y) ? r.observed_y : r.predicted_y; }
|
|
|
|
// Sum each rocking event into one full observation. A rotation reflection is integrated image by
|
|
// image, so it arrives here as a run of partials over consecutive frames; the run is cut where the
|
|
// 3D combine cuts it - same raw hkl, frames no further apart than RockingEventFrameGap - so the exported
|
|
// file and rugnux's own merge see exactly the same events.
|
|
// The parts are added, plainly, with their variances in quadrature, which is what every other
|
|
// rotation program writes as a full. Nothing is divided by the partiality: FRACTIONCALC carries the
|
|
// summed rocking-curve fraction, and an event the sweep cut short says so there rather than being
|
|
// scaled up or dropped here.
|
|
// The metadata a full carries is the partiality-weighted mean of its parts - the centroid of the
|
|
// rocking curve, which is what XD/YD/ZD mean in an XDS file and what a scale model downstream reads
|
|
// these columns for. An intensity weighting is the other candidate and it does not survive weak
|
|
// data: the parts of a weak reflection scatter about zero, so on the crystals this was measured on
|
|
// 8-41% of events have parts summing to zero or less and 10-21% of the intensity-weighted centroids
|
|
// fall outside the event's own frame range. Partialities are the rocking curve's own weights and are
|
|
// positive by construction.
|
|
// The floor is the SAME pair the 3D combine applies to the same event - --min-partiality (default
|
|
// 0.02) and --min-captured-fraction (default 0.7 on rotation) - so what is exported as a full is
|
|
// exactly what rugnux merges as a full. The captured-fraction floor is what makes the FULL
|
|
// declaration true: the batch header says LDTYPE=2 and M/ISYM carries M=0, AIMLESS reads that as
|
|
// "all runs have only fulls" and ignores FRACTIONCALC entirely, so an event that caught a twentieth
|
|
// of its rocking curve does not enter the merge at a twentieth of its intensity - it enters whole,
|
|
// with a small sigma, and is weighted heavily. There is no cut for the reader to make: the column
|
|
// that would let it make one is one the reading program has been told to ignore. Measured, dropping
|
|
// those events moves AIMLESS's Rmerge at 1.8 A from 1.353 to 0.694 and CC(1/2) from 0.985 to 0.993.
|
|
//
|
|
// FRACTIONCALC itself is written as summed, above 1 included. Each part's partiality is the erf pair
|
|
// the predictor computed on THAT frame, from that frame's own refined lattice and its own mosaicity
|
|
// estimate, so the parts of one event do not tile the rocking curve exactly - measured, the offset
|
|
// from one part to the next steps by the wedge to within 12% of it, and the sums land in a peak at
|
|
// 1.000 whose 95th percentile is 1.09. It is the honest estimate of a captured fraction and it is
|
|
// not the number rugnux scales on (the merge recomputes every partiality from a frame-order-smoothed
|
|
// mosaicity, RotationScaleMerge::SmoothMosaicity), so clamping it to 1 would hide the spread and buy
|
|
// a reading program nothing.
|
|
//
|
|
// The fulls come back in pieces, one per worker, which laid end to end are all of them in (h,k,l)
|
|
// order - joining them would be one more serial copy of the whole list.
|
|
std::vector<std::vector<Reflection>> SumRockingEvents(const std::vector<IntegrationOutcome> &outcomes,
|
|
double min_partiality, double min_captured_fraction,
|
|
float wedge_deg, size_t nthreads) {
|
|
const float max_frame_gap = RockingEventFrameGap(wedge_deg); // == RotationScaleMerge's
|
|
|
|
// The sort key travels with the part instead of being read back from the reflection, the way the
|
|
// merge's own ingest sort carries it (RotationScaleMerge's SortKey): there are millions of parts
|
|
// and an indirect compare is a cache miss on every one of them. Two parts can genuinely share
|
|
// (h,k,l) and image_number, and the event sums below are floating point, so the part's position
|
|
// in the input - its outcome, then its place there - is the last key: the order is total, and the
|
|
// sort gives the same one on any number of workers.
|
|
struct Part {
|
|
int32_t h, k, l;
|
|
float image_number;
|
|
uint32_t o, n; // the outcome, and the part's place in it
|
|
};
|
|
std::vector<size_t> first(outcomes.size() + 1, 0);
|
|
for (size_t o = 0; o < outcomes.size(); ++o)
|
|
first[o + 1] = first[o] + outcomes[o].reflections.size();
|
|
// Left uninitialised, so that the pages are first touched by the workers filling them.
|
|
const size_t n_parts = first.back();
|
|
const auto parts = std::make_unique_for_overwrite<Part[]>(n_parts);
|
|
ParallelFor(static_cast<int>(outcomes.size()), nthreads, [&](int o) {
|
|
for (size_t m = 0; m < outcomes[o].reflections.size(); ++m) {
|
|
const Reflection &r = outcomes[o].reflections[m];
|
|
parts[first[o] + m] = {r.h, r.k, r.l, r.image_number, static_cast<uint32_t>(o), static_cast<uint32_t>(m)};
|
|
}
|
|
});
|
|
ParallelSort(parts.get(), parts.get() + n_parts, nthreads, [](const Part &a, const Part &b) {
|
|
return std::tie(a.h, a.k, a.l, a.image_number, a.o, a.n) < std::tie(b.h, b.k, b.l, b.image_number, b.o, b.n);
|
|
});
|
|
|
|
// Part i continues the event of part i - 1.
|
|
const auto continues = [&](size_t i) {
|
|
return parts[i].h == parts[i - 1].h && parts[i].k == parts[i - 1].k && parts[i].l == parts[i - 1].l
|
|
&& parts[i].image_number - parts[i - 1].image_number <= max_frame_gap;
|
|
};
|
|
|
|
// The events are summed in chunks, each moved on to the next event start, so no event is split.
|
|
const size_t nchunks = ThreadsForWork(n_parts, nthreads);
|
|
std::vector<std::vector<Reflection>> chunk_fulls(nchunks);
|
|
ParallelFor(static_cast<int>(nchunks), nthreads, [&](int c) {
|
|
size_t i = c * n_parts / nchunks;
|
|
size_t end = (c + 1) * n_parts / nchunks;
|
|
while (i > 0 && i < n_parts && continues(i)) ++i;
|
|
while (end < n_parts && continues(end)) ++end;
|
|
while (i < end) {
|
|
size_t j = i + 1;
|
|
while (j < n_parts && continues(j))
|
|
++j;
|
|
|
|
double sum_p = 0.0, sum_I = 0.0, sum_var = 0.0, sum_var_bkg = 0.0;
|
|
double p_corr = 0.0, p_qe = 0.0, p_flight = 0.0, p_frame = 0.0, p_x = 0.0, p_y = 0.0,
|
|
p_delta_phi = 0.0, p_zeta = 0.0, p_bkg = 0.0;
|
|
bool overloaded = false;
|
|
for (size_t m = i; m < j; ++m) {
|
|
const Reflection &r = outcomes[parts[m].o].reflections[parts[m].n];
|
|
overloaded = overloaded || r.overloaded;
|
|
const double p = r.partiality;
|
|
const float pc = r.prescaling_corr * r.qe_corr * r.flight_corr; // LP x QE x flight path
|
|
sum_p += p;
|
|
sum_I += static_cast<double>(r.I) * pc;
|
|
sum_var += static_cast<double>(r.sigma) * r.sigma * pc * pc;
|
|
sum_var_bkg += static_cast<double>(r.var_bkg) * pc * pc;
|
|
p_corr += p * pc;
|
|
p_qe += p * r.qe_corr;
|
|
p_flight += p * r.flight_corr;
|
|
p_frame += p * r.image_number;
|
|
p_x += p * DetectorX(r);
|
|
p_y += p * DetectorY(r);
|
|
p_delta_phi += p * r.delta_phi_deg;
|
|
p_zeta += p * r.zeta;
|
|
p_bkg += p * r.bkg;
|
|
}
|
|
Reflection full = outcomes[parts[i].o].reflections[parts[i].n];
|
|
i = j;
|
|
// An event with a saturated pixel is not a measurement; the merge drops it too.
|
|
if (overloaded || sum_p < min_partiality || sum_p < min_captured_fraction)
|
|
continue;
|
|
|
|
// The prescaling factor is applied by the writer, which multiplies I by the whole correction, so
|
|
// divide the event's own factor back out of the sums here. It is the same geometry for every part
|
|
// of one event to a median 2e-4, so the file's I / LP * QE * FLIGHT is still the raw count sum. The
|
|
// event's mean is taken on the whole correction and on the two path terms; the
|
|
// Lorentz-polarization part is what is left when those are divided out of it.
|
|
const float mean_corr = static_cast<float>(p_corr / sum_p);
|
|
full.qe_corr = static_cast<float>(p_qe / sum_p);
|
|
full.flight_corr = static_cast<float>(p_flight / sum_p);
|
|
full.prescaling_corr = mean_corr / (full.qe_corr * full.flight_corr);
|
|
full.I = static_cast<float>(sum_I / mean_corr);
|
|
full.sigma = static_cast<float>(std::sqrt(sum_var) / mean_corr);
|
|
full.var_bkg = static_cast<float>(sum_var_bkg / (static_cast<double>(mean_corr) * mean_corr));
|
|
full.partiality = static_cast<float>(sum_p);
|
|
full.image_number = static_cast<float>(p_frame / sum_p);
|
|
full.observed_x = static_cast<float>(p_x / sum_p);
|
|
full.observed_y = static_cast<float>(p_y / sum_p);
|
|
full.delta_phi_deg = static_cast<float>(p_delta_phi / sum_p);
|
|
full.zeta = static_cast<float>(p_zeta / sum_p);
|
|
full.bkg = static_cast<float>(p_bkg / sum_p);
|
|
chunk_fulls[c].push_back(full);
|
|
}
|
|
});
|
|
return chunk_fulls;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
MtzTable UnmergedMtzTable(const std::vector<IntegrationOutcome> &outcomes,
|
|
const UnitCell &unitCell,
|
|
const DiffractionExperiment &experiment,
|
|
bool sum_partials,
|
|
const OutputProvenance &provenance,
|
|
size_t nthreads,
|
|
const std::vector<ScaledFull> *scaled_fulls) {
|
|
const gemmi::SpaceGroup *spacegroup = &experiment.GetSpaceGroupOrP1();
|
|
MtzTable mtz;
|
|
mtz.space_group_name = spacegroup->xhm();
|
|
mtz.unit_cell = unitCell;
|
|
mtz.title = "Rugnux unmerged observations";
|
|
mtz.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction");
|
|
AddProvenanceHistory(mtz.history, provenance);
|
|
mtz.datasets = {{"HKL_base", "HKL_base", "HKL_base", 0.0f},
|
|
{"unmerged", experiment.GetSampleName(), "unmerged", experiment.GetWavelength_A()}};
|
|
const int dataset_id = 1;
|
|
|
|
// Every data column belongs to HKL_base and only the batches carry a dataset id, as in a
|
|
// POINTLESS file. Down to FLAG this is POINTLESS's own column set, plus QE and FLIGHT; the four after it
|
|
// are what rugnux measures beyond it - the offset of the reflection from the centre of its
|
|
// rocking curve, the Lorentz geometry of that curve, and the background that was subtracted.
|
|
//
|
|
// QE is not POINTLESS's, it is DIALS's: the detector's angle-dependent quantum efficiency, kept
|
|
// OUT of LP and given its own column because that is what both reference implementations do.
|
|
// LP here is Lorentz x polarization and nothing else, as XDS's RLP is (measured: XDS's column is
|
|
// flat to 0.1% across a detector over which the efficiency term spans 7%) and as DIALS's is
|
|
// (its LP comes from Corrections::lp(), lorentz/polarization, and it writes QE separately -
|
|
// a column of ones when it has no correction to report). Folding the sensor response into LP
|
|
// made this file assert something false about itself: a program that divided LP out to recover
|
|
// raw counts stripped the efficiency with it, and on an ordinary 13 keV geometry that term
|
|
// spans 13% end to end.
|
|
// Direction: QE is a DIVISOR, as in DIALS, and is normalised to 1 at normal incidence, so it
|
|
// rises above 1 toward the detector edge where a longer crossing makes the sensor more
|
|
// efficient. DIALS stores the un-normalised absorbed fraction instead; the two differ by the
|
|
// per-dataset constant QE(0), which is an overall scale.
|
|
for (const auto &[label, type] : std::vector<std::pair<const char *, char>>{
|
|
{"H", 'H'}, {"K", 'H'}, {"L", 'H'}, {"M/ISYM", 'Y'}, {"BATCH", 'B'}, {"I", 'J'}, {"SIGI", 'Q'},
|
|
{"FRACTIONCALC", 'R'}, {"XDET", 'R'}, {"YDET", 'R'}, {"ROT", 'R'}, {"LP", 'R'}, {"QE", 'R'},
|
|
{"FLIGHT", 'R'}, {"FLAG", 'I'}, {"DELPHI", 'R'}, {"ZETA", 'R'}, {"BGMEAN", 'R'}, {"BGVAR", 'R'}})
|
|
mtz.columns.push_back({label, type, 0, {}});
|
|
|
|
const auto gon = experiment.GetGoniometer();
|
|
// A scanning goniometer means rugnux integrated a rocking curve image by image, so each
|
|
// observation is one PART of a reflection unless the parts are summed here; on a still each is a
|
|
// whole measurement either way, and consecutive stills are different crystals, so there is
|
|
// nothing to sum there. Which of the two the file holds has to be declared in two places, because
|
|
// that is where POINTLESS and AIMLESS look for it: the batch header's data type, and the M flag
|
|
// that M/ISYM packs as 256*M + ISYM. Told nothing, they take each partial for a whole reflection
|
|
// and neither the symmetry determination nor the scaling survives it.
|
|
const bool scanning = gon && gon->IsScanning();
|
|
const bool partials = scanning && !sum_partials;
|
|
const float wedge_deg = gon ? gon->GetWedge_deg() : 0.0f;
|
|
const auto phi_start_deg = [&](float image_number) {
|
|
return gon ? gon->GetAngle_deg(image_number) : 0.0f;
|
|
};
|
|
|
|
// H K L are the ASU indices and M/ISYM says which symmetry operation (and which Friedel hand) got
|
|
// them there, so the index the reflection was actually measured at is recoverable - that is the
|
|
// crystal-frame information careless scales on, and what makes the file unmerged rather than a
|
|
// merge waiting to happen.
|
|
// I and SIGI are the integrated intensity with every deterministic per-reflection correction
|
|
// applied and nothing else, which is what IOBS means in every unmerged format. Three columns
|
|
// record what was applied, so the raw count sum is recoverable exactly as I / LP * QE * FLIGHT:
|
|
// LP is Lorentz x polarization, QE the sensor efficiency and FLIGHT the attenuation of the
|
|
// flight-path medium (both divisors, hence the multiplies).
|
|
// None of the three is a scale, and a program that reads this file has no way to recompute them. The
|
|
// partiality is NOT divided out - that IS a scale, FRACTIONCALC carries it, and every program
|
|
// this file is for wants to handle it its own way. All of this is of the unscaled file: given the merge's
|
|
// scaled fulls, I and SIGI are theirs (below), every correction and scale applied, and the columns
|
|
// beside them describe the event they were measured in.
|
|
gemmi::UnmergedHklMover hkl_mover(spacegroup);
|
|
// The batch headers are written for every image the observations span, not only for the images
|
|
// that produced one. AIMLESS starts a new run wherever the phi series jumps, so an image that
|
|
// indexed nothing leaves a hole that cuts the sweep in two: measured on a 360 deg sweep with 67
|
|
// unindexed images, AIMLESS made 23 runs of it and its scale model diverged - CC(1/2) 0.811,
|
|
// Rmeas overflowed to -1266, and the result no longer correlated with rugnux's own merge
|
|
// (Pearson 0.0018). Filling the gaps gives one run and the numbers the sweep deserves.
|
|
int first_batch = 0, last_batch = -1;
|
|
|
|
// Lattice-centring absences are integrated on purpose - prediction runs in P so the space-group
|
|
// search can confirm or disprove the centring - but they are not written here. POINTLESS reads
|
|
// the interstitial nodes, takes the lattice for primitive and demotes the space group, and they
|
|
// roughly halve AIMLESS's <I/sigma>. Screw and glide absences ARE written: they are the evidence
|
|
// the space group was chosen on, so a program reading this file can check that choice, and
|
|
// removing them turns its test into an assumption. XDS and DIALS draw the line in the same place.
|
|
const char centering = spacegroup->hm[0];
|
|
|
|
// On stills only the crystals the merge takes are written: one whose cell is not the run's was indexed
|
|
// in another lattice (on few spots, typically an axis multiple), so its indices do not name the
|
|
// reflections they would be read as. A rotation run has the one lattice.
|
|
const auto &indexing = experiment.GetIndexingSettings();
|
|
const auto merged_crystal = [&](const IntegrationOutcome &outcome) {
|
|
return scanning || outcome.latt.GetUnitCell().is_close(unitCell, indexing.GetUnitCellDistTolerance(),
|
|
indexing.GetUnitCellAngleTolerance_deg());
|
|
};
|
|
|
|
std::vector<std::vector<Reflection>> fulls;
|
|
if (scanning && sum_partials)
|
|
fulls = SumRockingEvents(outcomes,
|
|
experiment.GetScalingSettings().GetMinPartiality(),
|
|
experiment.GetScalingSettings().GetMinCapturedFraction(),
|
|
wedge_deg, nthreads);
|
|
std::vector<const Reflection *> rows;
|
|
const auto add_row = [&](const Reflection &r) {
|
|
if (systematic_absence(r.h, r.k, r.l, centering))
|
|
return;
|
|
// A summed full's image_number is its rocking-curve centroid, so this is the batch the
|
|
// reflection is centred on - which is what a batch means for a full everywhere else.
|
|
const int batch = 1 + static_cast<int>(std::lround(r.image_number));
|
|
if (last_batch < first_batch) { first_batch = batch; last_batch = batch; }
|
|
first_batch = std::min(first_batch, batch);
|
|
last_batch = std::max(last_batch, batch);
|
|
rows.push_back(&r);
|
|
};
|
|
if (scanning && sum_partials) {
|
|
for (const auto &piece : fulls)
|
|
for (const auto &r : piece)
|
|
add_row(r);
|
|
} else {
|
|
for (const auto &outcome : outcomes)
|
|
if (merged_crystal(outcome))
|
|
for (const auto &r : outcome.reflections)
|
|
add_row(r);
|
|
}
|
|
|
|
// The merge's own observations, where it handed them over: each event takes the full of its index whose
|
|
// peak frame is nearest its centroid (within the bridge of RockingEventFrameGap and a frame of rounding),
|
|
// and one the merge did not keep - an outlier, past the resolution cut - is not written. Unscaled
|
|
// intensities fitted again by a reading program's own scale model lose the merge's corrections; these
|
|
// are the intensities the merged file was made from.
|
|
std::vector<const ScaledFull *> scaled;
|
|
if (scanning && sum_partials && scaled_fulls != nullptr && !scaled_fulls->empty()) {
|
|
std::vector<const ScaledFull *> by_hkl(scaled_fulls->size());
|
|
for (size_t i = 0; i < by_hkl.size(); ++i)
|
|
by_hkl[i] = &(*scaled_fulls)[i];
|
|
const auto key = [](const ScaledFull *f) { return std::tie(f->h, f->k, f->l, f->image_number); };
|
|
ParallelSort(by_hkl.begin(), by_hkl.end(), nthreads,
|
|
[&](const ScaledFull *a, const ScaledFull *b) { return key(a) < key(b); });
|
|
const float max_gap = RockingEventFrameGap(wedge_deg) + 1.0f;
|
|
scaled.resize(rows.size(), nullptr);
|
|
ParallelChunks(static_cast<int>(rows.size()), ThreadsForWork(rows.size(), nthreads), [&](int lo, int hi) {
|
|
for (int i = lo; i < hi; ++i) {
|
|
const Reflection &r = *rows[i];
|
|
auto it = std::lower_bound(by_hkl.begin(), by_hkl.end(), r, [](const ScaledFull *f, const Reflection &x) {
|
|
return std::tie(f->h, f->k, f->l) < std::tie(x.h, x.k, x.l);
|
|
});
|
|
for (; it != by_hkl.end() && (*it)->h == r.h && (*it)->k == r.k && (*it)->l == r.l; ++it)
|
|
if (std::fabs((*it)->image_number - r.image_number) <= max_gap
|
|
&& (scaled[i] == nullptr
|
|
|| std::fabs((*it)->image_number - r.image_number)
|
|
< std::fabs(scaled[i]->image_number - r.image_number)))
|
|
scaled[i] = *it;
|
|
}
|
|
});
|
|
size_t kept = 0;
|
|
for (size_t i = 0; i < rows.size(); ++i)
|
|
if (scaled[i] != nullptr) {
|
|
rows[kept] = rows[i];
|
|
scaled[kept++] = scaled[i];
|
|
}
|
|
rows.resize(kept);
|
|
scaled.resize(kept);
|
|
mtz.title = "Rugnux unmerged observations, scaled";
|
|
}
|
|
|
|
// Every row has its own slot in each column, so they are filled on several workers.
|
|
const size_t w = mtz.columns.size();
|
|
for (auto &c : mtz.columns)
|
|
c.values.resize(rows.size());
|
|
ParallelChunks(static_cast<int>(rows.size()), ThreadsForWork(rows.size(), nthreads), [&](int lo, int hi) {
|
|
std::vector<float> row(w);
|
|
for (int i = lo; i < hi; ++i) {
|
|
const Reflection &r = *rows[i];
|
|
std::array<int, 3> hkl{r.h, r.k, r.l};
|
|
const int isym = hkl_mover.move_to_asu(hkl);
|
|
row[0] = static_cast<float>(hkl[0]);
|
|
row[1] = static_cast<float>(hkl[1]);
|
|
row[2] = static_cast<float>(hkl[2]);
|
|
row[3] = static_cast<float>((partials ? 256 : 0) + isym);
|
|
row[4] = static_cast<float>(1 + static_cast<int>(std::lround(r.image_number))); // BATCH, as above
|
|
const float corr = r.prescaling_corr * r.qe_corr * r.flight_corr;
|
|
row[5] = scaled.empty() ? r.I * corr : scaled[i]->I;
|
|
row[6] = scaled.empty() ? r.sigma * corr : scaled[i]->sigma;
|
|
row[7] = r.partiality;
|
|
row[8] = DetectorX(r);
|
|
row[9] = DetectorY(r);
|
|
row[10] = phi_start_deg(r.image_number) + wedge_deg / 2.0f;
|
|
// prescaling_corr is Lorentz x polarization, which is what LP means; qe_corr is the sensor
|
|
// efficiency and flight_corr the flight path beside it, each written in the divide-by
|
|
// direction, so that I / LP * QE * FLIGHT is the raw count. FLIGHT is <= 1 where QE is >= 1:
|
|
// the sensor makes an oblique reflection read high and the medium makes it read low.
|
|
row[11] = r.prescaling_corr;
|
|
row[12] = 1.0f / r.qe_corr;
|
|
row[13] = 1.0f / r.flight_corr;
|
|
row[14] = 0.0f; // FLAG: nothing here is a rejected observation
|
|
row[15] = r.delta_phi_deg;
|
|
row[16] = r.zeta;
|
|
row[17] = r.bkg;
|
|
row[18] = r.var_bkg;
|
|
for (size_t c = 0; c < w; c++)
|
|
mtz.columns[c].values[i] = row[c];
|
|
}
|
|
});
|
|
|
|
// The batch header's orientation matrix is the crystal at rotation angle zero - each batch's own
|
|
// PHISTT is applied on top of it - but the lattice stored with an outcome is the crystal as it
|
|
// stood on that image. Turn the first indexed one back by its own angle to get the orientation of
|
|
// the sweep, which is the one matrix POINTLESS also writes into every batch.
|
|
std::optional<CrystalLattice> lattice_at_zero;
|
|
for (const auto &outcome : outcomes) {
|
|
if (outcome.reflections.empty() || outcome.latt.CalcVolume() <= 1.0f || !merged_crystal(outcome))
|
|
continue;
|
|
const float mid_deg = phi_start_deg(outcome.reflections.front().image_number) + wedge_deg / 2.0f;
|
|
lattice_at_zero = gon ? outcome.latt.Multiply(gon->GetTransformationAngle(mid_deg)) : outcome.latt;
|
|
break;
|
|
}
|
|
|
|
// The batch header is written in the "Cambridge" frame - z along the principal rotation axis, x
|
|
// along the beam - while the jfjoch lab frame has the beam along +z, so the two are related by a
|
|
// rotation. These three lab-frame vectors are the Cambridge axes; a lab vector's components in
|
|
// that frame are its dot products with them. A still has no rotation axis, and any axis across
|
|
// the beam then defines the frame just as consistently.
|
|
// The axis is NEGATED: rugnux turns an observation made at angle phi back to phi = 0 by rotating
|
|
// it by +phi about the goniometer axis, so the crystal itself turns by -phi about it, while the
|
|
// MTZ batch header's scan axis is the one a batch's own increasing PHI turns the crystal about.
|
|
const Coord beam = experiment.GetDiffractionGeometry().GetScatteringVector().Normalize();
|
|
const Coord z_cam = gon ? -gon->GetAxis().Normalize() : Coord(0, 1, 0);
|
|
const Coord y_cam = (z_cam % beam).Normalize();
|
|
const Coord x_cam = (y_cam % z_cam).Normalize();
|
|
|
|
gemmi::Mtz::Batch batch;
|
|
batch.title = "Rugnux";
|
|
batch.axes.emplace_back("PHI");
|
|
batch.set_dataset_id(dataset_id);
|
|
batch.ints[12] = 1; // ncryst
|
|
batch.ints[14] = partials ? 1 : 2; // ldtype: oscillation data (2D spots) / area detector (3D)
|
|
batch.ints[15] = 1; // jsaxs: the goniostat scan axis
|
|
batch.ints[17] = 1; // ngonax
|
|
batch.ints[19] = 1; // ndet
|
|
batch.set_cell(unitCell);
|
|
if (lattice_at_zero) {
|
|
// Orientation matrix U, built from the reciprocal axes and stored column by column in
|
|
// Cambridge components, as gemmi's XDS_ASCII converter builds it (gemmi/xds2mtz.hpp).
|
|
const Coord ar = lattice_at_zero->Astar().Normalize();
|
|
const Coord cr = (ar % lattice_at_zero->Bstar()).Normalize();
|
|
const Coord u[3] = {ar, cr % ar, cr};
|
|
for (int i = 0; i < 3; i++) {
|
|
batch.floats[6 + 3 * i] = u[i] * x_cam;
|
|
batch.floats[7 + 3 * i] = u[i] * y_cam;
|
|
batch.floats[8 + 3 * i] = u[i] * z_cam;
|
|
}
|
|
}
|
|
batch.floats[40] = 1.0f; // scanax = [0, 0, 1]: the rotation axis IS z in the Cambridge frame
|
|
batch.floats[47] = wedge_deg;
|
|
batch.floats[61] = 1.0f; // e1 = scanax, the only goniostat axis
|
|
batch.floats[80] = -1.0f; // idealised source vector, antiparallel to the beam
|
|
batch.floats[83] = -(beam * x_cam); // s0, the source vector of the geometry as it really stands
|
|
batch.floats[84] = -(beam * y_cam);
|
|
batch.floats[85] = -(beam * z_cam);
|
|
batch.set_wavelength(experiment.GetWavelength_A());
|
|
batch.floats[111] = experiment.GetDetectorDistance_mm();
|
|
batch.floats[113] = 1.0f; // detector limits, in pixels
|
|
batch.floats[114] = static_cast<float>(experiment.GetXPixelsNum());
|
|
batch.floats[115] = 1.0f;
|
|
batch.floats[116] = static_cast<float>(experiment.GetYPixelsNum());
|
|
|
|
for (int number = first_batch; number <= last_batch; number++) {
|
|
batch.number = number;
|
|
batch.floats[36] = phi_start_deg(static_cast<float>(number - 1)); // phistt
|
|
batch.floats[37] = batch.floats[36] + wedge_deg; // phiend
|
|
mtz.batches.push_back({batch.number, batch.title, {batch.ints.begin(), batch.ints.end()}, batch.floats,
|
|
batch.axes});
|
|
}
|
|
|
|
// By H K L M/ISYM BATCH, the order POINTLESS leaves an unmerged file in. This is Mtz::sort(5) -
|
|
// a stable sort on those five columns - done on several workers: rows equal in all five keep the
|
|
// order they were added in, so with the row number as the last key the order is total and the
|
|
// parallel sort returns the same rows in the same order.
|
|
if (rows.empty())
|
|
gemmi::fail("No data."); // as gemmi's own sort has it
|
|
{
|
|
std::vector<int> order(rows.size());
|
|
std::iota(order.begin(), order.end(), 0);
|
|
ParallelSort(order.begin(), order.end(), nthreads, [&](int i, int j) {
|
|
for (int n = 0; n < 5; ++n) {
|
|
const float a = mtz.columns[n].values[i], b = mtz.columns[n].values[j];
|
|
if (a != b)
|
|
return a < b;
|
|
}
|
|
return i < j;
|
|
});
|
|
for (auto &c : mtz.columns) {
|
|
std::vector<float> sorted(order.size());
|
|
ParallelChunks(static_cast<int>(order.size()), ThreadsForWork(order.size(), nthreads), [&](int lo, int hi) {
|
|
for (int i = lo; i < hi; ++i)
|
|
sorted[i] = c.values[order[i]];
|
|
});
|
|
c.values.swap(sorted);
|
|
}
|
|
mtz.sort_order = {1, 2, 3, 4, 5};
|
|
}
|
|
return mtz;
|
|
}
|
|
|
|
void SetUnmergedMtzMosaicity(MtzTable &mtz, const std::vector<IntegrationOutcome> &outcomes) {
|
|
// The reflecting range of the same image the orientation matrix is taken from.
|
|
std::optional<float> mosaicity_deg;
|
|
for (const auto &outcome : outcomes) {
|
|
if (outcome.reflections.empty() || outcome.latt.CalcVolume() <= 1.0f)
|
|
continue;
|
|
mosaicity_deg = outcome.mosaicity_deg;
|
|
break;
|
|
}
|
|
for (auto &batch : mtz.batches)
|
|
batch.floats[21] = mosaicity_deg.value_or(0.0f); // crydat(0), the reflecting range
|
|
}
|
|
|
|
void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes,
|
|
const UnitCell &unitCell,
|
|
const DiffractionExperiment &experiment,
|
|
bool sum_partials,
|
|
const OutputProvenance &provenance,
|
|
const std::string &filename,
|
|
size_t nthreads,
|
|
const std::vector<ScaledFull> *scaled_fulls) {
|
|
MtzTable mtz = UnmergedMtzTable(outcomes, unitCell, experiment, sum_partials, provenance, nthreads, scaled_fulls);
|
|
SetUnmergedMtzMosaicity(mtz, outcomes);
|
|
WriteMtz(mtz, filename);
|
|
}
|
|
|
|
void WriteReflections(const std::vector<MergedReflection> &reflections,
|
|
const UnitCell &unitCell,
|
|
const DiffractionExperiment &experiment,
|
|
const MergeStatistics &statistics,
|
|
const ErrorModelReport &error_model,
|
|
const TwinningAnalysisResult &twinning,
|
|
const OutputProvenance &provenance,
|
|
const std::string &filename,
|
|
size_t nthreads,
|
|
const std::vector<ScaledFull> *scaled_fulls) {
|
|
// Write an MTZ, an mmCIF and a SHELX HKLF-4 .hkl - each has its uses downstream (MTZ for the CCP4 /
|
|
// phenix reflection tools, mmCIF for deposition and as the self-describing native format, HKLF-4 for
|
|
// SHELXL refinement and SHELXC / ANODE substructure solution). The .hkl holds the unmerged scaled
|
|
// fulls where the merge provides them (rotation), the merged reflections otherwise (stills).
|
|
WriteMtzReflections(reflections, unitCell, experiment, provenance, filename + ".mtz");
|
|
WriteMmcifReflections(reflections, unitCell, experiment, statistics, error_model, twinning,
|
|
provenance,
|
|
filename + ".cif", nthreads);
|
|
if (scaled_fulls && !scaled_fulls->empty())
|
|
WriteShelxHklReflections(*scaled_fulls, filename + ".hkl", nthreads);
|
|
else
|
|
WriteShelxHklReflections(reflections, experiment, filename + ".hkl", nthreads);
|
|
WriteShelxIns(unitCell, nullptr, experiment, std::filesystem::path(filename).filename().string(),
|
|
filename + ".ins");
|
|
}
|
|
|
|
MtzTable ObservationMtzTable(const std::vector<IntegrationOutcome> &outcomes,
|
|
const std::vector<std::vector<ObservationVerdict>> &verdicts,
|
|
std::optional<double> written_d_min,
|
|
const std::optional<StillsRunWidth> &run_width,
|
|
const UnitCell &unitCell,
|
|
const DiffractionExperiment &experiment) {
|
|
constexpr uint32_t OBS_BEYOND_DMIN = 4096, OBS_NO_VERDICT = 2048;
|
|
MtzTable mtz;
|
|
mtz.space_group_name = experiment.GetSpaceGroupOrP1().xhm();
|
|
mtz.unit_cell = unitCell;
|
|
mtz.title = "Rugnux stills observations (diagnostic)";
|
|
mtz.datasets = {{"HKL_base", "HKL_base", "HKL_base", 0.0f},
|
|
{"obs", experiment.GetSampleName(), "obs", experiment.GetWavelength_A()}};
|
|
if (run_width) {
|
|
std::string t = "RUNWIDTH gamma0 " + std::to_string(run_width->gamma0) + " gamma_e " + std::to_string(run_width->gamma_e)
|
|
+ " bw " + std::to_string(run_width->bandwidth_sigma) + " radial_scale "
|
|
+ std::to_string(run_width->radial_scale);
|
|
mtz.history.push_back(t);
|
|
std::string x = "RUNWIDTH_DSTAR", y = "RUNWIDTH_S";
|
|
for (size_t i = 0; i < run_width->table_dstar.size(); ++i) {
|
|
x += " " + std::to_string(run_width->table_dstar[i]);
|
|
y += " " + std::to_string(run_width->table_s[i]);
|
|
}
|
|
mtz.history.push_back(x);
|
|
mtz.history.push_back(y);
|
|
}
|
|
if (written_d_min)
|
|
mtz.history.push_back("WRITTEN_DMIN " + std::to_string(*written_d_min));
|
|
|
|
const std::vector<std::string> labels = {
|
|
"H", "K", "L", "FRAME", "IMG", "XPRED", "YPRED", "XOBS", "YOBS", "PANEL", "D",
|
|
"EPS", "EPSR", "EPSABS", "QX", "QY", "QZ", "PART", "ROCKSIG", "ZETA", "DPHI",
|
|
"I", "SIGI", "BKG", "VARBKG", "LP", "QE", "FLIGHT", "SCALECORR", "ISCALED", "SIGSCALED", "SIGMERGE",
|
|
"G", "CCIMG", "MOSAIC", "ASX", "ASY", "ASZ", "BSX", "BSY", "BSZ", "CSX", "CSY", "CSZ",
|
|
"CELLA", "CELLB", "CELLC", "CELLAL", "CELLBE", "CELLGA", "FLAGS", "USED"};
|
|
for (const auto &l : labels)
|
|
mtz.columns.push_back({l, (l == "H" || l == "K" || l == "L") ? 'H' : (l == "FLAGS" || l == "USED" ? 'I' : 'R'), 0, {}});
|
|
|
|
// Predictions are on the detector plane, so the panel is the one whose placed footprint holds the
|
|
// point (-1 for a detector without panels, -2 in a gap or off every panel).
|
|
const auto panel_of = [&](float x, float y) {
|
|
const auto p = experiment.GetModules().PlaneToPanel(x, y);
|
|
return p ? static_cast<float>(p->panel) : -2.0f;
|
|
};
|
|
// The run's rocking width at d*, as StillsPartialityRefine evaluates it: its table, linear in d* and held
|
|
// beyond the end shells, else the three-term formula.
|
|
const auto rock_sigma = [&](double dstar, double qz) {
|
|
if (!run_width)
|
|
return NAN;
|
|
const auto &x = run_width->table_dstar;
|
|
const auto &y = run_width->table_s;
|
|
if (!x.empty()) {
|
|
if (dstar <= x.front()) return static_cast<float>(y.front());
|
|
if (dstar >= x.back()) return static_cast<float>(y.back());
|
|
const size_t i = std::upper_bound(x.begin(), x.end(), dstar) - x.begin();
|
|
return static_cast<float>(y[i - 1] + (dstar - x[i - 1]) / (x[i] - x[i - 1]) * (y[i] - y[i - 1]));
|
|
}
|
|
const double ang = run_width->gamma_e * dstar, bw = run_width->bandwidth_sigma * qz;
|
|
return static_cast<float>(std::sqrt(run_width->gamma0 * run_width->gamma0 + ang * ang + bw * bw));
|
|
};
|
|
const double radial = run_width ? run_width->radial_scale : 1.0;
|
|
|
|
std::vector<float> row(labels.size());
|
|
const auto push = [&]() {
|
|
for (size_t c = 0; c < row.size(); ++c)
|
|
mtz.columns[c].values.push_back(row[c]);
|
|
};
|
|
for (size_t i = 0; i < outcomes.size(); ++i) {
|
|
const auto &o = outcomes[i];
|
|
if (o.reflections.empty() || o.latt.CalcVolume() <= 1.0f) // an unindexed frame: no lattice
|
|
continue;
|
|
const bool have_verdicts = i < verdicts.size() && verdicts[i].size() == o.reflections.size();
|
|
const Coord A = o.latt.Astar(), B = o.latt.Bstar(), C = o.latt.Cstar();
|
|
const Coord S0 = o.geom.GetScatteringVector();
|
|
const double inv_lambda = 1.0 / o.geom.GetWavelength_A();
|
|
const UnitCell cell = o.latt.GetUnitCell();
|
|
for (size_t j = 0; j < o.reflections.size(); ++j) {
|
|
const Reflection &r = o.reflections[j];
|
|
const Coord q = A * static_cast<float>(r.h) + B * static_cast<float>(r.k) + C * static_cast<float>(r.l);
|
|
const auto eps_at = [&](double s) {
|
|
const double x = s * q.x + S0.x, y = s * q.y + S0.y, z = s * q.z + S0.z;
|
|
return static_cast<float>(std::sqrt(x * x + y * y + z * z) - inv_lambda);
|
|
};
|
|
const double dstar = std::sqrt(q.x * q.x + q.y * q.y + q.z * q.z);
|
|
uint32_t flags = have_verdicts ? verdicts[i][j].flags : OBS_NO_VERDICT;
|
|
if (written_d_min && r.d < *written_d_min)
|
|
flags |= OBS_BEYOND_DMIN;
|
|
row = {static_cast<float>(r.h), static_cast<float>(r.k), static_cast<float>(r.l), r.image_number,
|
|
static_cast<float>(i), r.predicted_x, r.predicted_y, r.observed_x, r.observed_y,
|
|
panel_of(r.predicted_x, r.predicted_y), r.d,
|
|
eps_at(1.0), eps_at(radial), r.dist_ewald, q.x, q.y, q.z, r.partiality,
|
|
rock_sigma(radial * dstar, std::fabs(radial * q.z)), r.zeta, r.delta_phi_deg,
|
|
r.I, r.sigma, r.bkg, r.var_bkg, r.prescaling_corr, r.qe_corr, r.flight_corr,
|
|
r.image_scale_corr, r.I * r.image_scale_corr, r.sigma * r.image_scale_corr,
|
|
have_verdicts ? verdicts[i][j].sigma : NAN,
|
|
o.image_scale_g.value_or(NAN), o.image_scale_cc.value_or(NAN), o.mosaicity_deg.value_or(NAN),
|
|
A.x, A.y, A.z, B.x, B.y, B.z, C.x, C.y, C.z,
|
|
cell.a, cell.b, cell.c, cell.alpha, cell.beta, cell.gamma,
|
|
static_cast<float>(flags), flags == 0 ? 1.0f : 0.0f};
|
|
push();
|
|
}
|
|
}
|
|
return mtz;
|
|
}
|