diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index f198ef5e5..28cc8328f 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -74,12 +74,13 @@ that, or to lift the per-image loop's default ceiling of 16 workers per GPU. Nothing more is needed to pick the workflow: a dataset carrying a **goniometer axis** that turns is processed as a rotation sweep, one without (or with an axis that stands still) as **independent stills**, and scaling and merging run by default -in both. A rotation run that merges — the default — leaves eight files next to each other: +in both. A rotation run that merges — the default — leaves nine files next to each other: ``` myrun.mtz merged intensities + French-Wilson amplitudes, for CCP4 / phenix myrun.cif the same, as mmCIF - the self-describing format, and what to deposit myrun.hkl the scaled observations unmerged, as SHELX HKLF 4 - for SHELXL, SHELXC / ANODE +myrun.ins the SHELX instructions to go with it: CELL, ZERR (the cell's e.s.d.s), LATT, SYMM myrun_unmerged.mtz every observation before scaling, for pointless / aimless / careless - the largest file of the run (--no-export-unmerged skips it) myrun_report.txt what the run determined: cell, space group, statistics, warnings diff --git a/docs/RUGNUX_INTEGRATION.md b/docs/RUGNUX_INTEGRATION.md index a295624e9..0edab83f3 100644 --- a/docs/RUGNUX_INTEGRATION.md +++ b/docs/RUGNUX_INTEGRATION.md @@ -13,7 +13,10 @@ downstream suite running on Rugnux output. **mmCIF.** Standard items carry their standard meanings — `_refln.intensity_meas` / `_intensity_sigma`, the `pdbx_I_plus`/`pdbx_I_minus` and `pdbx_F_plus`/`pdbx_F_minus` anomalous pairs, `_reflns.*` and `_reflns_shell.*` for the merging statistics, `_reflns.B_iso_Wilson_estimate` for the Wilson B, and -`_cell.*` / `_diffrn_radiation_wavelength.wavelength` for the geometry. +`_cell.*` / `_diffrn_radiation_wavelength.wavelength` for the geometry. On rotation data `_cell.*_esd` +carries the cell's standard uncertainties (`.` for a parameter the symmetry fixes) and +`_cell_measurement.reflns_used` / `theta_min` / `theta_max` the reflections they were measured on (see +[Cell e.s.d.s](RUGNUX_REPORT.md#cell-esds)). Anything Rugnux reports that has no standard item is written under a **`jfjoch_`** prefix, inside the standard category it belongs to. That is a deliberate choice: a reader that does not know these items @@ -372,35 +375,24 @@ arrived without its `.mtz`. ## Small-molecule structures with SHELXT and SHELXL A small-molecule sweep (see [Small-molecule data](RUGNUX_TUTORIAL.md#small-molecule-data)) goes to -SHELX as **`xtal.hkl`**, the rotation run's scaled reflections unmerged, in HKLF 4. `shelxt` and -`shelxl` come with CCP4 as well as with the SHELX distribution. HKLF 4 carries no metadata, so the -instruction file is written by hand: the wavelength and the cell from the report, the lattice and -symmetry of the space group the report names, and the cell contents, which only the user knows. +SHELX as **`xtal.hkl`**, the rotation run's scaled reflections unmerged, in HKLF 4, with the instruction +file **`xtal.ins`** beside it. `shelxt` and `shelxl` come with CCP4 as well as with the SHELX +distribution. The `.ins` holds what the data determine: `CELL` at the wavelength, `ZERR` with the +cell's standard uncertainties (see [Cell e.s.d.s](RUGNUX_REPORT.md#cell-esds); its Z assumes one +formula unit per general position), and the `LATT` and `SYMM` of the space group the report names. +`SFAC` and `UNIT` are a placeholder - an organic composition at 18 ų per non-hydrogen atom - since +the cell contents are what only the user knows; edit them before solving: ``` cp xtal.hkl struct.hkl -w=$(grep '^WAVELENGTH= ' xtal_report.txt | cut -d' ' -f2) -cell=$(grep '^UNIT_CELL_CONSTANTS= ' xtal_report.txt | cut -d' ' -f2-) -cat > struct.ins < 3σ, +unclipped, unsaturated, off the ice rings) gives 1/d at its observed centroid, which depends on the +detector geometry but not on the crystal orientation, and q² = hᵀG\*h is linear in the reciprocal +metric G\*. The components of G\* the point group leaves free are fitted by least squares (one +e.s.d. on q for every observation, those more than 4σ off rejected), the covariance is scaled by the +goodness of fit and carried onto a, b, c, α, β, γ - **`CELL_ESD_LSQ=`**. Least-squares cell e.s.d.s +are known to be too small (F. H. Herbstein, "How precise are measurements of unit-cell dimensions +from single crystals?" (2000), Acta Cryst. B56, 547-557 +[doi:10.1107/S010876810000269X](https://doi.org/10.1107/S010876810000269X)), as they see the scatter of the centroids and nothing that +changes over the sweep, so the same fit is made on eight blocks of consecutive frames and the +scatter of their cells, over √blocks, is **`CELL_ESD_BLOCKS=`** (where at least four blocks could be +fitted; **`CELL_ESD_BLOCK_COUNT=`** says how many). `CELL_ESD=` is, per parameter, the larger of the +two. **`CELL_ESD_REFLECTIONS=`** is the number of distinct reflections used and +**`CELL_ESD_THETA_RANGE_DEG=`** their θ range, written as `_cell_measurement.reflns_used`, +`theta_min` and `theta_max` in the mmCIF. + +Neither estimate includes what is common to every observation - the wavelength and the detector +distance scale the whole cell - so the e.s.d. is a precision, not an accuracy. Stills runs have no +cell e.s.d.: the crystals' cells differ from one another by more than any one of them is measured to. + ## Mosaicity **`MOSAICITY_DEG=`** (section 2, and a `Mosaicity` line in the summary) is the mosaic spread the merge diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index fdfa111a8..3383f3c3f 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -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 #include #include +#include #include #include #include @@ -155,18 +157,38 @@ void AddProvenanceHistory(std::vector &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 &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 &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 &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(latt_code) + 1 : 1; + std::vector 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(std::lround(gemmi::UnitCell(cell).volume / 18.0))); + const int n_c = std::max(1, static_cast(std::lround(0.7 * nonh))); + const int n_n = std::max(1, static_cast(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(std::toupper(static_cast(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 &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"); } diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 1d9053fd8..b48b05018 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -14,6 +14,7 @@ struct MergeStatistics; struct TwinningAnalysisResult; +struct CellEsd; // geom_refinement/CellEsd.h // The error model as it is reported, already formatted. `isa` is the whole-range 1/sqrt(a*b), the // same quantity XDS's ISa denotes, so a file written here is directly comparable with a CORRECT.LP; @@ -36,8 +37,12 @@ struct OutputProvenance { std::string collection_date; // as the input states it (NXmx /entry/start_time), ISO-8601 }; +// An e.s.d. below this is a parameter the symmetry fixes: written "." in the mmCIF. +constexpr double CELL_ESD_FIXED = 1e-9; + // The merged reflections, their statistics and the experiment, as one PDBx/mmCIF data block. // nthreads: workers for the per-reflection row formatting, which is the bulk of the file. +// cell_esd: the cell's e.s.d.s (_cell.*_esd, _cell_measurement), where there are any. void WriteMmcifReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, @@ -46,7 +51,8 @@ void WriteMmcifReflections(const std::vector &reflections, const TwinningAnalysisResult &twinning, const OutputProvenance &provenance, const std::string &filename, - size_t nthreads); + size_t nthreads, + const CellEsd *cell_esd = nullptr); // The merged MTZ: H K L IMEAN SIGIMEAN [I(+) SIGI(+) I(-) SIGI(-)] F SIGF [F(+) SIGF(+) F(-) SIGF(-)] // FreeR_flag, one row per unique reflection in the CCP4 anomalous layout, the bracketed columns only where @@ -75,6 +81,12 @@ void WriteShelxHklReflections(const std::vector &fulls, const std::string &filename, size_t nthreads); +// The SHELX instruction file to go with the .hkl: TITL, CELL (at the wavelength), ZERR where the +// e.s.d.s are known, LATT, SYMM, a placeholder SFAC/UNIT, HKLF 4, END - what SHELXT reads, and what +// SHELXL's CIF takes the cell and its e.s.d.s from. +void WriteShelxIns(const UnitCell &cell, const CellEsd *esd, const DiffractionExperiment &experiment, + const std::string &title, const std::string &filename); + // Unmerged observations in the column and batch-header layout POINTLESS writes: aimless, pointless, // careless and iotbx.merging_statistics all read that layout. H K L are the ASU indices and M/ISYM // recovers the index the reflection was measured at (which is what careless needs to see the crystal diff --git a/image_analysis/geom_refinement/CMakeLists.txt b/image_analysis/geom_refinement/CMakeLists.txt index c655cf25f..397034851 100644 --- a/image_analysis/geom_refinement/CMakeLists.txt +++ b/image_analysis/geom_refinement/CMakeLists.txt @@ -34,6 +34,8 @@ ADD_LIBRARY(JFJochGeomRefinement STATIC GeometryRefiner.h LatticeReduction.cpp LatticeReduction.h + CellEsd.cpp + CellEsd.h ) # fftw3f: BeamCenterFFT scores every candidate centre with one transform set (always fetched, diff --git a/image_analysis/geom_refinement/CellEsd.cpp b/image_analysis/geom_refinement/CellEsd.cpp new file mode 100644 index 000000000..23340f18b --- /dev/null +++ b/image_analysis/geom_refinement/CellEsd.cpp @@ -0,0 +1,257 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "CellEsd.h" + +#include +#include + +#include "../../common/ParallelFor.h" + +namespace { + using Vec6 = Eigen::Matrix; + // At most six parameters: fixed-capacity, so the per-observation products never touch the heap. + using VecK = Eigen::Matrix; + using MatK = Eigen::Matrix; + + // A symmetric matrix as (g11, g22, g33, g12, g13, g23) and back. + Vec6 ToVec6(const Eigen::Matrix3d &m) { + Vec6 v; + v << m(0, 0), m(1, 1), m(2, 2), m(0, 1), m(0, 2), m(1, 2); + return v; + } + + Eigen::Matrix3d FromVec6(const Vec6 &v) { + Eigen::Matrix3d m; + m << v(0), v(3), v(4), + v(3), v(1), v(5), + v(4), v(5), v(2); + return m; + } + + // h^T G h = hh . (g11, g22, g33, g12, g13, g23) + Vec6 IndexProducts(const CellObservation &o) { + const double h = o.h, k = o.k, l = o.l; + Vec6 v; + v << h * h, k * k, l * l, 2 * h * k, 2 * h * l, 2 * k * l; + return v; + } + + double ToDeg(double rad) { return rad * 180.0 / M_PI; } +} + +std::vector MetricBasis(const gemmi::SpaceGroup &sg) { + // Project each of the six symmetric unit matrices onto the metrics the group leaves invariant, + // averaging R G R^T over the rotation parts, and keep an orthonormal basis of what comes out. + const gemmi::GroupOps ops = sg.operations(); + std::vector rots; + for (const gemmi::Op &op : ops.sym_ops) { + Eigen::Matrix3d r; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + r(i, j) = static_cast(op.rot[i][j]) / gemmi::Op::DEN; + rots.push_back(r); + } + Eigen::Matrix projected; + for (int m = 0; m < 6; m++) { + Vec6 e = Vec6::Zero(); + e(m) = 1.0; + const Eigen::Matrix3d g = FromVec6(e); + Eigen::Matrix3d sum = Eigen::Matrix3d::Zero(); + for (const auto &r : rots) + sum += r * g * r.transpose(); + projected.col(m) = ToVec6(sum / static_cast(rots.size())); + } + Eigen::JacobiSVD> svd(projected, Eigen::ComputeFullU); + std::vector basis; + for (int m = 0; m < 6; m++) + if (svd.singularValues()(m) > 1e-8 * svd.singularValues()(0)) + basis.push_back(FromVec6(svd.matrixU().col(m))); + return basis; +} + +std::array CellFromReciprocalMetric(const Eigen::Matrix3d &gstar) { + const Eigen::Matrix3d g = gstar.inverse(); + const double a = std::sqrt(g(0, 0)), b = std::sqrt(g(1, 1)), c = std::sqrt(g(2, 2)); + return {a, b, c, ToDeg(std::acos(g(1, 2) / (b * c))), ToDeg(std::acos(g(0, 2) / (a * c))), + ToDeg(std::acos(g(0, 1) / (a * b)))}; +} + +std::array CellEsdFromCovariance(const std::vector &basis, const Eigen::VectorXd &p, + const Eigen::MatrixXd &cov) { + Eigen::Matrix3d gstar = Eigen::Matrix3d::Zero(); + for (size_t k = 0; k < basis.size(); k++) + gstar += p(k) * basis[k]; + const Eigen::Matrix3d g = gstar.inverse(); + const double len[3] = {std::sqrt(g(0, 0)), std::sqrt(g(1, 1)), std::sqrt(g(2, 2))}; + // The angle opposite axis i lies between the other two, j and l: cos = G_jl / (len_j len_l). + const int pair[3][2] = {{1, 2}, {0, 2}, {0, 1}}; + + // Analytic derivatives of the six parameters with respect to each p_k: dG = -G dG* G. + Eigen::MatrixXd jac(6, static_cast(basis.size())); + for (size_t k = 0; k < basis.size(); k++) { + const Eigen::Matrix3d dg = -g * basis[k] * g; + double dlen[3]; + for (int i = 0; i < 3; i++) { + dlen[i] = dg(i, i) / (2.0 * len[i]); + jac(i, k) = dlen[i]; + } + for (int i = 0; i < 3; i++) { + const int j = pair[i][0], l = pair[i][1]; + const double cosv = g(j, l) / (len[j] * len[l]); + const double dcos = dg(j, l) / (len[j] * len[l]) - cosv * (dlen[j] / len[j] + dlen[l] / len[l]); + jac(3 + i, k) = ToDeg(-dcos / std::sqrt(1.0 - cosv * cosv)); + } + } + const Eigen::MatrixXd cell_cov = jac * cov * jac.transpose(); + std::array esd{}; + // A parameter the symmetry fixes comes out at rounding level (1e-19), not at zero: make it zero. + for (int i = 0; i < 6; i++) { + esd[i] = std::sqrt(std::max(0.0, cell_cov(i, i))); + if (esd[i] < 1e-12) + esd[i] = 0.0; + } + return esd; +} + +MetricFit FitReciprocalMetric(std::span obs, const std::vector &basis) { + const int nparam = static_cast(basis.size()); + Eigen::Matrix v(6, nparam); + for (int k = 0; k < nparam; k++) + v.col(k) = ToVec6(basis[k]); + + MetricFit fit; + fit.used.assign(obs.size(), true); + MatK normal; + double s2 = 0.0; + // Fit, reject the observations further than 4 sigma off it, fit again - until nothing changes. + for (int iter = 0; iter < 5; iter++) { + normal = MatK::Zero(nparam, nparam); + VecK rhs = VecK::Zero(nparam); + int64_t n = 0; + for (size_t i = 0; i < obs.size(); i++) { + if (!fit.used[i]) + continue; + const double q2 = static_cast(obs[i].q) * obs[i].q; + const VecK a = v.transpose() * IndexProducts(obs[i]); + normal.noalias() += a * a.transpose() / q2; // weight 1/q^2 on q^2: one e.s.d. on q + rhs.noalias() += a; // a * q^2 / q^2 + n++; + } + if (n < 10 * nparam + 10) + return {}; + fit.p = normal.ldlt().solve(rhs); + double wr2 = 0.0; + std::vector z(obs.size()); + for (size_t i = 0; i < obs.size(); i++) { + const double q2 = static_cast(obs[i].q) * obs[i].q; + const VecK a = v.transpose() * IndexProducts(obs[i]); + const double r = q2 - a.dot(fit.p); + z[i] = r * r / q2; + if (fit.used[i]) + wr2 += z[i]; + } + s2 = wr2 / static_cast(n - nparam); + bool changed = false; + for (size_t i = 0; i < obs.size(); i++) { + const bool keep = z[i] <= 16.0 * s2; + changed |= keep != fit.used[i]; + fit.used[i] = keep; + } + if (!changed) + break; + } + fit.cov = s2 * Eigen::MatrixXd(normal).inverse(); + fit.ok = fit.p.allFinite() && fit.cov.allFinite(); + return fit; +} + +std::optional EstimateCellEsd(std::span obs, const gemmi::SpaceGroup &sg, + double wavelength_A, int nblocks) { + const auto basis = MetricBasis(sg); + const MetricFit fit = FitReciprocalMetric(obs, basis); + if (!fit.ok) + return {}; + const auto metric = [&basis](const Eigen::VectorXd &p) { + Eigen::Matrix3d g = Eigen::Matrix3d::Zero(); + for (size_t k = 0; k < basis.size(); k++) + g += p(k) * basis[k]; + return g; + }; + + CellEsd out; + out.esd_lsq = CellEsdFromCovariance(basis, fit.p, fit.cov); + const auto cell = CellFromReciprocalMetric(metric(fit.p)); + out.fitted = UnitCell{static_cast(cell[0]), static_cast(cell[1]), static_cast(cell[2]), + static_cast(cell[3]), static_cast(cell[4]), static_cast(cell[5])}; + + std::vector keys; + double q_min = INFINITY, q_max = 0.0; + int32_t frame_min = INT32_MAX, frame_max = INT32_MIN; + for (size_t i = 0; i < obs.size(); i++) { + if (!fit.used[i]) + continue; + const auto &o = obs[i]; + keys.push_back(((static_cast(o.h) + 1024) << 22) | ((static_cast(o.k) + 1024) << 11) + | (static_cast(o.l) + 1024)); + q_min = std::min(q_min, static_cast(o.q)); + q_max = std::max(q_max, static_cast(o.q)); + frame_min = std::min(frame_min, o.frame); + frame_max = std::max(frame_max, o.frame); + } + out.observations = static_cast(keys.size()); + std::sort(keys.begin(), keys.end()); + out.reflections = std::unique(keys.begin(), keys.end()) - keys.begin(); + out.theta_min_deg = ToDeg(std::asin(std::min(1.0, wavelength_A * q_min / 2.0))); + out.theta_max_deg = ToDeg(std::asin(std::min(1.0, wavelength_A * q_max / 2.0))); + + // The same fit on blocks of consecutive frames: the scatter of their cells, over sqrt(blocks). + std::vector> block_cells; + const int64_t span = static_cast(frame_max) - frame_min + 1; + std::vector> block_obs(nblocks); + for (size_t i = 0; i < obs.size(); i++) + if (fit.used[i]) + block_obs[(obs[i].frame - frame_min) * nblocks / span].push_back(obs[i]); + for (const auto &b : block_obs) { + const MetricFit bf = FitReciprocalMetric(b, basis); + if (bf.ok) + block_cells.push_back(CellFromReciprocalMetric(metric(bf.p))); + } + out.blocks = static_cast(block_cells.size()); + out.esd = out.esd_lsq; + if (out.blocks >= 4) { + std::array esd{}; + for (int j = 0; j < 6; j++) { + double mean = 0.0; + for (const auto &c : block_cells) + mean += c[j] / out.blocks; + double var = 0.0; + for (const auto &c : block_cells) + var += (c[j] - mean) * (c[j] - mean) / (out.blocks - 1); + esd[j] = std::sqrt(var / out.blocks); + out.esd[j] = std::max(out.esd[j], esd[j]); + } + out.esd_blocks = esd; + } + return out; +} + +std::vector GatherCellObservations(const std::vector &outcomes, + size_t nthreads) { + std::vector> per_frame(outcomes.size()); + ParallelFor(static_cast(outcomes.size()), nthreads, [&](int f) { + const IntegrationOutcome &out = outcomes[f]; + for (const auto &r : out.reflections) { + if (!r.observed || r.clipped || r.overloaded || r.on_ice_ring || !(r.sigma > 0.0f) + || !(r.I > 3.0f * r.sigma) || !std::isfinite(r.observed_x) || !std::isfinite(r.observed_y)) + continue; + const float q = out.geom.DetectorToRecip(r.observed_x, r.observed_y).Length(); + if (q > 0.0f) + per_frame[f].push_back(CellObservation{r.h, r.k, r.l, q, f}); + } + }); + std::vector obs; + for (const auto &v : per_frame) + obs.insert(obs.end(), v.begin(), v.end()); + return obs; +} diff --git a/image_analysis/geom_refinement/CellEsd.h b/image_analysis/geom_refinement/CellEsd.h new file mode 100644 index 000000000..8a112d290 --- /dev/null +++ b/image_analysis/geom_refinement/CellEsd.h @@ -0,0 +1,84 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include +#include + +#include + +#include "../../common/UnitCell.h" +#include "../IntegrationOutcome.h" +#include "gemmi/symmetry.hpp" + +// The e.s.d.s of the unit cell, as SHELXL's ZERR and a CIF's _cell.*_esd ask for them. +// +// The cell is fitted to d-spacings: each strong integrated observation gives |q| = 1/d at its +// observed centroid, which depends on the detector geometry alone (not on the orientation, nor on the +// rotation angle), and q^2 = h^T G* h is LINEAR in the reciprocal metric tensor G*. So the fit is a +// plain linear least squares in the components of G* the point group leaves free, and its covariance, +// scaled by the goodness of fit, maps onto a, b, c, alpha, beta, gamma through the derivatives of the +// cell with respect to G*. A length the group ties to another gets the same e.s.d., an angle it fixes +// gets zero. +// +// Least-squares cell e.s.d.s are known to be too small (Herbstein (2000) Acta Cryst. B56, 547-557): +// they see only the scatter of the centroids, not what changes over the sweep. The same fit is +// therefore also made on blocks of consecutive frames, and the scatter of the block cells, over the +// square root of their number, is the second estimate. The e.s.d. written is the larger of the two. + +// One observation: the Miller index in the setting of the cell, and |q| = 1/d at the observed +// centroid, in 1/A. +struct CellObservation { + int32_t h = 0, k = 0, l = 0; + float q = 0.0f; + int32_t frame = 0; +}; + +struct CellEsd { + std::array esd{}; // written: per parameter the larger of the two below + std::array esd_lsq{}; // least squares, scaled by the goodness of fit + std::optional> esd_blocks; // block scatter / sqrt(blocks) + int blocks = 0; + int64_t observations = 0; // observations in the fit, after rejection + int64_t reflections = 0; // distinct h k l among them (_cell_measurement.reflns_used) + double theta_min_deg = 0.0; + double theta_max_deg = 0.0; + UnitCell fitted{}; // the d-spacing fit's own cell, for comparison with the reported one +}; + +// Symmetric 3x3 matrices spanning the reciprocal metrics G* the point group of `sg` allows +// (G* = R G* R^T for every rotation part R): 6 for triclinic down to 1 for cubic. Orthonormal as +// 6-vectors (g11, g22, g33, g12, g13, g23). +std::vector MetricBasis(const gemmi::SpaceGroup &sg); + +// a, b, c (A) and alpha, beta, gamma (deg) of the reciprocal metric G*. +std::array CellFromReciprocalMetric(const Eigen::Matrix3d &gstar); + +// The e.s.d.s of a, b, c, alpha, beta, gamma for G* = sum_k p_k basis_k with covariance cov(p). +std::array CellEsdFromCovariance(const std::vector &basis, const Eigen::VectorXd &p, + const Eigen::MatrixXd &cov); + +struct MetricFit { + bool ok = false; + Eigen::VectorXd p; + Eigen::MatrixXd cov; // scaled by the goodness of fit + std::vector used; +}; + +// Weighted linear least squares of q^2 = h^T G* h, |q| taken at one constant e.s.d. (weight 1/q^2 on +// q^2), with the observations further than 4 sigma off the fit rejected and the fit repeated. +MetricFit FitReciprocalMetric(std::span obs, const std::vector &basis); + +// Both estimates for observations in the setting of `sg`; blocks are cut by frame. Empty with too +// few observations to fit. +std::optional EstimateCellEsd(std::span obs, const gemmi::SpaceGroup &sg, + double wavelength_A, int nblocks = 8); + +// The observations of a rotation run's integrated frames: observed, unclipped, unsaturated, off the +// ice rings and I > 3 sigma. +std::vector GatherCellObservations(const std::vector &outcomes, + size_t nthreads); diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 06fe848e8..0596612a9 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -886,6 +886,24 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, } if (result.consensus_cell.has_value()) Add(s, KeyText("UNIT_CELL_CONSTANTS", CellString(*result.consensus_cell))); + // The cell's e.s.d.s as XDS's E.S.D. OF CELL PARAMETERS (CellEsd.h): the written one (the .ins's + // ZERR, the mmCIF's _cell.*_esd) is per parameter the larger of the least-squares e.s.d. and the + // scatter of the cells fitted on blocks of frames. + if (result.cell_esd.has_value()) { + const CellEsd &e = *result.cell_esd; + const auto esd_string = [](const std::array &v) { + return fmt::format("{:.2E} {:.2E} {:.2E} {:.2E} {:.2E} {:.2E}", v[0], v[1], v[2], v[3], v[4], v[5]); + }; + Add(s, KeyText("CELL_ESD", esd_string(e.esd))); + Add(s, KeyText("CELL_ESD_LSQ", esd_string(e.esd_lsq))); + if (e.esd_blocks) + Add(s, KeyText("CELL_ESD_BLOCKS", esd_string(*e.esd_blocks))); + Add(s, KeyInt("CELL_ESD_BLOCK_COUNT", e.blocks)); + Add(s, KeyInt("CELL_ESD_REFLECTIONS", e.reflections)); + Add(s, KeyInt("CELL_ESD_OBSERVATIONS", e.observations, /*dev=*/true)); + Add(s, KeyText("CELL_ESD_THETA_RANGE_DEG", fmt::format("{:.2f} {:.2f}", e.theta_min_deg, e.theta_max_deg))); + Add(s, KeyText("CELL_DSPACING_FIT", CellString(e.fitted), /*dev=*/true)); + } // Kabsch's sigma_M, the Gaussian e.s.d. of the rocking width in degrees - XDS's // REFLECTING_RANGE_E.S.D. - that the merge computed partiality from. if (std::isfinite(result.mosaicity_deg)) { diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 0e4501b1f..2389bc1ac 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -33,6 +33,7 @@ #include "../image_analysis/scale_merge/StillsPartialityRefine.h" // StillsPartialityRefine::RunWidth #include "../image_analysis/scale_merge/TranslationalNCS.h" // TranslationalNCSResult #include "../image_analysis/indexing/HarmonicContamination.h" // HarmonicContaminationResult +#include "../image_analysis/geom_refinement/CellEsd.h" // CellEsd #include "../image_analysis/structure_refinement/ModelValidation.h" // ModelValidationResult #include "../image_analysis/structure_refinement/MolecularReplacement.h" // MolecularReplacementResult #include "../image_analysis/structure_refinement/Refinement.h" // RefinementResult @@ -393,6 +394,7 @@ struct MergedOutput { ErrorModelReport error_model; TwinningAnalysisResult twinning; // a twin hypothesis' law and fraction, for its mmCIF std::vector scaled_fulls; // the .hkl's records, where the merge made them (rotation) + std::optional cell_esd; // the cell's e.s.d.s: the mmCIF's _cell.*_esd and the .ins's ZERR }; // An atomic model as placed against the data (_model.cif, WritePlacedModel): `cell` and @@ -482,6 +484,8 @@ struct ProcessResult { // file's. The report's XDS keyword block reads it. Empty on stills. std::optional used_rotation_axis; std::optional consensus_cell; + // The e.s.d.s of consensus_cell (CellEsd.h); rotation only. + std::optional cell_esd; // The centring `consensus_cell` is written in; empty when nothing named one (the stills path // averages per-image cells and has no lattice search behind them). A conventional cell says // nothing about size on its own - a centred one is an exact integer multiple of its primitive diff --git a/rugnux/RugnuxOutputFiles.cpp b/rugnux/RugnuxOutputFiles.cpp index f5650f5c8..121c67270 100644 --- a/rugnux/RugnuxOutputFiles.cpp +++ b/rugnux/RugnuxOutputFiles.cpp @@ -137,20 +137,24 @@ void WriteRugnuxOutputs(const ProcessResult &result, const std::string &prefix, void WriteMergedOutput(const MergedOutput &merged, const std::string &prefix, const OutputProvenance &provenance, size_t nthreads, bool hkl) { - // An MTZ, an mmCIF and a SHELX HKLF-4 .hkl - each has its uses downstream (MTZ for the CCP4 / phenix + // An MTZ, an mmCIF and a SHELX HKLF-4 .hkl with its .ins - 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). const std::string file = prefix + merged.suffix; WriteMtz(merged.mtz, file + ".mtz"); + const CellEsd *esd = merged.cell_esd ? &*merged.cell_esd : nullptr; WriteMmcifReflections(merged.reflections, merged.cell, merged.experiment, merged.statistics, merged.error_model, - merged.twinning, provenance, file + ".cif", nthreads); + merged.twinning, provenance, file + ".cif", nthreads, esd); if (!hkl) return; if (!merged.scaled_fulls.empty()) WriteShelxHklReflections(merged.scaled_fulls, file + ".hkl", nthreads); else WriteShelxHklReflections(merged.reflections, merged.experiment, file + ".hkl", nthreads); + // The SHELX instructions beside it, with the cell's e.s.d.s as ZERR where the run measured them. + WriteShelxIns(merged.cell, esd, merged.experiment, std::filesystem::path(file).filename().string(), + file + ".ins"); } void WriteModelMaps(const ModelMaps &maps, const std::string &prefix, Logger &logger) { diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 3e5c3ca89..4d9eee3c0 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -4826,6 +4826,23 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { } } + // The cell's e.s.d.s (CellEsd.h), from the same integrated observations and at the same point as the + // unmerged MTZ below: the indexing and the group are final here. Rotation only - a stills run keeps no + // per-image reflections to the end, and its crystals' cells differ by more than any one is measured to. + if (is_rotation && !geometry_prepass && !superseded && result.consensus_cell) { + TimingMark mark(logger, "cell e.s.d."); + const auto obs = GatherCellObservations(indexer->GetIntegrationOutcome(), config_.nthreads); + result.cell_esd = EstimateCellEsd(obs, experiment_.GetSpaceGroupOrP1(), experiment_.GetWavelength_A()); + if (result.cell_esd) { + const CellEsd &e = *result.cell_esd; + logger.Info("Cell e.s.d.s from {} observations of {} reflections: {:.5f} {:.5f} {:.5f} {:.4f} {:.4f} {:.4f}" + " (least squares {:.5f} {:.5f} {:.5f} {:.4f} {:.4f} {:.4f}, {} frame blocks)", + e.observations, e.reflections, e.esd[0], e.esd[1], e.esd[2], e.esd[3], e.esd[4], e.esd[5], + e.esd_lsq[0], e.esd_lsq[1], e.esd_lsq[2], e.esd_lsq[3], e.esd_lsq[4], e.esd_lsq[5], + e.blocks); + } + } + // The unmerged MTZ (kept with the other outputs, see there) is built from here on, beside the // reflection files: the indexing and the group it is written in are final once the model has been // read, and nothing after this changes the integrated observations it reads. It gets its own copy of @@ -4867,6 +4884,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { config_.provenance), sm.merged, *result.consensus_cell, experiment_, sm.statistics, em_report, result.twinning, std::move(sm.scaled_fulls)}; + result.outputs->merged->cell_esd = result.cell_esd; if (keep_merged_reflections_) result.merged_reflections = sm.merged; // The twin hypothesis beside it (see prepare_twin_file): the reflection files only, with the diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index e1e3288bb..2f50cb2ab 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -137,6 +137,7 @@ ADD_EXECUTABLE(jfjoch_test ResultReportTest.cpp DiagnosticOutputTest.cpp PostRefineTest.cpp + CellEsdTest.cpp RotationScaleWalkTest.cpp AndersonMixingTest.cpp RugnuxLargeTest.cpp diff --git a/tests/CellEsdTest.cpp b/tests/CellEsdTest.cpp new file mode 100644 index 000000000..f993a2149 --- /dev/null +++ b/tests/CellEsdTest.cpp @@ -0,0 +1,264 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include +#include +#include + +#include + +#include "../image_analysis/geom_refinement/CellEsd.h" +#include "../image_analysis/WriteReflections.h" +#include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/TwinningAnalysis.h" + +namespace { + Eigen::Matrix3d ReciprocalMetric(double a, double b, double c, double al, double be, double ga) { + const gemmi::UnitCell cell(a, b, c, al, be, ga); + Eigen::Matrix3d g; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + g(i, j) = cell.frac.mat[i][0] * cell.frac.mat[j][0] + cell.frac.mat[i][1] * cell.frac.mat[j][1] + + cell.frac.mat[i][2] * cell.frac.mat[j][2]; // a*_i . a*_j, a*_i a row + return g; + } + + // The basis coordinates of G* (the basis is orthonormal over g11, g22, g33, g12, g13, g23), and + // how far G* lies outside the basis' span. + Eigen::VectorXd Coordinates(const std::vector &basis, const Eigen::Matrix3d &g, + double &residual) { + const auto vec = [](const Eigen::Matrix3d &m) { + Eigen::Matrix v; + v << m(0, 0), m(1, 1), m(2, 2), m(0, 1), m(0, 2), m(1, 2); + return v; + }; + Eigen::VectorXd p(basis.size()); + Eigen::Matrix rest = vec(g); + for (size_t k = 0; k < basis.size(); k++) { + p(k) = vec(basis[k]).dot(vec(g)); + rest -= p(k) * vec(basis[k]); + } + residual = rest.norm() / vec(g).norm(); + return p; + } + + const gemmi::SpaceGroup &Sg(const char *name) { return *gemmi::find_spacegroup_by_name(name); } + + // The e.s.d.s of the cell under a covariance of the metric's coordinates, against the scatter of + // the cells of metrics drawn from it. + void CheckAgainstSampling(const char *sg, double a, double b, double c, double al, double be, double ga) { + const auto basis = MetricBasis(Sg(sg)); + double residual = 0.0; + const Eigen::VectorXd p0 = Coordinates(basis, ReciprocalMetric(a, b, c, al, be, ga), residual); + REQUIRE(residual < 1e-12); + const int n = static_cast(basis.size()); + Eigen::MatrixXd cov = Eigen::MatrixXd::Zero(n, n); + for (int k = 0; k < n; k++) + cov(k, k) = std::pow(1e-5 * p0.norm() * (k + 1), 2); + if (n > 1) + cov(0, 1) = cov(1, 0) = 0.3 * std::sqrt(cov(0, 0) * cov(1, 1)); + const auto esd = CellEsdFromCovariance(basis, p0, cov); + + const Eigen::MatrixXd chol = cov.llt().matrixL(); + std::mt19937 rng(7); + std::normal_distribution normal; + constexpr int SAMPLES = 20000; + std::array sum{}, sum2{}; + for (int s = 0; s < SAMPLES; s++) { + Eigen::VectorXd z(n); + for (int k = 0; k < n; k++) + z(k) = normal(rng); + const Eigen::VectorXd p = p0 + chol * z; + Eigen::Matrix3d g = Eigen::Matrix3d::Zero(); + for (int k = 0; k < n; k++) + g += p(k) * basis[k]; + const auto cell = CellFromReciprocalMetric(g); + for (int j = 0; j < 6; j++) { + sum[j] += cell[j]; + sum2[j] += cell[j] * cell[j]; + } + } + for (int j = 0; j < 6; j++) { + const double mean = sum[j] / SAMPLES; + const double sd = std::sqrt(std::max(0.0, sum2[j] / SAMPLES - mean * mean)); + if (esd[j] < CELL_ESD_FIXED) + CHECK(sd < 1e-7); + else + CHECK(sd == Catch::Approx(esd[j]).epsilon(0.03)); + } + } +} + +TEST_CASE("CellEsd: the metric basis has the free parameters of each crystal system", "[cell_esd]") { + CHECK(MetricBasis(Sg("P 1")).size() == 6); + CHECK(MetricBasis(Sg("P 1 21 1")).size() == 4); + CHECK(MetricBasis(Sg("C 1 2/c 1")).size() == 4); + CHECK(MetricBasis(Sg("P 21 21 21")).size() == 3); + CHECK(MetricBasis(Sg("I 41/a")).size() == 2); + CHECK(MetricBasis(Sg("P 61 2 2")).size() == 2); + CHECK(MetricBasis(Sg("R 3 :H")).size() == 2); + CHECK(MetricBasis(Sg("R 3 :R")).size() == 2); + CHECK(MetricBasis(Sg("I a -3 d")).size() == 1); +} + +TEST_CASE("CellEsd: the metric basis spans the reciprocal metric of a cell of the group", "[cell_esd]") { + // R G* R^T, not R^T G* R: only the right one keeps a hexagonal G* (g12* = +a*^2/2) in the span. + double residual = 1.0; + Coordinates(MetricBasis(Sg("P 61 2 2")), ReciprocalMetric(5.4, 5.4, 56.3, 90, 90, 120), residual); + CHECK(residual < 1e-12); + Coordinates(MetricBasis(Sg("R 3 :R")), ReciprocalMetric(7.0, 7.0, 7.0, 80, 80, 80), residual); + CHECK(residual < 1e-12); + Coordinates(MetricBasis(Sg("P 1 21 1")), ReciprocalMetric(11.0, 6.5, 11.3, 90, 96, 90), residual); + CHECK(residual < 1e-12); + // and nothing outside it + Coordinates(MetricBasis(Sg("P 1 21 1")), ReciprocalMetric(11.0, 6.5, 11.3, 91, 96, 90), residual); + CHECK(residual > 1e-4); +} + +TEST_CASE("CellEsd: e.s.d.s from a covariance follow the constraints and match sampling", "[cell_esd]") { + CheckAgainstSampling("P 1", 7.1, 8.3, 9.2, 81.0, 76.0, 69.0); + CheckAgainstSampling("P 1 21 1", 11.26, 6.54, 11.26, 90.0, 95.9, 90.0); + CheckAgainstSampling("P 21 21 21", 8.35, 9.58, 27.07, 90.0, 90.0, 90.0); + CheckAgainstSampling("P 61 2 2", 5.42, 5.42, 56.28, 90.0, 90.0, 120.0); + CheckAgainstSampling("R 3 :R", 7.0, 7.0, 7.0, 80.0, 80.0, 80.0); + + // Tetragonal: a and b share one e.s.d., the angles have none. + const auto basis = MetricBasis(Sg("I -4 2 d")); + double residual = 0.0; + const Eigen::VectorXd p = Coordinates(basis, ReciprocalMetric(7.4, 7.4, 6.9, 90, 90, 90), residual); + const auto esd = CellEsdFromCovariance(basis, p, Eigen::MatrixXd::Identity(2, 2) * 1e-12); + CHECK(esd[0] > 0.0); + CHECK(esd[0] == Catch::Approx(esd[1]).epsilon(1e-9)); + for (int j = 3; j < 6; j++) + CHECK(esd[j] < CELL_ESD_FIXED); +} + +TEST_CASE("CellEsd: the d-spacing fit recovers the cell, and its e.s.d. is its scatter", "[cell_esd]") { + const double a = 11.26, b = 6.54, c = 11.26, beta = 95.9; + const Eigen::Matrix3d gstar = ReciprocalMetric(a, b, c, 90, beta, 90); + std::vector clean; + for (int h = -8; h <= 8; h++) + for (int k = -5; k <= 5; k++) + for (int l = -8; l <= 8; l++) { + const Eigen::Vector3d v(h, k, l); + const double q = std::sqrt(v.dot(gstar * v)); + if (q > 0.0 && q < 1.0 / 1.0) + clean.push_back(CellObservation{h, k, l, static_cast(q), + static_cast(clean.size() * 37 % 100)}); + } + REQUIRE(clean.size() > 2500); + + const auto &sg = Sg("P 1 21 1"); + std::mt19937 rng(11); + std::normal_distribution noise(0.0, 2e-4); // e.s.d. on q, 1/A + std::vector fitted_a, fitted_beta; + std::vector fits; + for (int trial = 0; trial < 40; trial++) { + std::vector obs = clean; + for (auto &o : obs) + o.q = static_cast(o.q + noise(rng)); + // One gross outlier, which the 4-sigma rejection takes out. + obs[17].q *= 1.05f; + const auto e = EstimateCellEsd(obs, sg, 0.62); + REQUIRE(e.has_value()); + // a 4-sigma tail of 3000 Gaussian residuals may lose one or two more + CHECK(e->observations <= static_cast(obs.size()) - 1); + CHECK(e->observations >= static_cast(obs.size()) - 4); + fitted_a.push_back(e->fitted.a); + fitted_beta.push_back(e->fitted.beta); + fits.push_back(*e); + } + const auto sd = [](const std::vector &v) { + double m = 0, s = 0; + for (double x : v) m += x / v.size(); + for (double x : v) s += (x - m) * (x - m) / (v.size() - 1); + return std::sqrt(s); + }; + const CellEsd &e = fits.front(); + CHECK(sd(fitted_a) == Catch::Approx(e.esd_lsq[0]).epsilon(0.3)); + CHECK(sd(fitted_beta) == Catch::Approx(e.esd_lsq[4]).epsilon(0.3)); + CHECK(std::abs(e.fitted.a - a) < 5 * e.esd_lsq[0] + 1e-4); // float cell + CHECK(e.esd_lsq[3] < CELL_ESD_FIXED); + CHECK(e.esd_lsq[5] < CELL_ESD_FIXED); + // Pure centroid noise: the blocks agree with the least squares, and the larger of the two is written. + REQUIRE(e.esd_blocks.has_value()); + CHECK(e.blocks == 8); + CHECK((*e.esd_blocks)[0] == Catch::Approx(e.esd_lsq[0]).epsilon(0.8)); + for (int j = 0; j < 6; j++) + CHECK(e.esd[j] == std::max(e.esd_lsq[j], (*e.esd_blocks)[j])); + CHECK(e.theta_max_deg == Catch::Approx(std::asin(0.62 / 2.0) * 180.0 / M_PI).margin(0.1)); +} + +TEST_CASE("CellEsd: the SHELX .ins carries CELL, ZERR, LATT and SYMM", "[cell_esd][write_reflections]") { + CellEsd esd; + esd.esd = {0.0014, 0.0013, 0.0012, 0.0, 0.0034, 0.0}; + const UnitCell cell{11.262f, 6.543f, 11.257f, 90.0f, 95.916f, 90.0f}; + const auto path = (std::filesystem::temp_directory_path() / "rugnux_cell_esd_test.ins").string(); + const auto read_ins = [&](int sg_number, const CellEsd *e) { + DiffractionExperiment x; + x.IncidentEnergy_keV(20.0f); + x.SpaceGroupNumber(sg_number); + WriteShelxIns(cell, e, x, "p", path); + std::ifstream in(path); + std::stringstream ss; + ss << in.rdbuf(); + std::filesystem::remove(path); + return ss.str(); + }; + + // P 21/c: centrosymmetric with -1 at the origin, so LATT +1 and one operator of each pair. + std::string ins = read_ins(14, &esd); + CHECK(ins.find("TITL p in P 1 21/c 1\n") == 0); + CHECK(ins.find(" 11.26200 6.54300 11.25700 90.0000 95.9160 90.0000\n") != std::string::npos); + CHECK(ins.find("\nCELL 0.6199") != std::string::npos); + CHECK(ins.find("\nZERR 4 0.00140 0.00130 0.00120 0.0000 0.0034 0.0000\n") != std::string::npos); + CHECK(ins.find("\nLATT 1\n") != std::string::npos); + CHECK(ins.find("\nSYMM -X,Y+1/2,-Z+1/2\n") != std::string::npos); + CHECK(ins.find("\nSYMM X,-Y+1/2,Z+1/2\n") == std::string::npos); + CHECK(ins.find("\nHKLF 4\nEND\n") != std::string::npos); + + // P 21 21 21: no inversion, LATT -1 and three operators. + ins = read_ins(19, nullptr); + CHECK(ins.find("\nLATT -1\n") != std::string::npos); + CHECK(ins.find("\nZERR") == std::string::npos); + size_t symm = 0; + for (size_t pos = ins.find("SYMM"); pos != std::string::npos; pos = ins.find("SYMM", pos + 1)) + symm++; + CHECK(symm == 3); + + // C 2/c: the C centring, centrosymmetric, Z = 8. + ins = read_ins(15, &esd); + CHECK(ins.find("\nLATT 7\n") != std::string::npos); + CHECK(ins.find("\nZERR 8 ") != std::string::npos); +} + +TEST_CASE("CellEsd: the mmCIF carries the e.s.d.s, '.' where the symmetry fixes the parameter", + "[cell_esd][write_reflections]") { + DiffractionExperiment x; + x.IncidentEnergy_keV(20.0f); + x.SpaceGroupNumber(4); + const UnitCell cell{11.262f, 6.543f, 11.257f, 90.0f, 95.916f, 90.0f}; + x.SetUnitCell(cell); + CellEsd esd; + esd.esd = {0.0014, 0.0013, 0.0012, 0.0, 0.0034, 0.0}; + esd.reflections = 1234; + esd.theta_min_deg = 1.5; + esd.theta_max_deg = 18.0; + std::vector refl(1); + refl[0].h = 1; refl[0].I = 10.0f; refl[0].sigma = 1.0f; refl[0].d = 5.0f; + const auto path = (std::filesystem::temp_directory_path() / "rugnux_cell_esd_test.cif").string(); + WriteMmcifReflections(refl, cell, x, MergeStatistics{}, ErrorModelReport{}, TwinningAnalysisResult{}, + OutputProvenance{}, path, 1, &esd); + gemmi::cif::Document doc = gemmi::read_cif_gz(path); + std::filesystem::remove(path); + gemmi::cif::Block &block = doc.blocks[0]; + CHECK(*block.find_value("_cell.length_a") == "11.2620"); + CHECK(*block.find_value("_cell.length_a_esd") == "0.00140"); + CHECK(*block.find_value("_cell.angle_beta_esd") == "0.0034"); + CHECK(*block.find_value("_cell.angle_alpha_esd") == "."); + CHECK(*block.find_value("_cell_measurement.reflns_used") == "1234"); + CHECK(*block.find_value("_cell_measurement.theta_max") == "18.00"); +}