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 <prefix>_refined_model.cif (+ .pdb), <prefix>_refined_maps.mtz and <prefix>_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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit c65220d96a96d2b1897f3c6632c94f9379b670a3)
This commit is contained in:
leonarski_f committed 2026-10-26 15:50:40 +01:00
1 parent 2461d67e14
commit 2e296e6bed
13 files changed
+282 -36

No files matched your search

+17
View File
@@ -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 <prefix>_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<std::string> kind;
std::vector<double> x, y, z;
std::vector<double> peak_sigma, volume, score;
std::vector<std::string> nearest_atom;
std::vector<double> nearest_atom_distance;
std::vector<int32_t> polar_contacts, metal_contacts;
std::vector<double> 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<MtzTable> merged; // MergedMtzTable (image_analysis/WriteReflections.h)
std::optional<MtzTable> map_coefficients; // ModelMaps::coefficients (ModelValidation.h)
std::optional<PlacedModelMessage> placed_model;
// And of its refinement, where one ran: <file_prefix>_refined_model.cif, _refined_maps.mtz and _blobs.json.
std::optional<PlacedModelMessage> refined_model;
std::optional<MtzTable> refined_map_coefficients;
std::optional<DensityBlobsMessage> blobs;
};
struct MetadataMessage {
+4 -1
View File
@@ -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 `<file_prefix>.mtz`. An MTZ table (layout below) | |
| map_coefficients | object (optional) | Map coefficients of an in-broker analysis run, written as `<file_prefix>_maps.mtz`. An MTZ table (layout below) | |
| placed_model | object (optional) | Atomic model placed against the data, written as `<file_prefix>_model.cif` and `<file_prefix>_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 `<file_prefix>_refined_model.cif` and `<file_prefix>_refined_model.pdb`; same layout as `placed_model` | |
| refined_map_coefficients | object (optional) | Map coefficients of the refined model, written as `<file_prefix>_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 `<file_prefix>_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 |
|------------------|---------------|------------------------------------------------------------------------------------------------------|
+5 -3
View File
@@ -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 `<file_prefix>.mtz`,
`<file_prefix>_maps.mtz`, `<file_prefix>_model.cif` and `<file_prefix>_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
`<file_prefix>.mtz`, `<file_prefix>_maps.mtz`, `<file_prefix>_model.cif` and `<file_prefix>_model.pdb`,
and `<file_prefix>_refined_model.cif` (and `.pdb`), `<file_prefix>_refined_maps.mtz` and
`<file_prefix>_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.
@@ -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<float> tmp;
GetCBORFloatArray(value, tmp);
+22
View File
@@ -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));
@@ -42,6 +42,25 @@ const char *DensityBlobKindName(DensityBlob::Kind kind) {
return "UNASSIGNED";
}
DensityBlobsMessage DensityBlobsToMessage(const std::vector<DensityBlob> &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<DensityBlob> FindDensityBlobs(const gemmi::Grid<float> &difference, double rms, const gemmi::Structure &model,
double d_min, const gemmi::Grid<float> *anomalous, double anomalous_rms) {
@@ -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 <prefix>_blobs.json.
DensityBlobsMessage DensityBlobsToMessage(const std::vector<DensityBlob> &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
+3
View File
@@ -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);
+15 -1
View File
@@ -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();
}
}
+1 -15
View File
@@ -6,8 +6,6 @@
#include <filesystem>
#include <fstream>
#include <nlohmann/json.hpp>
#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<DensityBlob> &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);
}
+85
View File
@@ -11,11 +11,13 @@
#include <gemmi/mtz.hpp>
#include <gemmi/mmread_gz.hpp>
#include <gemmi/calculate.hpp> // count_atom_sites
#include <nlohmann/json.hpp>
#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 <prefix>_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<std::pair<const char *, char>>{
{"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<float>(i % 5), static_cast<float>(i % 7), static_cast<float>(i),
10.0f + i, static_cast<float>((i * 37) % 360), 0.5f * i, static_cast<float>(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<DensityBlob> 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<uint8_t> 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<std::string>{"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<double>() == 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);
}
+68 -12
View File
@@ -6,10 +6,12 @@
#include <chrono>
#include <stdexcept>
#include <filesystem>
#include <fstream>
#include <gemmi/mtz.hpp>
#include <gemmi/read_cif.hpp> // cif::read_string
#include <gemmi/mmcif.hpp> // make_structure_from_block
#include <nlohmann/json.hpp>
#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());
}
// <file_prefix><name>_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());
}
}
+8 -4
View File
@@ -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 <file_prefix>.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 <prefix>_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 <file_prefix>.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);