diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 819a2b11..2cf690d1 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -32,7 +32,9 @@ ADD_LIBRARY(JFJochImageAnalysis STATIC RotationParameters.cpp RotationParameters.h RotationSpotAccumulator.cpp - RotationSpotAccumulator.h) + RotationSpotAccumulator.h + WriteMmcif.cpp + WriteMmcif.h) FIND_PACKAGE(Eigen3 3.4 REQUIRED NO_MODULE) # provides Eigen3::Eigen diff --git a/image_analysis/WriteMmcif.cpp b/image_analysis/WriteMmcif.cpp new file mode 100644 index 00000000..28b4d994 --- /dev/null +++ b/image_analysis/WriteMmcif.cpp @@ -0,0 +1,163 @@ +// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "WriteMmcif.h" + +#include +#include +#include +#include +#include +#include +#include + +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; +} + +/// Format a double with given decimal places; returns "?" for non-finite. +std::string Fmt(double val, int decimals = 4) { + if (!std::isfinite(val)) + return "?"; + std::ostringstream ss; + ss << std::fixed << std::setprecision(decimals) << val; + return ss.str(); +} + +/// 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; +} + +} // namespace + +void WriteMmcifReflections(std::ostream& out, + const std::vector& reflections, + const MmcifMetadata& meta) { + + out << std::fixed; + + // ---------- data block ---------- + out << "data_" << meta.data_block_name << "\n"; + out << "#\n"; + + // ---------- _audit ---------- + out << "_audit.revision_id 1\n"; + out << "_audit.creation_date " << CurrentDateISO() << "\n"; + out << "_audit.update_record 'Initial release'\n"; + out << "#\n"; + + // ---------- _software ---------- + out << "_software.name 'Jungfraujoch'\n"; + if (meta.software_version.has_value()) + out << "_software.version " << CifStr(meta.software_version.value()) << "\n"; + out << "_software.classification reduction\n"; + out << "#\n"; + + // ---------- _cell ---------- + out << "_cell.length_a " << Fmt(meta.unit_cell.a, 3) << "\n"; + out << "_cell.length_b " << Fmt(meta.unit_cell.b, 3) << "\n"; + out << "_cell.length_c " << Fmt(meta.unit_cell.c, 3) << "\n"; + out << "_cell.angle_alpha " << Fmt(meta.unit_cell.alpha, 2) << "\n"; + out << "_cell.angle_beta " << Fmt(meta.unit_cell.beta, 2) << "\n"; + out << "_cell.angle_gamma " << Fmt(meta.unit_cell.gamma, 2) << "\n"; + + // ---------- _symmetry ---------- + out << "_symmetry.space_group_name_H-M " << CifStr(meta.space_group_name) << "\n"; + out << "_symmetry.Int_Tables_number " << meta.space_group_number << "\n"; + out << "#\n"; + + // ---------- _diffrn_source / _diffrn_detector ---------- + if (meta.source.has_value()) { + out << "_diffrn_source.pdbx_synchrotron_site" + << CifStr(meta.source.value()) << "\n"; + } + + if (meta.beamline.has_value()) { + out << "_diffrn_source.pdbx_synchrotron_beamline" + << CifStr(meta.beamline.value()) << "\n"; + } + + if (meta.wavelength_A.has_value()) { + out << "_diffrn_radiation_wavelength.wavelength " + << Fmt(meta.wavelength_A.value(), 5) << "\n"; + } + + if (!meta.detector_name.empty()) { + out << "_diffrn_detector.detector " << CifStr(meta.detector_name) << "\n"; + } + + if (meta.detector_distance_mm.has_value()) { + out << "_diffrn_detector.distance " + << Fmt(meta.detector_distance_mm.value(), 2) << "\n"; + } + + if (meta.sample_temperature_K.has_value()) { + out << "_diffrn.ambient_temp " + << Fmt(meta.sample_temperature_K.value(), 1) << "\n"; + } + + if (meta.sample_name.has_value() && !meta.sample_name->empty()) { + out << "_entity.id 1\n"; + out << "_entity.type polymer\n"; + out << "_entity.pdbx_description " << CifStr(meta.sample_name.value()) << "\n"; + } + + // ---------- _refln loop ---------- + out << "loop_\n"; + out << "_refln.index_h\n"; + out << "_refln.index_k\n"; + out << "_refln.index_l\n"; + out << "_refln.F_meas_au\n"; + out << "_refln.F_meas_sigma_au\n"; + out << "_refln.intensity_meas\n"; + out << "_refln.intensity_sigma\n"; + out << "_refln.status\n"; + + for (const auto& r : reflections) { + out << std::setw(5) << r.h << " " + << std::setw(5) << r.k << " " + << std::setw(5) << r.l << " " + << std::setw(12) << Fmt(r.F, 4) << " " + << std::setw(12) << Fmt(r.sigmaF, 4) << " " + << std::setw(14) << Fmt(r.I, 4) << " " + << std::setw(14) << Fmt(r.sigmaI, 4) << " " + << "o" // 'o' = observed + << "\n"; + } + + out << "#\n"; + out << "# End of reflections\n"; +} + +void WriteMmcifReflections(const std::string& path, + const std::vector& reflections, + const MmcifMetadata& meta) { + std::ofstream file(path); + if (!file) + throw std::runtime_error("WriteMmcifReflections: cannot open " + path); + WriteMmcifReflections(file, reflections, meta); + file.close(); + if (!file) + throw std::runtime_error("WriteMmcifReflections: I/O error writing " + path); +} \ No newline at end of file diff --git a/image_analysis/WriteMmcif.h b/image_analysis/WriteMmcif.h new file mode 100644 index 00000000..10e05b22 --- /dev/null +++ b/image_analysis/WriteMmcif.h @@ -0,0 +1,48 @@ +// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include + +#include "scale_merge/FrenchWilson.h" +#include "../common/UnitCell.h" +#include "../symmetry/gemmi/symmetry.hpp" + +/// Metadata needed to write a meaningful mmCIF reflection file. +struct MmcifMetadata { + // Required + UnitCell unit_cell{}; + std::string space_group_name; // e.g. "P 21 21 21" + int space_group_number = 1; + + // Optional but recommended + std::string data_block_name = "jfjoch"; // CIF data_ + std::string detector_name; // e.g. "JUNGFRAU 4M" + std::optional wavelength_A; // incident wavelength + std::optional detector_distance_mm; + std::optional sample_temperature_K; + std::optional sample_name; + std::optional software_version; // jfjoch version string + + std::optional source; + std::optional beamline; +}; + +/// Write a PDBx/mmCIF reflection file (containing _refln loop with F/sigmaF/I/sigmaI) +/// to the given output stream. +/// +/// The file follows the conventions expected by CCP4/CCTBX/Phenix for +/// structure-factor mmCIF files. +void WriteMmcifReflections(std::ostream& out, + const std::vector& reflections, + const MmcifMetadata& meta); + +/// Convenience overload that writes to a file path. +/// Throws std::runtime_error on I/O failure. +void WriteMmcifReflections(const std::string& path, + const std::vector& reflections, + const MmcifMetadata& meta); \ No newline at end of file diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 41fb0de6..9510388d 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -25,6 +25,7 @@ #include "../receiver/JFJochReceiverPlots.h" #include "../compression/JFJochCompressor.h" #include "../image_analysis/scale_merge/FrenchWilson.h" +#include "../image_analysis/WriteMmcif.h" // (actually include at top of file) void print_usage(Logger &logger) { logger.Info("Usage ./jfjoch_analysis {} "); @@ -489,19 +490,44 @@ int main(int argc, char **argv) { auto fw = FrenchWilson(scale_result->merged, fw_opts); - const std::string fw_path = output_prefix + "_amplitudes.hkl"; - std::ofstream fw_file(fw_path); - if (!fw_file) { - logger.Error("Cannot open {} for writing", fw_path); - } else { - fw_file << "# h k l F sigmaF I_fw sigmaI\n"; - for (const auto& r : fw) { - fw_file << r.h << " " << r.k << " " << r.l << " " - << r.F << " " << r.sigmaF << " " - << r.I << " " << r.sigmaI << "\n"; + { + + MmcifMetadata cif_meta; + + // Unit cell — from rotation indexing result or experiment setting + if (rotation_indexer_ret.has_value()) { + cif_meta.unit_cell = rotation_indexer_ret->lattice.GetUnitCell(); + } else if (experiment.GetUnitCell().has_value()) { + cif_meta.unit_cell = experiment.GetUnitCell().value(); + } + + // Space group + if (auto sg = experiment.GetGemmiSpaceGroup(); sg.has_value()) { + cif_meta.space_group_name = sg->hm; + cif_meta.space_group_number = sg->number; + } else if (space_group) { + cif_meta.space_group_name = space_group->hm; + cif_meta.space_group_number = space_group->number; + } + + // Detector & experiment info + cif_meta.detector_name = experiment.GetDetectorDescription(); + cif_meta.wavelength_A = experiment.GetWavelength_A(); + cif_meta.detector_distance_mm = experiment.GetDetectorDistance_mm(); + cif_meta.sample_temperature_K = experiment.GetSampleTemperature_K(); + cif_meta.sample_name = experiment.GetSampleName(); + cif_meta.data_block_name = output_prefix; + + cif_meta.beamline = experiment.GetInstrumentName(); + cif_meta.source = experiment.GetSourceName(); + + const std::string cif_path = output_prefix + "_amplitudes.cif"; + try { + WriteMmcifReflections(cif_path, fw, cif_meta); + logger.Info("Wrote mmCIF reflections to {}", cif_path); + } catch (const std::exception& e) { + logger.Error("Failed to write mmCIF: {}", e.what()); } - fw_file.close(); - logger.Info("French-Wilson: wrote {} amplitudes to {}", fw.size(), fw_path); } } } else {