From 2e296e6bed37b21d76b1c8a0d5d0abaa182520c2 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 26 Oct 2026 15:48:40 +0100 Subject: [PATCH] jfjoch_writer: the live refinement's model, maps and blobs are written beside the master file In broker live mode the refined model, its map coefficients and the difference-map blobs stayed in the broker for the viewer routes and were lost at the next run. They now go to the writer in the END message as the placed model and maps already do - EndMessage::refined_model (PlacedModelMessage), refined_map_coefficients (MtzTable) and blobs (DensityBlobsMessage, one vector per field) - and the writer writes _refined_model.cif (+ .pdb), _refined_maps.mtz and _blobs.json, the names the CLI gives them. The blobs' JSON is made by one function (DensityBlobsJson) that the CLI now writes through as well, so the two files are byte-identical. Optional CBOR fields, no version bump. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit c65220d96a96d2b1897f3c6632c94f9379b670a3) --- common/JFJochMessages.h | 17 ++++ docs/CBOR.md | 5 +- docs/JFJOCH_BROKER.md | 8 +- frame_serialize/CBORStream2Deserializer.cpp | 30 +++++++ frame_serialize/CBORStream2Serializer.cpp | 22 +++++ .../structure_refinement/DensityBlobs.cpp | 19 +++++ .../structure_refinement/DensityBlobs.h | 5 ++ receiver/JFJochReceiver.cpp | 3 + receiver/JFJochReceiverRugnux.cpp | 16 +++- rugnux/RugnuxOutputFiles.cpp | 16 +--- tests/RugnuxResultFilesTest.cpp | 85 +++++++++++++++++++ writer/WriteRugnuxFiles.cpp | 80 ++++++++++++++--- writer/WriteRugnuxFiles.h | 12 ++- 13 files changed, 282 insertions(+), 36 deletions(-) diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index 70f762546..390707bd1 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -406,6 +406,19 @@ struct PlacedModelMessage { std::string metadata_mmcif; }; +// The unexplained blobs of the refined model's mFo-DFc map (DensityBlob, image_analysis/structure_refinement/ +// DensityBlobs.h), strongest first, one vector per field, for the writer to write _blobs.json. kind is +// DensityBlobKindName's; centroid x, y, z in A (orthogonal), volume in A^3, score in sigma A^3. +struct DensityBlobsMessage { + std::vector kind; + std::vector x, y, z; + std::vector peak_sigma, volume, score; + std::vector nearest_atom; + std::vector nearest_atom_distance; + std::vector polar_contacts, metal_contacts; + std::vector anomalous_sigma; +}; + struct EndMessage { // The sample transformation chain, base first (see DetectorTransformation). Optional: when it is // empty the writer builds the same chain itself from the start message, so a producer that does @@ -531,6 +544,10 @@ struct EndMessage { std::optional merged; // MergedMtzTable (image_analysis/WriteReflections.h) std::optional map_coefficients; // ModelMaps::coefficients (ModelValidation.h) std::optional placed_model; + // And of its refinement, where one ran: _refined_model.cif, _refined_maps.mtz and _blobs.json. + std::optional refined_model; + std::optional refined_map_coefficients; + std::optional blobs; }; struct MetadataMessage { diff --git a/docs/CBOR.md b/docs/CBOR.md index 20f6ff90a..ccaaef192 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -348,10 +348,13 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | merged | object (optional) | Merged reflections of an in-broker analysis run, written by the writer as `.mtz`. An MTZ table (layout below) | | | map_coefficients | object (optional) | Map coefficients of an in-broker analysis run, written as `_maps.mtz`. An MTZ table (layout below) | | | placed_model | object (optional) | Atomic model placed against the data, written as `_model.cif` and `_model.pdb`. `name` (string), `unit_cell` (object, as above), `space_group_name` (the data's, gemmi extended H-M), `model_space_group_name` (the coordinates'), then one array per atom: `model_number`, `seq_num`, `charge` (int32), `chain`, `residue_name`, `atom_name`, `element` (Array(string)), `icode` (uint8, space where none), `het_flag` (uint8, `A`, `H` or 0 where unknown), `altloc` (uint8, 0 where none), `x`, `y`, `z` (double, Angstrom), `b_iso`, `occupancy` (float), `aniso` (float, 6 per atom: U11 U22 U33 U12 U13 U23; empty where no atom has one), `subchain`, `entity_id` (Array(string), mmCIF label_asym_id / label_entity_id), `label_seq` (int32, label_seq_id, INT32_MIN where none) - these three optional; `metadata_mmcif` (string, optional): the input model's metadata (entities, sequences, references, authors, software, secondary structure) as mmCIF categories without coordinates, parsed by the writer with gemmi | | +| refined_model | object (optional) | The placed model after refinement, written as `_refined_model.cif` and `_refined_model.pdb`; same layout as `placed_model` | | +| refined_map_coefficients | object (optional) | Map coefficients of the refined model, written as `_refined_maps.mtz`. An MTZ table (layout below) | | +| blobs | object (optional) | Unexplained blobs of the refined model's mFo-DFc map, strongest first, written as `_blobs.json`. One array per field, one entry per blob: `kind` (Array(string): `WATER`, `ION`, `LIGAND`, `MODEL_ERROR`, `UNASSIGNED`), `x`, `y`, `z` (double, centroid, Angstrom), `peak_sigma` (double, map rms), `volume` (double, Angstrom^3), `score` (double, sigma Angstrom^3), `nearest_atom` (Array(string)), `nearest_atom_distance` (double, Angstrom), `polar_contacts`, `metal_contacts` (int32), `anomalous_sigma` (double, 0 where there is no anomalous map) | | End-message vector fields are optional. When present, they provide master-file summary data so readers can inspect scan-level and per-image analysis results without opening every linked data file. Missing optional per-image values are encoded by the producer as zero unless otherwise noted. -An MTZ table (`merged`, `map_coefficients`) carries an MTZ file's content as the writer writes it, column by column: +An MTZ table (`merged`, `map_coefficients`, `refined_map_coefficients`) carries an MTZ file's content as the writer writes it, column by column: | Field name | Type | Description | |------------------|---------------|------------------------------------------------------------------------------------------------------| diff --git a/docs/JFJOCH_BROKER.md b/docs/JFJOCH_BROKER.md index 10d2ee14b..9e48cedb8 100644 --- a/docs/JFJOCH_BROKER.md +++ b/docs/JFJOCH_BROKER.md @@ -182,9 +182,11 @@ collection - one slot per image. Where it holds fewer, the collection is only fo as with the mode off; the live status reports it as skipped and the collection's statistics carry a warning saying how many images the buffer holds against how many the collection has. -The results - the merged data, the map coefficients and, with a model, the placed model - travel to -the writer in the end message, which writes them beside the master file as `.mtz`, -`_maps.mtz`, `_model.cif` and `_model.pdb`. An end message +The results - the merged data, the map coefficients and, with a model, the placed model and its +refinement - travel to the writer in the end message, which writes them beside the master file as +`.mtz`, `_maps.mtz`, `_model.cif` and `_model.pdb`, +and `_refined_model.cif` (and `.pdb`), `_refined_maps.mtz` and +`_blobs.json`, as rugnux names them. An end message that would be too large for the pushers' buffer (256 MiB) goes without them, with an error in the log and a warning in the live status; the collection itself is not affected. diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index 066bad575..f5e56cc52 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -1679,6 +1679,30 @@ namespace { return m; } + DensityBlobsMessage ProcessDensityBlobsElement(CborValue &value) { + DensityBlobsMessage b; + CborValue map_value; + cborErr(cbor_value_enter_container(&value, &map_value)); + while (!cbor_value_at_end(&map_value)) { + auto key = GetCBORString(map_value); + if (key == "kind") GetCBORStringArray(map_value, b.kind); + else if (key == "x") GetCBORDoubleArray(map_value, b.x); + else if (key == "y") GetCBORDoubleArray(map_value, b.y); + else if (key == "z") GetCBORDoubleArray(map_value, b.z); + else if (key == "peak_sigma") GetCBORDoubleArray(map_value, b.peak_sigma); + else if (key == "volume") GetCBORDoubleArray(map_value, b.volume); + else if (key == "score") GetCBORDoubleArray(map_value, b.score); + else if (key == "nearest_atom") GetCBORStringArray(map_value, b.nearest_atom); + else if (key == "nearest_atom_distance") GetCBORDoubleArray(map_value, b.nearest_atom_distance); + else if (key == "polar_contacts") GetCBORInt32Array(map_value, b.polar_contacts); + else if (key == "metal_contacts") GetCBORInt32Array(map_value, b.metal_contacts); + else if (key == "anomalous_sigma") GetCBORDoubleArray(map_value, b.anomalous_sigma); + else SkipCBORValue(map_value); + } + cborErr(cbor_value_leave_container(&value, &map_value)); + return b; + } + bool ProcessEndMessageElement(EndMessage &message, CborValue &value) { if (cbor_value_at_end(&value)) return false; @@ -1783,6 +1807,12 @@ namespace { message.map_coefficients = ProcessMtzTableElement(value); else if (key == "placed_model") message.placed_model = ProcessPlacedModelElement(value); + else if (key == "refined_model") + message.refined_model = ProcessPlacedModelElement(value); + else if (key == "refined_map_coefficients") + message.refined_map_coefficients = ProcessMtzTableElement(value); + else if (key == "blobs") + message.blobs = ProcessDensityBlobsElement(value); else if (key == "rotation_lattice") { std::vector tmp; GetCBORFloatArray(value, tmp); diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 46e7510dc..ec07c26f3 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -630,6 +630,25 @@ inline void CBOR_ENC(CborEncoder &encoder, const char* key, const PlacedModelMes cborErr(cbor_encoder_close_container(&encoder, &mapEncoder)); } +inline void CBOR_ENC(CborEncoder &encoder, const char* key, const DensityBlobsMessage &b) { + CborEncoder mapEncoder; + cborErr(cbor_encode_text_stringz(&encoder, key)); + cborErr(cbor_encoder_create_map(&encoder, &mapEncoder, CborIndefiniteLength)); + CBOR_ENC(mapEncoder, "kind", b.kind); + CBOR_ENC(mapEncoder, "x", b.x); + CBOR_ENC(mapEncoder, "y", b.y); + CBOR_ENC(mapEncoder, "z", b.z); + CBOR_ENC(mapEncoder, "peak_sigma", b.peak_sigma); + CBOR_ENC(mapEncoder, "volume", b.volume); + CBOR_ENC(mapEncoder, "score", b.score); + CBOR_ENC(mapEncoder, "nearest_atom", b.nearest_atom); + CBOR_ENC(mapEncoder, "nearest_atom_distance", b.nearest_atom_distance); + CBOR_ENC(mapEncoder, "polar_contacts", b.polar_contacts); + CBOR_ENC(mapEncoder, "metal_contacts", b.metal_contacts); + CBOR_ENC(mapEncoder, "anomalous_sigma", b.anomalous_sigma); + cborErr(cbor_encoder_close_container(&encoder, &mapEncoder)); +} + // Add encoder for LatticeMessage inline void CBOR_ENC(CborEncoder &encoder, const char* key, const LatticeMessage &val) { CborEncoder mapEncoder; @@ -977,6 +996,9 @@ void CBORStream2Serializer::SerializeSequenceEnd(const EndMessage& message) { CBOR_ENC(mapEncoder, "merged", message.merged); CBOR_ENC(mapEncoder, "map_coefficients", message.map_coefficients); CBOR_ENC(mapEncoder, "placed_model", message.placed_model); + CBOR_ENC(mapEncoder, "refined_model", message.refined_model); + CBOR_ENC(mapEncoder, "refined_map_coefficients", message.refined_map_coefficients); + CBOR_ENC(mapEncoder, "blobs", message.blobs); cborErr(cbor_encoder_close_container(&encoder, &mapEncoder)); diff --git a/image_analysis/structure_refinement/DensityBlobs.cpp b/image_analysis/structure_refinement/DensityBlobs.cpp index 67a76d84e..b2a9b7fb9 100644 --- a/image_analysis/structure_refinement/DensityBlobs.cpp +++ b/image_analysis/structure_refinement/DensityBlobs.cpp @@ -42,6 +42,25 @@ const char *DensityBlobKindName(DensityBlob::Kind kind) { return "UNASSIGNED"; } +DensityBlobsMessage DensityBlobsToMessage(const std::vector &blobs) { + DensityBlobsMessage m; + for (const DensityBlob &b : blobs) { + m.kind.emplace_back(DensityBlobKindName(b.kind)); + m.x.push_back(b.centroid.x); + m.y.push_back(b.centroid.y); + m.z.push_back(b.centroid.z); + m.peak_sigma.push_back(b.peak_sigma); + m.volume.push_back(b.volume); + m.score.push_back(b.score); + m.nearest_atom.push_back(b.nearest_atom); + m.nearest_atom_distance.push_back(b.nearest_atom_distance); + m.polar_contacts.push_back(b.polar_contacts); + m.metal_contacts.push_back(b.metal_contacts); + m.anomalous_sigma.push_back(b.anomalous_sigma); + } + return m; +} + // Following Coot's unmodelled blobs, Emsley et al. (2010) Acta Cryst. D66, 486-501 std::vector FindDensityBlobs(const gemmi::Grid &difference, double rms, const gemmi::Structure &model, double d_min, const gemmi::Grid *anomalous, double anomalous_rms) { diff --git a/image_analysis/structure_refinement/DensityBlobs.h b/image_analysis/structure_refinement/DensityBlobs.h index 657c22f51..499aa37af 100644 --- a/image_analysis/structure_refinement/DensityBlobs.h +++ b/image_analysis/structure_refinement/DensityBlobs.h @@ -9,6 +9,8 @@ #include "gemmi/grid.hpp" #include "gemmi/model.hpp" +#include "../../common/JFJochMessages.h" + // What a model leaves unexplained: the connected regions of positive mFo-DFc density above a cutoff ("blobs"), // each placed beside the model, measured and given a likely reading - a water, an ion, something ligand-sized, // or an error of the model itself (a misplaced or missing atom, an alternative conformation). @@ -32,6 +34,9 @@ struct DensityBlob { const char *DensityBlobKindName(DensityBlob::Kind kind); +// The blobs as the end message carries them to the writer, for _blobs.json. +DensityBlobsMessage DensityBlobsToMessage(const std::vector &blobs); + // The blobs of `difference` (an mFo-DFc map of the whole cell, rms `rms`) above 3 sigma that peak above 4 sigma, // strongest (by score) first. `model`: what the map was made from, in the map's cell. `anomalous`: the anomalous // difference map on any grid of the same cell, or null. Following the flood fill of Coot's unmodelled blobs as diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp index 516f48900..8ecc77425 100644 --- a/receiver/JFJochReceiver.cpp +++ b/receiver/JFJochReceiver.cpp @@ -210,6 +210,9 @@ void JFJochReceiver::SendEndMessage() { message.merged.reset(); message.map_coefficients.reset(); message.placed_model.reset(); + message.refined_model.reset(); + message.refined_map_coefficients.reset(); + message.blobs.reset(); if (zmq_metadata_socket != nullptr) zmq_metadata_socket->EndDataCollection(message); diff --git a/receiver/JFJochReceiverRugnux.cpp b/receiver/JFJochReceiverRugnux.cpp index 12948715f..78c8a1de4 100644 --- a/receiver/JFJochReceiverRugnux.cpp +++ b/receiver/JFJochReceiverRugnux.cpp @@ -540,6 +540,14 @@ void JFJochReceiverRugnux::AddToEndMessage(EndMessage &message) { if (o.placed_model) message.placed_model = PlacedModelFromStructure(*o.placed_model->model, o.placed_model->cell, o.placed_model->space_group.xhm()); + // And the refinement's, as the CLI writes them (_refined_model.cif, _refined_maps.mtz, _blobs.json). + if (o.refined_model) + message.refined_model = PlacedModelFromStructure(*o.refined_model->model, o.refined_model->cell, + o.refined_model->space_group.xhm()); + if (o.refined_maps) + message.refined_map_coefficients = std::move(o.refined_maps->coefficients); + if (o.blobs) + message.blobs = DensityBlobsToMessage(*o.blobs); // Every pusher serialises END into a buffer of MESSAGE_SIZE_FOR_START_END, and an END that does // not fit there fails the collection. So it is serialised here first, into a buffer of the same @@ -550,7 +558,7 @@ void JFJochReceiverRugnux::AddToEndMessage(EndMessage &message) { try { serializer.SerializeSequenceEnd(message); } catch (const JFJochException &e) { - const std::string text = fmt::format("the merged data, map coefficients and placed model do not " + const std::string text = fmt::format("the merged data, map coefficients and models do not " "fit in the end message ({} MiB) and are not written", MESSAGE_SIZE_FOR_START_END / (1024 * 1024)); logger.Error("rugnux live: {} - {}", text, e.what()); @@ -558,12 +566,18 @@ void JFJochReceiverRugnux::AddToEndMessage(EndMessage &message) { message.merged.reset(); message.map_coefficients.reset(); message.placed_model.reset(); + message.refined_model.reset(); + message.refined_map_coefficients.reset(); + message.blobs.reset(); } } catch (const std::exception &e) { logger.Error("rugnux live: the end message goes without the analysis results - {}", e.what()); message.merged.reset(); message.map_coefficients.reset(); message.placed_model.reset(); + message.refined_model.reset(); + message.refined_map_coefficients.reset(); + message.blobs.reset(); } } diff --git a/rugnux/RugnuxOutputFiles.cpp b/rugnux/RugnuxOutputFiles.cpp index 121c67270..e7d3e7a4e 100644 --- a/rugnux/RugnuxOutputFiles.cpp +++ b/rugnux/RugnuxOutputFiles.cpp @@ -6,8 +6,6 @@ #include #include -#include - #include "../writer/FileWriter.h" #include "../writer/WriteModel.h" #include "../writer/WriteRugnuxFiles.h" @@ -16,24 +14,12 @@ namespace { // The unexplained blobs of the refined model's difference map (DensityBlobs.h), strongest first. void WriteDensityBlobs(const std::vector &blobs, const std::string &path, Logger &logger) { - nlohmann::json list = nlohmann::json::array(); - for (const DensityBlob &b : blobs) - list.push_back({{"kind", DensityBlobKindName(b.kind)}, - {"centroid", {b.centroid.x, b.centroid.y, b.centroid.z}}, - {"peak_sigma", b.peak_sigma}, - {"volume_A3", b.volume}, - {"score_sigma_A3", b.score}, - {"nearest_atom", b.nearest_atom}, - {"nearest_atom_distance_A", b.nearest_atom_distance}, - {"polar_contacts", b.polar_contacts}, - {"metal_contacts", b.metal_contacts}, - {"anomalous_sigma", b.anomalous_sigma}}); std::ofstream f(path, std::ofstream::out | std::ofstream::trunc); if (!f) { logger.Warning("Cannot open {} for writing", path); return; } - f << nlohmann::json{{"blobs", list}}.dump(1) << std::endl; + f << DensityBlobsJson(DensityBlobsToMessage(blobs)); logger.Info("Difference-map blobs written to {}", path); } diff --git a/tests/RugnuxResultFilesTest.cpp b/tests/RugnuxResultFilesTest.cpp index 8792a465d..a21df3a2c 100644 --- a/tests/RugnuxResultFilesTest.cpp +++ b/tests/RugnuxResultFilesTest.cpp @@ -11,11 +11,13 @@ #include #include #include // count_atom_sites +#include #include "../common/Definitions.h" #include "../common/DiffractionExperiment.h" #include "../common/Logger.h" #include "../image_analysis/WriteReflections.h" +#include "../image_analysis/structure_refinement/DensityBlobs.h" #include "../frame_serialize/CBORStream2Serializer.h" #include "../frame_serialize/CBORStream2Deserializer.h" #include "../writer/WriteModel.h" @@ -307,3 +309,86 @@ TEST_CASE("RugnuxResultFiles: the input model's metadata reaches _model.cif thro std::filesystem::remove_all(dir); } + +TEST_CASE("RugnuxResultFiles: the refinement's model, maps and blobs reach the writer through the END block", + "[rugnux_result_files]") { + Logger logger("RugnuxResultFiles"); + const std::filesystem::path dir = std::filesystem::temp_directory_path() / "rugnux_result_files_refined"; + std::filesystem::remove_all(dir); + std::filesystem::create_directories(dir); + const std::string original = (dir / "rugnux").string(); + const std::string written = (dir / "writer").string(); + const gemmi::SpaceGroup &sg = *gemmi::find_spacegroup_by_name("P 43 21 2"); + + // As rugnux writes them (RugnuxOutputFiles.cpp): WritePlacedModel at _refined, the maps' MTZ, and + // the blobs' JSON. + WritePlacedModel(TestModel(), TEST_CELL, sg, original + "_refined", logger); + MtzTable maps; + maps.unit_cell = TEST_CELL; + maps.space_group_name = sg.xhm(); + maps.datasets = {{"HKL_base", "HKL_base", "HKL_base", 0.0f}, {"refined", "refined", "refined", 1.0f}}; + for (const auto &[label, type] : std::vector>{ + {"H", 'H'}, {"K", 'H'}, {"L", 'H'}, {"FWT", 'F'}, {"PHWT", 'P'}, {"DELFWT", 'F'}, {"PHDELWT", 'P'}}) + maps.columns.push_back({label, type, type == 'H' ? 0 : 1, {}}); + for (int i = 1; i <= 30; i++) { + const float row[] = {static_cast(i % 5), static_cast(i % 7), static_cast(i), + 10.0f + i, static_cast((i * 37) % 360), 0.5f * i, static_cast(i)}; + for (size_t c = 0; c < maps.columns.size(); c++) + maps.columns[c].values.push_back(row[c]); + } + WriteMtz(maps, original + "_refined_maps.mtz"); + std::vector blobs(2); + blobs[0] = {.kind = DensityBlob::Kind::Ligand, .centroid = {12.25, -3.5, 40.125}, .peak_sigma = 7.3, + .volume = 41.2, .score = 160.4, .nearest_atom_distance = 3.1, .nearest_atom = "OG1 THR A42", + .polar_contacts = 2, .metal_contacts = 0, .anomalous_sigma = 0.0}; + blobs[1] = {.kind = DensityBlob::Kind::Water, .centroid = {1.0 / 3.0, 2.0, 3.0}, .peak_sigma = 4.4, + .volume = 2.5, .score = 9.1, .nearest_atom_distance = 2.8, .nearest_atom = "O HOH A101", + .polar_contacts = 3, .metal_contacts = 1, .anomalous_sigma = 1.25}; + { std::ofstream f(original + "_blobs.json"); f << DensityBlobsJson(DensityBlobsToMessage(blobs)); } + + // As the broker puts them in the END message (JFJochReceiverRugnux::AddToEndMessage). + EndMessage end{}; + end.refined_model = PlacedModelFromStructure(TestModel(), TEST_CELL, sg.xhm()); + end.refined_map_coefficients = maps; + end.blobs = DensityBlobsToMessage(blobs); + std::vector buffer(16 * 1024 * 1024); + CBORStream2Serializer serializer(buffer.data(), buffer.size()); + REQUIRE_NOTHROW(serializer.SerializeSequenceEnd(end)); + auto deserialized = CBORStream2Deserialize(buffer.data(), serializer.GetBufferSize()); + REQUIRE(deserialized); + REQUIRE(deserialized->end_message); + const EndMessage &out = *deserialized->end_message; + REQUIRE(out.refined_model); + REQUIRE(out.refined_map_coefficients); + REQUIRE(out.blobs); + CHECK(!out.placed_model); + CHECK(out.refined_model->x == end.refined_model->x); + CHECK(out.blobs->kind == std::vector{"LIGAND", "WATER"}); + CHECK(out.blobs->x == end.blobs->x); + CHECK(out.blobs->nearest_atom == end.blobs->nearest_atom); + CHECK(out.blobs->metal_contacts == end.blobs->metal_contacts); + CHECK(out.blobs->anomalous_sigma == end.blobs->anomalous_sigma); + + WriteRugnuxResultFiles(out, written, false, logger); + CHECK(ReadFile(written + "_refined_model.cif") == ReadFile(original + "_refined_model.cif")); + CHECK(ReadFile(written + "_refined_model.pdb") == ReadFile(original + "_refined_model.pdb")); + CHECK(ReadFile(written + "_refined_maps.mtz") == ReadFile(original + "_refined_maps.mtz")); + CHECK(ReadFile(written + "_blobs.json") == ReadFile(original + "_blobs.json")); + CHECK(!std::filesystem::exists(written + "_model.cif")); + CHECK(!std::filesystem::exists(written + "_blobs.json.tmp")); + const auto json = nlohmann::json::parse(ReadFile(written + "_blobs.json")); + REQUIRE(json["blobs"].size() == 2); + CHECK(json["blobs"][0]["kind"] == "LIGAND"); + CHECK(json["blobs"][1]["centroid"][0].get() == 1.0 / 3.0); + + // A blobs block whose fields disagree is reported, not written. + DensityBlobsMessage bad = *end.blobs; + bad.score.pop_back(); + CHECK_THROWS(DensityBlobsJson(bad)); + EndMessage bad_end{}; + bad_end.blobs = bad; + REQUIRE_NOTHROW(WriteRugnuxResultFiles(bad_end, written + "_bad", false, logger)); + CHECK(!std::filesystem::exists(written + "_bad_blobs.json")); + + std::filesystem::remove_all(dir); +} diff --git a/writer/WriteRugnuxFiles.cpp b/writer/WriteRugnuxFiles.cpp index a2e16b515..a9516a81c 100644 --- a/writer/WriteRugnuxFiles.cpp +++ b/writer/WriteRugnuxFiles.cpp @@ -6,10 +6,12 @@ #include #include #include +#include #include #include // cif::read_string #include // make_structure_from_block +#include #include "WriteModel.h" #include "../common/Logger.h" @@ -142,6 +144,28 @@ gemmi::Structure PlacedModelStructure(const PlacedModelMessage &model) { return st; } +std::string DensityBlobsJson(const DensityBlobsMessage &b) { + const size_t n = b.kind.size(); + for (const size_t size : {b.x.size(), b.y.size(), b.z.size(), b.peak_sigma.size(), b.volume.size(), b.score.size(), + b.nearest_atom.size(), b.nearest_atom_distance.size(), b.polar_contacts.size(), + b.metal_contacts.size(), b.anomalous_sigma.size()}) + if (size != n) + throw std::runtime_error(fmt::format("blobs: a field has {} entries for {} blobs", size, n)); + nlohmann::json list = nlohmann::json::array(); + for (size_t i = 0; i < n; i++) + list.push_back({{"kind", b.kind[i]}, + {"centroid", {b.x[i], b.y[i], b.z[i]}}, + {"peak_sigma", b.peak_sigma[i]}, + {"volume_A3", b.volume[i]}, + {"score_sigma_A3", b.score[i]}, + {"nearest_atom", b.nearest_atom[i]}, + {"nearest_atom_distance_A", b.nearest_atom_distance[i]}, + {"polar_contacts", b.polar_contacts[i]}, + {"metal_contacts", b.metal_contacts[i]}, + {"anomalous_sigma", b.anomalous_sigma[i]}}); + return nlohmann::json{{"blobs", list}}.dump(1) + "\n"; +} + namespace { double WallClockSeconds() { @@ -159,6 +183,20 @@ void RenameIntoPlace(const std::string &tmp, const std::string &final_name, bool logger.Info("{} written, renamed into place at wall clock {:.3f} s", final_name, WallClockSeconds()); } +// _model.cif (and .pdb) from `model`. +void WriteModelFiles(const PlacedModelMessage &model, const std::string &file_prefix, const std::string &name, + bool overwrite, Logger &logger) { + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(model.space_group_name); + if (sg == nullptr) + throw std::runtime_error("unknown space group '" + model.space_group_name + "'"); + // WritePlacedModel appends _model.cif and _model.pdb to the prefix it is given. + const std::string tmp_prefix = file_prefix + name + ".tmp"; + WritePlacedModel(PlacedModelStructure(model), model.unit_cell, *sg, tmp_prefix, logger); + for (const std::string suffix : {"_model.cif", "_model.pdb"}) + if (std::filesystem::exists(tmp_prefix + suffix)) + RenameIntoPlace(tmp_prefix + suffix, file_prefix + name + suffix, overwrite, logger); +} + } // namespace void WriteRugnuxResultFiles(const EndMessage &end, const std::string &file_prefix, bool overwrite, @@ -181,19 +219,37 @@ void WriteRugnuxResultFiles(const EndMessage &end, const std::string &file_prefi logger.Error("Cannot write {}_maps.mtz: {}", file_prefix, e.what()); } try { - if (end.placed_model) { - const PlacedModelMessage &model = *end.placed_model; - const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(model.space_group_name); - if (sg == nullptr) - throw std::runtime_error("unknown space group '" + model.space_group_name + "'"); - // WritePlacedModel appends _model.cif and _model.pdb to the prefix it is given. - const std::string tmp_prefix = file_prefix + ".tmp"; - WritePlacedModel(PlacedModelStructure(model), model.unit_cell, *sg, tmp_prefix, logger); - for (const std::string suffix : {"_model.cif", "_model.pdb"}) - if (std::filesystem::exists(tmp_prefix + suffix)) - RenameIntoPlace(tmp_prefix + suffix, file_prefix + suffix, overwrite, logger); - } + if (end.placed_model) + WriteModelFiles(*end.placed_model, file_prefix, "", overwrite, logger); } catch (const std::exception &e) { logger.Error("Cannot write {}_model.cif: {}", file_prefix, e.what()); } + try { + if (end.refined_model) + WriteModelFiles(*end.refined_model, file_prefix, "_refined", overwrite, logger); + } catch (const std::exception &e) { + logger.Error("Cannot write {}_refined_model.cif: {}", file_prefix, e.what()); + } + try { + if (end.refined_map_coefficients) { + WriteMtz(*end.refined_map_coefficients, file_prefix + "_refined_maps.mtz.tmp"); + RenameIntoPlace(file_prefix + "_refined_maps.mtz.tmp", file_prefix + "_refined_maps.mtz", overwrite, + logger); + } + } catch (const std::exception &e) { + logger.Error("Cannot write {}_refined_maps.mtz: {}", file_prefix, e.what()); + } + try { + if (end.blobs) { + const std::string text = DensityBlobsJson(*end.blobs); + { + std::ofstream f(file_prefix + "_blobs.json.tmp", std::ofstream::out | std::ofstream::trunc); + f.exceptions(std::ios::failbit | std::ios::badbit); + f << text; + } + RenameIntoPlace(file_prefix + "_blobs.json.tmp", file_prefix + "_blobs.json", overwrite, logger); + } + } catch (const std::exception &e) { + logger.Error("Cannot write {}_blobs.json: {}", file_prefix, e.what()); + } } diff --git a/writer/WriteRugnuxFiles.h b/writer/WriteRugnuxFiles.h index 6d0bdca38..9496d64d8 100644 --- a/writer/WriteRugnuxFiles.h +++ b/writer/WriteRugnuxFiles.h @@ -25,9 +25,13 @@ void WriteMtz(const MtzTable &table, const std::string &filename); // nor six per atom. gemmi::Structure PlacedModelStructure(const PlacedModelMessage &model); -// The three files above, from whichever blocks `end` carries, at .mtz, _maps.mtz and -// _model.cif (WritePlacedModel adds _model.pdb where a PDB can hold the cell). Each is written under a -// temporary name and renamed into place. An existing file is replaced only with `overwrite`; otherwise -// it is kept and the new one is not written. +// The text of _blobs.json. Throws where a field has another number of entries than kind. +std::string DensityBlobsJson(const DensityBlobsMessage &blobs); + +// The files above, from whichever blocks `end` carries, at .mtz, _maps.mtz and _model.cif +// (WritePlacedModel adds _model.pdb where a PDB can hold the cell), and the refinement's +// _refined_model.cif (.pdb), _refined_maps.mtz and _blobs.json. Each is written under a temporary name and +// renamed into place. An existing file is replaced only with `overwrite`; otherwise it is kept and the +// new one is not written. void WriteRugnuxResultFiles(const EndMessage &end, const std::string &file_prefix, bool overwrite, Logger &logger);