diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index d2ee03df..b6fe74f3 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,5 +1,10 @@ # Changelog ## 1.0.0 +### 1.0.0-rc.164 +This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. + +* rugnux: `--export-unmerged ` writes the integrated observations as an unmerged MTZ - readable by aimless, pointless, careless and `iotbx.merging_statistics` - in `--mode mx` and `--mode scale` alike. + ### 1.0.0-rc.163 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 20012b5e..cd72a274 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -537,6 +538,168 @@ void WriteShelxHklReflections(const std::vector &reflections, out.close(); } +void WriteUnmergedMtzReflections(const std::vector &outcomes, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + const std::string &filename) { + gemmi::Mtz mtz; + mtz.spacegroup = gemmi::find_spacegroup_by_number(experiment.GetSpaceGroupNumber().value_or(1)); + mtz.set_cell_for_all(unitCell); + mtz.title = "Rugnux unmerged observations"; + mtz.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction"); + mtz.add_base(); // the HKL_base dataset and the H K L columns + gemmi::Mtz::Dataset &ds = mtz.add_dataset("unmerged"); + ds.crystal_name = experiment.GetSampleName(); + ds.wavelength = experiment.GetWavelength_A(); + + // 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; 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. + mtz.add_column("M/ISYM", 'Y', 0, -1, false); + mtz.add_column("BATCH", 'B', 0, -1, false); + mtz.add_column("I", 'J', 0, -1, false); + mtz.add_column("SIGI", 'Q', 0, -1, false); + mtz.add_column("FRACTIONCALC", 'R', 0, -1, false); + mtz.add_column("XDET", 'R', 0, -1, false); + mtz.add_column("YDET", 'R', 0, -1, false); + mtz.add_column("ROT", 'R', 0, -1, false); + mtz.add_column("LP", 'R', 0, -1, false); + mtz.add_column("FLAG", 'I', 0, -1, false); + mtz.add_column("DELPHI", 'R', 0, -1, false); + mtz.add_column("ZETA", 'R', 0, -1, false); + mtz.add_column("BGMEAN", 'R', 0, -1, false); + mtz.add_column("BGVAR", 'R', 0, -1, false); + + const auto gon = experiment.GetGoniometer(); + // A scanning goniometer means rugnux integrated a rocking curve image by image, so every + // observation here is one PART of a reflection; on a still each is a whole measurement. The + // difference 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 + // neither, they take each partial for a whole reflection and neither the symmetry determination + // nor the scaling survives it. + const bool partials = gon && gon->IsScanning(); + 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 the Lorentz-polarization factor applied and + // nothing else, which is what IOBS means in every unmerged format (LP records the factor, so the + // raw counts are I/LP). LP is geometry, not a scale, and a program that reads this file has no + // way to recover it. 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. + gemmi::UnmergedHklMover hkl_mover(mtz.spacegroup); + std::set batch_numbers; + for (const auto &outcome : outcomes) { + for (const auto &r : outcome.reflections) { + std::array hkl{r.h, r.k, r.l}; + const int isym = hkl_mover.move_to_asu(hkl); + const int batch = 1 + static_cast(std::lround(r.image_number)); + batch_numbers.insert(batch); + mtz.data.push_back(static_cast(hkl[0])); + mtz.data.push_back(static_cast(hkl[1])); + mtz.data.push_back(static_cast(hkl[2])); + mtz.data.push_back(static_cast((partials ? 256 : 0) + isym)); + mtz.data.push_back(static_cast(batch)); + mtz.data.push_back(r.I * r.rlp); + mtz.data.push_back(r.sigma * r.rlp); + mtz.data.push_back(r.partiality); + // A reflection too weak to have a measured centroid still has a predicted position, and a + // detector coordinate is what the scale models downstream read this column for. + mtz.data.push_back(std::isfinite(r.observed_x) ? r.observed_x : r.predicted_x); + mtz.data.push_back(std::isfinite(r.observed_y) ? r.observed_y : r.predicted_y); + mtz.data.push_back(phi_start_deg(r.image_number) + wedge_deg / 2.0f); + mtz.data.push_back(r.rlp); + mtz.data.push_back(0.0f); // FLAG: nothing here is a rejected observation + mtz.data.push_back(r.delta_phi_deg); + mtz.data.push_back(r.zeta); + mtz.data.push_back(r.bkg); + mtz.data.push_back(r.var_bkg); + } + } + mtz.nreflections = static_cast(mtz.data.size() / mtz.columns.size()); + + // 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 lattice_at_zero; + std::optional mosaicity_deg; + for (const auto &outcome : outcomes) { + if (outcome.reflections.empty() || outcome.latt.CalcVolume() <= 1.0f) + 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; + mosaicity_deg = outcome.mosaicity_deg; + 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(ds.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[21] = mosaicity_deg.value_or(0.0f); // crydat(0), the reflecting range + 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(experiment.GetXPixelsNum()); + batch.floats[115] = 1.0f; + batch.floats[116] = static_cast(experiment.GetYPixelsNum()); + + for (const int number : batch_numbers) { + batch.number = number; + batch.floats[36] = phi_start_deg(static_cast(number - 1)); // phistt + batch.floats[37] = batch.floats[36] + wedge_deg; // phiend + mtz.batches.push_back(batch); + } + + mtz.sort(5); // by H K L M/ISYM BATCH, the order POINTLESS leaves an unmerged file in + mtz.write_to_file(filename); +} + void WriteReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index fb3f86c5..997b4d97 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -9,6 +9,7 @@ #include "../common/Reflection.h" #include "../common/UnitCell.h" #include "../common/DiffractionExperiment.h" +#include "IntegrationOutcome.h" struct MergeStatistics; struct TwinningAnalysisResult; @@ -46,6 +47,19 @@ void WriteShelxHklReflections(const std::vector &reflections, const std::string &filename, size_t nthreads); +// Unmerged observations, one row per integrated reflection, 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 frame) and says whether the observation is a partial. +// One observation per integrated box, so on a rotation run a reflection arrives as a run of partials +// over consecutive batches, for the reader to sum. The intensities carry the Lorentz-polarization +// factor and nothing else - the partiality and the per-image scale are left for the reading program +// to fit, since every program this file is for fits a scale model of its own. +void WriteUnmergedMtzReflections(const std::vector &outcomes, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + const std::string &filename); + void WriteReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index fff06f75..00d827c1 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2986,6 +2986,18 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } } + // Unmerged observations (--export-unmerged). After the merge, because that is where the space + // group is settled and the file has to be written in it, but from the integrated observations + // rather than the merged ones: the partiality and the per-image scale are left for the reading + // program, which fits a scale model of its own. + if (full && !cancelled_ && write_output && !geometry_prepass + && !config_.unmerged_export_path.empty() && result.consensus_cell) { + if (observer) observer->OnPhase("Writing unmerged reflections"); + WriteUnmergedMtzReflections(indexer->GetIntegrationOutcome(), *result.consensus_cell, + experiment_, config_.unmerged_export_path); + logger.Info("Unmerged observations written to {}", config_.unmerged_export_path); + } + // When merging is on, the _process.h5 is skipped because the merged reflections are the wanted // output and the per-image file is large. If nothing indexed there are no merged reflections, so // the run would finish successfully having written no file at all. Write one now: the per-image diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 2e6c40fe..23343f65 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -119,6 +119,10 @@ struct ProcessConfig { // Diagnostic: if set, the -P rot3d combine writes the unmerged fulls here (for comparison vs XDS). std::string observation_dump_path; + // If set, the integrated observations are written here as an unmerged MTZ (--export-unmerged), + // for aimless / pointless / careless and anything else that scales unmerged data itself. + std::string unmerged_export_path; + // If set, after merging, validate the merged intensities against this atomic model (PDB): // compute R-free with an optimized bulk solvent and write 2Fo-Fc / Fo-Fc maps (--model). std::string model_path; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 8217882e..ee8bc783 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -162,6 +162,7 @@ void print_usage() { std::cout << " -z, --reference-mtz Reference MTZ file" << std::endl; std::cout << " --reference-column