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

jfjoch_process: Save reflections as mmCIF (experimental)

Browse Source
This commit is contained in:
leonarski_f
2026-02-12 18:09:35 +01:00
parent 1dedd3de33
commit 38199c480f
4 changed files with 252 additions and 13 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/CMakeLists.txt
+3 -1
View File
@@ -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
image_analysis/WriteMmcif.cpp
+163
View File
@@ -0,0 +1,163 @@
// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute
// SPDX-License-Identifier: GPL-3.0-only
#include "WriteMmcif.h"
#include <cmath>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <stdexcept>
#include <ctime>
#include <chrono>
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<FrenchWilsonReflection>& 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<FrenchWilsonReflection>& 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);
}
image_analysis/WriteMmcif.h
+48
View File
@@ -0,0 +1,48 @@
// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
#include <ostream>
#include <string>
#include <vector>
#include <optional>
#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_<name>
std::string detector_name; // e.g. "JUNGFRAU 4M"
std::optional<float> wavelength_A; // incident wavelength
std::optional<float> detector_distance_mm;
std::optional<float> sample_temperature_K;
std::optional<std::string> sample_name;
std::optional<std::string> software_version; // jfjoch version string
std::optional<std::string> source;
std::optional<std::string> 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<FrenchWilsonReflection>& 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<FrenchWilsonReflection>& reflections,
const MmcifMetadata& meta);
tools/jfjoch_process.cpp
+38 -12
View File
@@ -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 {<options>} <input.h5>");
@@ -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 {