rugnux: read mmCIF models with --model, and say when a model could not be used
--model called gemmi::read_pdb unconditionally, so a coordinate file in mmCIF - which is how the PDB serves coordinates by default - was refused outright: Model validation: cannot read model 6G8A.cif: Incorrect file format (perhaps it is cif not pdb?) Worse than the refusal was what followed it. ValidateAgainstModel logged, returned an empty result, and the run went on to finish with exit code 0, no R-free, no maps and nothing in the report - indistinguishable from a run that was never given --model at all. A script that passed a .cif, checked the exit code and grepped for R-free simply got nothing back. Both read sites now go through read_structure_gz with CoorFormat::Detect, so PDB, mmCIF and mmJSON are all read, gzipped or not, and the format comes from the file's content rather than from its name. Detect is passed explicitly: GEMMI otherwise takes the extension and only falls back to the content when it does not recognise one, and a model arrives named however whoever produced it named it. That needed a part of GEMMI the vendored subset had trimmed away. mmread.hpp was already here but its CIF parser was not, so this restores what upstream GEMMI ships: read_cif/mmcif/json/mmread_gz, plus GEMMI's own copies of PEGTL (MIT, Colin Hirsch and Daniel Frey; 155 headers, 672 kB) and sajson (MIT, Chad Austin; one header). Both are listed in THIRD_PARTY_NOTICES.md with their own licence texts, PEGTL's collected by COLLECT.sh and sajson's kept by hand, its terms being a comment block rather than a file. json.cpp carries a one-line change: upstream keeps sajson at the root of its repository, outside the include tree this subset copies, so the include points at where the other bundled headers live. ModelValidationResult gains failure_reason, and a model that was asked for and could not be used now reaches the results report as WARNING: Model validation did not run: model bogus.pdb has no atoms or no unit cell in --mode mx and --mode scale alike. The run still succeeds - a merge that is good is not made bad by a model that was not - but it no longer does so quietly. One log line was making a claim it could not support: the indexing-ambiguity message said the ambiguity was "resolved against the supplied model" while being printed before the model is read, so it announced a resolution that a failed model never performed. It now says the model will be used; the reference branch, where the work really has already happened, keeps the past tense. Verified on the rotation test dataset with a deposited mmCIF, the same file gzipped, an mmCIF under a .pdb name, a PDB, and an unreadable file. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_016WmryXe8ASbNi632sUMfsa
This commit is contained in:
Claude Opus 5
@@ -0,0 +1,104 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <catch2/catch_all.hpp>
|
||||
|
||||
#include <filesystem>
|
||||
#include <fstream>
|
||||
|
||||
#include "../common/Logger.h"
|
||||
#include "../rugnux/ModelValidation.h"
|
||||
|
||||
namespace {
|
||||
// A synthetic P1 cell with two carbon atoms - enough for a reader to produce a Structure with
|
||||
// atoms, a cell and a space group, which is all these tests are about. Neutral by construction:
|
||||
// no real specimen's parameters are involved.
|
||||
const char *kMmcif = R"(data_test
|
||||
_cell.length_a 40.000
|
||||
_cell.length_b 50.000
|
||||
_cell.length_c 60.000
|
||||
_cell.angle_alpha 90.000
|
||||
_cell.angle_beta 90.000
|
||||
_cell.angle_gamma 90.000
|
||||
_symmetry.space_group_name_H-M 'P 1'
|
||||
loop_
|
||||
_atom_site.group_PDB
|
||||
_atom_site.id
|
||||
_atom_site.type_symbol
|
||||
_atom_site.label_alt_id
|
||||
_atom_site.label_atom_id
|
||||
_atom_site.label_comp_id
|
||||
_atom_site.label_asym_id
|
||||
_atom_site.label_seq_id
|
||||
_atom_site.Cartn_x
|
||||
_atom_site.Cartn_y
|
||||
_atom_site.Cartn_z
|
||||
_atom_site.occupancy
|
||||
_atom_site.B_iso_or_equiv
|
||||
ATOM 1 C . CA GLY A 1 10.000 12.000 14.000 1.00 20.00
|
||||
ATOM 2 C . CB GLY A 1 12.000 14.000 16.000 1.00 20.00
|
||||
)";
|
||||
|
||||
const char *kPdb =
|
||||
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n"
|
||||
"ATOM 1 CA GLY A 1 10.000 12.000 14.000 1.00 20.00 C\n"
|
||||
"ATOM 2 CB GLY A 1 12.000 14.000 16.000 1.00 20.00 C\n"
|
||||
"END\n";
|
||||
|
||||
std::string WriteTemp(const std::string &name, const char *content) {
|
||||
std::ofstream f(name);
|
||||
f << content;
|
||||
f.close();
|
||||
return name;
|
||||
}
|
||||
}
|
||||
|
||||
// --model used to call gemmi::read_pdb unconditionally, so a deposited model handed over as mmCIF -
|
||||
// which is how the PDB serves coordinates by default - was refused outright. Both formats now go
|
||||
// through the same reader, which decides on the file's content.
|
||||
TEST_CASE("ModelValidation_ReadsPdbAndMmcif", "[ModelValidation]") {
|
||||
Logger logger("ModelValidation_ReadsPdbAndMmcif");
|
||||
|
||||
const auto pdb = WriteTemp("model_validation_test.pdb", kPdb);
|
||||
const auto cif = WriteTemp("model_validation_test.cif", kMmcif);
|
||||
|
||||
const auto from_pdb = ModelReferenceIntensities(pdb, {}, {}, 4.0, logger);
|
||||
const auto from_cif = ModelReferenceIntensities(cif, {}, {}, 4.0, logger);
|
||||
|
||||
REQUIRE_FALSE(from_pdb.empty());
|
||||
REQUIRE_FALSE(from_cif.empty());
|
||||
// The same structure either way, so the same reflections come out of it.
|
||||
CHECK(from_cif.size() == from_pdb.size());
|
||||
|
||||
// The extension is not what decides: the same mmCIF under a .pdb name still reads.
|
||||
const auto misnamed = WriteTemp("model_validation_test_misnamed.pdb", kMmcif);
|
||||
CHECK_FALSE(ModelReferenceIntensities(misnamed, {}, {}, 4.0, logger).empty());
|
||||
|
||||
std::filesystem::remove(pdb);
|
||||
std::filesystem::remove(cif);
|
||||
std::filesystem::remove(misnamed);
|
||||
}
|
||||
|
||||
// A model that cannot be used must say why. Returning an empty result and logging was enough to lose
|
||||
// the fact entirely: the run finished successfully with no R-free and no maps, which is exactly what
|
||||
// a run without --model looks like.
|
||||
TEST_CASE("ModelValidation_UnusableModelGivesAReason", "[ModelValidation]") {
|
||||
Logger logger("ModelValidation_UnusableModelGivesAReason");
|
||||
|
||||
const auto empty_file = WriteTemp("model_validation_test_bogus.pdb", "not a coordinate file\n");
|
||||
|
||||
const auto result = ValidateAgainstModel({}, UnitCell{.a = 40, .b = 50, .c = 60,
|
||||
.alpha = 90, .beta = 90, .gamma = 90},
|
||||
empty_file, "", logger);
|
||||
CHECK_FALSE(result.ok);
|
||||
CHECK_FALSE(result.failure_reason.empty());
|
||||
CHECK_THAT(result.failure_reason, Catch::Matchers::ContainsSubstring(empty_file));
|
||||
|
||||
const auto missing = ValidateAgainstModel({}, UnitCell{.a = 40, .b = 50, .c = 60,
|
||||
.alpha = 90, .beta = 90, .gamma = 90},
|
||||
"model_validation_test_does_not_exist.pdb", "", logger);
|
||||
CHECK_FALSE(missing.ok);
|
||||
CHECK_FALSE(missing.failure_reason.empty());
|
||||
|
||||
std::filesystem::remove(empty_file);
|
||||
}
|
||||
Reference in New Issue
Block a user