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

rugnux: unit-cell e.s.d.s - ZERR in a new SHELX .ins beside the .hkl, _cell.*_esd in the mmCIF, CELL_ESD_* in the report

Browse Source 
Rotation runs now estimate the standard uncertainties of the cell. Every strong integrated
observation (I > 3 sigma, unclipped, unsaturated, off ice rings) gives |q| = 1/d at its observed
centroid, which depends on the detector geometry only; q^2 = h^T G* h is linear in the reciprocal
metric, so the components of G* the point group leaves free (projected from the group's rotation
parts, any setting) are fitted by weighted linear least squares with 4-sigma rejection, and the
GoF-scaled covariance is carried analytically onto a, b, c, alpha, beta, gamma (ties give equal
e.s.d.s, fixed angles zero). Least-squares cell e.s.d.s are known to be too small (Herbstein 2000),
so the same fit is made on 8 blocks of consecutive frames and their scatter / sqrt(8) is the second
estimate; the written e.s.d. is the larger of the two per parameter. Stills get none.

Outputs: <prefix>.ins next to every .hkl (TITL, CELL at the wavelength, ZERR with Z = order of the
group, LATT/SYMM from the space group - positive LATT only with -1 at the origin - a placeholder
organic SFAC/UNIT at 18 A^3 per non-H atom, HKLF 4, END); the mmCIF gains _cell.*_esd ("." where
fixed) and _cell_measurement.reflns_used/theta_min/theta_max, and writes the cell to 4/3 decimals;
report keys CELL_ESD, CELL_ESD_LSQ, CELL_ESD_BLOCKS, CELL_ESD_BLOCK_COUNT, CELL_ESD_REFLECTIONS,
CELL_ESD_THETA_RANGE_DEG (+ dev CELL_ESD_OBSERVATIONS, CELL_DSPACING_FIT).

On 13 small-molecule sets the block scatter is 10-40x the least-squares e.s.d. and is what gets
written; it lands within ~0.3-5x of XDS's CORRECT.LP e.s.d.s. SHELXT solves from the generated
.ins and SHELXL carries ZERR into its CIF su. Costs 12 ms on a small molecule, 0.16 s on 690k
lysozyme observations.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit e74a3be76261fe1b6c0ff7b80c21afc4768c46cd)
This commit is contained in:
leonarski_f committed 2026-10-26 14:33:38 +01:00
1 parent 22debb1be6
commit d8c14046b2
14 files changed
+805 -38
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

image_analysis/WriteReflections.cpp
+92 -10
View File
@@ -2,6 +2,7 @@
// 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"
@@ -19,6 +20,7 @@
#include <map>
#include <memory>
#include <tuple>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <sstream>
@@ -155,18 +157,38 @@ void AddProvenanceHistory(std::vector<std::string> &history, const OutputProvena
add("Collected: ", provenance.collection_date);
}
// _entry, _cell and _symmetry.
// 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 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, 3) << "\n";
out << "_cell.length_b " << Fmt(unitCell.b, 3) << "\n";
out << "_cell.length_c " << Fmt(unitCell.c, 3) << "\n";
out << "_cell.angle_alpha " << Fmt(unitCell.alpha, 2) << "\n";
out << "_cell.angle_beta " << Fmt(unitCell.beta, 2) << "\n";
out << "_cell.angle_gamma " << Fmt(unitCell.gamma, 2) << "\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";
@@ -270,7 +292,8 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
const TwinningAnalysisResult &twinning,
const OutputProvenance &provenance,
const std::string &filename,
size_t nthreads) {
size_t nthreads,
const CellEsd *cell_esd) {
std::ofstream out(filename);
if (!out)
@@ -287,7 +310,7 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
// ---------- data block ----------
out << "data_sample" << "\n";
out << "#\n";
WriteBlockHeader(out, entry_id, unitCell, sg);
WriteBlockHeader(out, entry_id, unitCell, sg, cell_esd);
// ---------- _audit ----------
out << "_audit.revision_id 1\n";
@@ -774,6 +797,63 @@ void WriteShelxHklReflections(const std::vector<ScaledFull> &fulls,
});
}
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);
char line[256];
out << "TITL " << title << " in " << sg.xhm() << "\n";
out << "REM Rugnux " << jfjoch_version() << ": CELL and ZERR from the data, SFAC/UNIT a placeholder"
" to replace with the formula\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's first number is Z, the formula units in the cell: one per general position is assumed.
if (esd) {
std::snprintf(line, sizeof(line), "ZERR %d %.5f %.5f %.5f %.4f %.4f %.4f\n", ops.order(), esd->esd[0],
esd->esd[1], esd->esd[2], esd->esd[3], esd->esd[4], esd->esd[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";
out << "HKLF 4\n";
out << "END\n";
}
namespace {
// The detector position an observation is written at: a reflection too weak to have a measured
@@ -1233,4 +1313,6 @@ void WriteReflections(const std::vector<MergedReflection> &reflections,
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");
}