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);