mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 81 Wiki Activity

rugnux: .hkl holds the unmerged scaled fulls on rotation data

Browse Source 
The SHELX HKLF 4 file now has one record per full reflection - its partials summed,
the per-frame scale and every correction applied, sigma(I) as the merge weighted it -
at the index it was measured at, not averaged with its equivalents: the chemical
crystallographer's convention, so SHELXL computes Rint and Rsigma itself. Outliers
the merge rejected and fulls beyond its resolution cut are left out; no batch column
(it would select a BASF scale in SHELXL). The engine hands the fulls back only for a
merge that may be written (RotationScaleMerge::SetExportScaledFulls), so the search
merges, the pre-pass and the P1 cross-check carry no copy; the fulls follow the same
relabelling as the merged reflections (merge_to_written). Stills keep the merged
file. --mode scale writes the unmerged form too.

Validation: p.mtz md5 unchanged on myob/cytc/thau (GPU). SHELXL on the same runs,
merged-old vs unmerged-new (COD models, harness /data/tmp/sm_shared):
aspirin 20 keV Rint 0 -> 0.071, Rsigma 0.036 -> 0.043, R1 0.0964 -> 0.0969,
wR2 0.312 -> 0.310, GooF 1.53 -> 1.47; HEPES 20 keV Rint 0 -> 0.163, Rsigma 0.057
-> 0.068, R1 0.0903 -> 0.0897, wR2 0.318 -> 0.263, GooF 1.57 -> 1.10; the "input
data appear to be merged" warning is gone.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB
This commit is contained in:
leonarski_f
2026-10-04 01:16:50 +02:00
co-authored by Claude Opus 5.5
parent 359e749074
commit 9874db5d7f
12 changed files with 221 additions and 88 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/WriteReflections.cpp
+97 -65
View File
@@ -532,55 +532,35 @@ void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
mtz.write_to_file(filename);
}
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);
// SHELX HKLF 4 (SHELXC / ANODE input): fixed FORMAT(3I4,2F8.2), one record per reflection as
// h k l I sigma(I). 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.
// Intensities are put on a common scale so the largest value fits the F8.2 field (the absolute scale
// is irrelevant to SHELXC / ANODE, which use only ratios); I and sigma share the scale, so the
// anomalous signal is untouched. The file ends with a 0 0 0 terminator record.
const auto usable = [](float v, float s) { return std::isfinite(v) && std::isfinite(s) && s > 0.0f; };
double max_abs = 0.0;
for (const auto& r : rows) {
if (usable(r.Ip, r.sIp)) max_abs = std::max({max_abs, std::fabs(double(r.Ip)), double(r.sIp)});
if (usable(r.Im, r.sIm)) max_abs = std::max({max_abs, std::fabs(double(r.Im)), double(r.sIm)});
if (!usable(r.Ip, r.sIp) && !usable(r.Im, r.sIm) && usable(r.Imean, r.sImean))
max_abs = std::max({max_abs, std::fabs(double(r.Imean)), double(r.sImean)});
}
const double scale = (std::isfinite(max_abs) && max_abs > 0.0) ? 9999.0 / max_abs : 1.0;
std::ofstream out(filename);
if (!out)
throw std::runtime_error("WriteShelxHklReflections: cannot open " + filename);
// Up to two records per reflection, five formatted numbers each. Built in parallel into per-worker
// blocks and handed to the file in order, exactly as the mmCIF rows are; "%.2f" right-aligned in
// the fixed field is what `fixed` + `setprecision(2)` + `setw` made the stream write.
const auto column = [](std::string &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) {
char b[64];
const int n = std::snprintf(b, sizeof b, "%.2f", v);
return std::string(b, static_cast<size_t>(std::clamp(n, 0, static_cast<int>(sizeof b) - 1)));
};
const auto emit = [&column, &num2, scale](std::string &s, int h, int k, int l, float I, float sigma) {
column(s, std::to_string(h), 4);
column(s, std::to_string(k), 4);
column(s, std::to_string(l), 4);
column(s, num2(scale * I), 8);
column(s, num2(scale * sigma), 8);
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) {
char b[64];
const int n = std::snprintf(b, sizeof b, "%.2f", v);
return std::string(b, static_cast<size_t>(std::clamp(n, 0, static_cast<int>(sizeof b) - 1)));
};
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');
};
{
const size_t nrow = rows.size();
}
// 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);
@@ -589,24 +569,71 @@ void WriteShelxHklReflections(const std::vector<MergedReflection> &reflections,
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) {
const auto &r = rows[i];
const bool plus = usable(r.Ip, r.sIp);
const bool minus = usable(r.Im, r.sIm);
if (plus) emit(s, r.h, r.k, r.l, r.Ip, r.sIp);
if (minus) emit(s, -r.h, -r.k, -r.l, r.Im, r.sIm);
if (!plus && !minus && usable(r.Imean, r.sImean))
emit(s, r.h, r.k, r.l, r.Imean, r.sImean);
}
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()));
}
std::string tail;
emit(tail, 0, 0, 0, 0.0f, 0.0f); // HKLF-4 end-of-data marker
out.write(tail.data(), static_cast<std::streamsize>(tail.size()));
out.close();
// 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);
});
}
namespace {
@@ -1043,12 +1070,17 @@ void WriteReflections(const std::vector<MergedReflection> &reflections,
const ErrorModelReport &error_model,
const TwinningAnalysisResult &twinning,
const std::string &filename,
size_t nthreads) {
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 as
// the SHELXC / ANODE substructure-solution input).
// 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, filename + ".mtz");
WriteMmcifReflections(reflections, unitCell, experiment, statistics, error_model, twinning,
filename + ".cif", nthreads);
WriteShelxHklReflections(reflections, experiment, filename + ".hkl", nthreads);
if (scaled_fulls && !scaled_fulls->empty())
WriteShelxHklReflections(*scaled_fulls, filename + ".hkl", nthreads);
else
WriteShelxHklReflections(reflections, experiment, filename + ".hkl", nthreads);
}