mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases 81 Wiki Activity
Files
40fcaa985cd3068cf0f3139ecdad353e18245938
Jungfraujoch/tests/ModelValidationTest.cpp
T
leonarski_fandClaude Opus 5.5 1f27fba1f3 Model validation expands strict NCS; the written cell follows the model's indexing 
- ModelValidation reads the model and, where it carries non-given NCS operators
  (_struct_ncs_oper 'generate' / MTRIX without the given flag), expands them
  (gemmi expand_ncs, vendored assembly.cpp from GEMMI 0.7.5, MPL-2.0 as the rest
  of gemmi_gph). Same in the battery's model_check.py and depdata_check.py.
  An icosahedral-virus entry deposited as 3 of the 45 chains of its ASU was
  scored against a fifteenth of its atoms: R-free 0.545 -> 0.297.
- The indexing the model chooses (AdoptModelFrame) now moves the cell with the
  indices (reseat_integrated in the main path; CellInBasis in the forecast,
  twin-file and stills/scale paths). A near-symmetry law (b within 2% of c, beta
  a degree off 90, a~b in P1) relabelled the reflections only, so the files were
  written with two axes swapped or beta flipped relative to the indices. Where
  the law moved the cell, the validation is remade on the moved cell.
  Measured on 18 open/in-house sets with models: p.hkl and MTZ reflection data
  bit-identical on all; cell and model R unchanged where the law is exact; the
  cell now matches the deposited one on the 6 pseudo-merohedral cases (orthorhombic
  a/c swap, three monoclinic beta flips, a P1 alpha/beta swap), R-free within
  +-0.001 there.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3
2026-10-10 13:35:49 +02:00

1553 lines
78 KiB
C++
Raw Blame History

// 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 <cmath>
#include <cstdio>
#include <cstring>
#include <filesystem>
#include <fstream>
#include <map>
#include <random>
#include <set>
#include <sstream>
#include <gemmi/mmread_gz.hpp>
#include <gemmi/fourier.hpp>
#include <gemmi/solmask.hpp>
#include "../common/Logger.h"
#include "../image_analysis/structure_refinement/ModelFFT.h"
#include "../image_analysis/structure_refinement/ModelGrid.h"
#include "../image_analysis/structure_refinement/ModelValidation.h"
#include "../image_analysis/structure_refinement/RigidBodyRefine.h"
#ifdef JFJOCH_USE_CUDA
#include "../image_analysis/structure_refinement/ModelStructureFactorsGPU.h"
#include "../image_analysis/structure_refinement/RigidBodyGPU.h"
#include "../image_analysis/structure_refinement/RigidBodyGPUEngine.h"
#include "../common/CUDAWrapper.h"
#endif
#include "../image_analysis/structure_refinement/SigmaA.h"
#include "../writer/WriteModel.h"
#include "../image_analysis/scale_merge/ReindexAmbiguity.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";
// A synthetic "protein": carbons filling one asymmetric unit of a small cell. The space group
// matters - in P1 the origin is free in all three directions, so |F| does not change when the
// whole content is translated and there is no translation to recover. No specimen is involved;
// the positions come out of a fixed seed.
const char *kCryst = "CRYST1 30.000 34.000 38.000 90.00 90.00 90.00 P 21 21 21 4\n";
// The same cell in a polar group: b is the unique axis, so the origin is free along y alone.
const char *kPolarCryst = "CRYST1 30.000 34.000 38.000 90.00 90.00 90.00 P 1 2 1 2\n";
std::string ClusterPdb(const char *cryst = kCryst) {
std::string pdb = cryst;
std::mt19937 rng(20260902);
std::uniform_real_distribution<double> x(2, 14), y(2, 16), z(2, 18);
char line[96];
for (int i = 1; i <= 150; i++) {
std::snprintf(line, sizeof line,
"ATOM %5d C UNK A 1 %8.3f%8.3f%8.3f 1.00 20.00 C\n",
i, x(rng), y(rng), z(rng));
pdb += line;
}
return pdb + "END\n";
}
// The same two atoms as kPdb, plus an anisotropic U, a partial occupancy and a water: a written
// model has to carry what the input carried, not only its coordinates.
const char *kPdbRich =
"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"
"ANISOU 1 CA GLY A 1 1000 1200 1400 100 200 300 C\n"
"ATOM 2 CB GLY A 1 12.000 14.000 16.000 0.60 25.00 C\n"
"HETATM 3 O HOH A 101 20.000 22.000 24.000 1.00 30.00 O\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 deposition can hold part of the asymmetric unit and the strict-NCS operators that generate the rest
// (MTRIX without the given flag). The model it describes is the expanded one: here the two atoms and
// their copy half a cell along a, which the explicit four-atom model spells out.
TEST_CASE("ModelValidation_ExpandsStrictNcs", "[ModelValidation]") {
Logger logger("ModelValidation_ExpandsStrictNcs");
const char *with_ncs =
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n"
"MTRIX1 1 1.000000 0.000000 0.000000 0.00000 1\n"
"MTRIX2 1 0.000000 1.000000 0.000000 0.00000 1\n"
"MTRIX3 1 0.000000 0.000000 1.000000 0.00000 1\n"
"MTRIX1 2 1.000000 0.000000 0.000000 20.00000\n"
"MTRIX2 2 0.000000 1.000000 0.000000 0.00000\n"
"MTRIX3 2 0.000000 0.000000 1.000000 0.00000\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";
const char *expanded =
"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"
"ATOM 3 CA GLY B 1 30.000 12.000 14.000 1.00 20.00 C\n"
"ATOM 4 CB GLY B 1 32.000 14.000 16.000 1.00 20.00 C\n"
"END\n";
const auto ncs_path = WriteTemp("model_validation_test_ncs.pdb", with_ncs);
const auto full_path = WriteTemp("model_validation_test_ncs_full.pdb", expanded);
const auto from_ncs = ModelReferenceIntensities(ncs_path, {}, {}, 4.0, logger);
const auto from_full = ModelReferenceIntensities(full_path, {}, {}, 4.0, logger);
REQUIRE(from_ncs.size() == from_full.size());
REQUIRE_FALSE(from_ncs.empty());
for (size_t i = 0; i < from_ncs.size(); i++) {
CHECK(from_ncs[i].h == from_full[i].h);
CHECK(from_ncs[i].I == Catch::Approx(from_full[i].I).margin(1e-3 * from_full[i].I + 1e-6));
}
std::filesystem::remove(ncs_path);
std::filesystem::remove(full_path);
}
// 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);
}
// A rigid-body step is only worth taking if it can find a shift it was not told about. The check is
// closed: the "observed" amplitudes are the model's own, so the answer is known to be zero shift.
TEST_CASE("ModelValidation_RigidBodyRecoversASmallShift", "[ModelValidation]") {
Logger logger("ModelValidation_RigidBodyRecoversASmallShift");
const auto path = WriteTemp("rigid_body_test.pdb", ClusterPdb().c_str());
gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
const gemmi::SpaceGroup *sg = st.find_spacegroup();
REQUIRE(sg != nullptr);
st.setup_cell_images();
// "Observed" amplitudes: the model's own structure factors, so the target's minimum is exactly
// where the model started.
const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger);
REQUIRE_FALSE(ref.empty());
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
fobs.unit_cell_ = st.cell;
fobs.spacegroup_ = sg;
for (const auto &r : ref)
fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}});
fobs.ensure_sorted();
const std::vector<gemmi::Position> original = ModelPositions(st.models[0]);
std::vector<gemmi::Position> displaced;
for (const gemmi::Position &p : original)
displaced.emplace_back(p.x + 0.40, p.y - 0.30, p.z + 0.20); // 0.54 A off
SetModelPositions(st.models[0], displaced);
const RigidBodyRefineResult result =
RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger);
CHECK(result.converged);
const std::vector<gemmi::Position> refined = ModelPositions(st.models[0]);
double before = 0, after = 0;
for (size_t i = 0; i < original.size(); i++) {
before += original[i].dist_sq(displaced[i]);
after += original[i].dist_sq(refined[i]);
}
before = std::sqrt(before / original.size());
after = std::sqrt(after / original.size());
logger.Info("Rigid-body test: rmsd from the truth {:.3f} A -> {:.3f} A", before, after);
CHECK(after < 0.2 * before);
std::filesystem::remove(path);
}
// A polar space group leaves the origin free along one direction: moving the whole cell content
// along it multiplies every structure factor by a phase and changes no amplitude, so the data
// cannot say where the body sits along it and the refinement must not pretend otherwise. The check
// is closed - the "observed" amplitudes are the model's own - and the model is displaced in all
// three directions at once, so the same run says both what is recovered and what is left alone.
TEST_CASE("ModelValidation_RigidBodyLeavesThePolarDirectionAlone", "[ModelValidation]") {
Logger logger("ModelValidation_RigidBodyLeavesThePolarDirectionAlone");
const auto path = WriteTemp("rigid_body_polar_test.pdb", ClusterPdb(kPolarCryst).c_str());
gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
const gemmi::SpaceGroup *sg = st.find_spacegroup();
REQUIRE(sg != nullptr);
st.setup_cell_images();
const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger);
REQUIRE_FALSE(ref.empty());
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
fobs.unit_cell_ = st.cell;
fobs.spacegroup_ = sg;
for (const auto &r : ref)
fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}});
fobs.ensure_sorted();
// b is the unique axis of P 1 2 1, so y is the free direction and x and z are determined.
const std::vector<gemmi::Position> original = ModelPositions(st.models[0]);
std::vector<gemmi::Position> displaced;
for (const gemmi::Position &p : original)
displaced.emplace_back(p.x + 0.35, p.y + 0.50, p.z - 0.30);
SetModelPositions(st.models[0], displaced);
const RigidBodyRefineResult result =
RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger);
CHECK(result.converged);
const std::vector<gemmi::Position> refined = ModelPositions(st.models[0]);
gemmi::Vec3 left;
for (size_t i = 0; i < original.size(); i++)
left += gemmi::Vec3(refined[i]) - gemmi::Vec3(original[i]);
left *= 1.0 / static_cast<double>(original.size());
logger.Info("Rigid-body polar test: left over ({:.3f}, {:.3f}, {:.3f}) A", left.x, left.y, left.z);
CHECK(std::fabs(left.x) < 0.10);
CHECK(std::fabs(left.z) < 0.10);
// Along b nothing was refined away, because there is nothing there to refine.
CHECK(left.y == Catch::Approx(0.50).margin(0.02));
std::filesystem::remove(path);
}
// The bulk solvent the placement is scored through has to stay inside the range a flat solvent model
// means anything in. gemmi's own scaler is an unbounded Levenberg-Marquardt and reaches b_sol of
// hundreds or thousands of A^2, which does not corrupt a reported number here but distorts the target
// that decides where the model goes, at every one of the hundreds of evaluations.
TEST_CASE("ModelValidation_RigidBodySolventStaysPhysical", "[ModelValidation]") {
Logger logger("ModelValidation_RigidBodySolventStaysPhysical");
const auto path = WriteTemp("rigid_body_solvent_test.pdb", ClusterPdb().c_str());
gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
const gemmi::SpaceGroup *sg = st.find_spacegroup();
REQUIRE(sg != nullptr);
st.setup_cell_images();
const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger);
REQUIRE_FALSE(ref.empty());
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
fobs.unit_cell_ = st.cell;
fobs.spacegroup_ = sg;
for (const auto &r : ref)
fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}});
fobs.ensure_sorted();
// Turned right around, so the model explains nothing: that is where an unbounded solvent fit
// has nothing to hold it and runs away.
std::vector<gemmi::Position> turned = ModelPositions(st.models[0]);
gemmi::Vec3 centre;
for (const gemmi::Position &p : turned)
centre += p;
centre *= 1.0 / static_cast<double>(turned.size());
for (gemmi::Position &p : turned)
p = gemmi::Position(centre - (gemmi::Vec3(p) - centre));
SetModelPositions(st.models[0], turned);
const RigidBodyRefineResult result =
RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger);
logger.Info("Rigid-body solvent test: k_sol {:.3f}, b_sol {:.1f} A^2", result.k_sol, result.b_sol);
CHECK(result.k_sol >= 0.10);
CHECK(result.k_sol <= 0.60);
CHECK(result.b_sol >= 10.0);
CHECK(result.b_sol <= 80.0);
std::filesystem::remove(path);
}
// sigma_A is what says how much of the model to believe, so the two ends of its range are what the
// weighting has to get right: a model that explains the data completely, and one that explains none
// of it.
TEST_CASE("ModelValidation_SigmaAWeightsFollowTheModelsAgreement", "[ModelValidation]") {
gemmi::UnitCell cell(40, 50, 60, 90, 90, 90);
std::mt19937 rng(12345);
std::normal_distribution<double> normal(0.0, 1.0);
auto weights = [&](bool agreeing) {
std::vector<SigmaAReflection> refl;
for (int i = 0; i < 2000; i++) {
SigmaAReflection r;
r.f_calc = std::fabs(normal(rng)) * 100;
r.f_obs = agreeing ? r.f_calc : std::fabs(normal(rng)) * 100;
r.inv_d2 = 0.01 + 0.2 * (i / 2000.0);
r.free = (i % 20) == 0; // the usual 5 %
refl.push_back(r);
}
return EstimateSigmaA(refl, cell);
};
const SigmaAResult perfect = weights(true);
const SigmaAResult useless = weights(false);
CHECK(perfect.mean_fom > 0.85);
CHECK(useless.mean_fom < 0.2);
CHECK(perfect.shells == 2); // 100 free reflections, 50 to a shell
// No free reflections to estimate on: the coefficients are left alone rather than weighted by a
// number that was never measured.
std::vector<SigmaAReflection> no_free;
for (int i = 0; i < 100; i++)
no_free.push_back({100.0, 100.0, 0.05, 1, false, false});
const SigmaAResult unweighted = EstimateSigmaA(no_free, cell);
CHECK(unweighted.weight.size() == no_free.size());
CHECK(unweighted.weight[0].m == 1.0);
CHECK(unweighted.weight[0].d == 1.0);
}
// The model rugnux scored has to reach disk, or a user overlaying their input model on rugnux's maps
// is wrong by the whole rigid-body shift. Two things have to hold: the file carries what the input
// carried, and it carries the cell and space group the reflection files beside it are written in -
// which, once --model has adopted the model's enantiomorph, is neither the data's original label nor
// necessarily the input model's.
// A depositor's model and rugnux's own reduction routinely describe the SAME lattice on different
// axes - I-centred against C-centred, unique axis b against c, a cyclic permutation of an
// orthorhombic cell - and the space-group NUMBER is identical in every one of those, so nothing that
// compares numbers can notice. Re-fractionalizing straight across such a pair scrambles the model,
// and the rigid body cannot undo it: six parameters about a centroid are not a change of basis. The
// check is closed - the data are the model's own structure factors - so a model that reaches the
// data's description of the lattice comes out at R near zero, and one that does not comes out near
// the 0.6 an unrelated structure gives.
TEST_CASE("ModelValidation_FindsTheDatasDescriptionOfTheLattice", "[ModelValidation]") {
Logger logger("ModelValidation_FindsTheDatasDescriptionOfTheLattice");
// The data: the cluster's own |F|^2 in its own frame, 30 x 34 x 38, P 21 21 21.
const auto data_model = WriteTemp("cob_data_test.pdb", ClusterPdb().c_str());
auto obs = ModelReferenceIntensities(data_model, {}, {}, 2.5, logger);
REQUIRE(obs.size() > 500);
for (size_t i = 0; i < obs.size(); i++) {
obs[i].F = std::sqrt(std::max(0.0f, obs[i].I));
obs[i].rfree_flag = (i % 20) == 0;
}
const UnitCell data_cell{.a = 30, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90};
const gemmi::SpaceGroup *data_sg = gemmi::find_spacegroup_by_name("P 21 21 21");
const std::string prefix = (std::filesystem::temp_directory_path() / "cob_test").string();
// The same structure as a deposition on permuted axes, a' = b, b' = c, c' = a. The cell is
// orthogonal, so the Cartesian coordinates permute with the axes.
std::string permuted = "CRYST1 34.000 38.000 30.000 90.00 90.00 90.00 P 21 21 21 4\n";
{
std::istringstream in(ClusterPdb());
std::string line;
char buf[96];
while (std::getline(in, line)) {
if (line.compare(0, 4, "ATOM") != 0)
continue;
std::snprintf(buf, sizeof buf, "%s%8.3f%8.3f%8.3f%s",
line.substr(0, 30).c_str(), std::stod(line.substr(38, 8)),
std::stod(line.substr(46, 8)), std::stod(line.substr(30, 8)),
line.substr(54).c_str());
permuted += buf;
permuted += "\n";
}
permuted += "END\n";
}
const auto other_setting = WriteTemp("cob_model_test.pdb", permuted.c_str());
const auto moved = ValidateAgainstModel(obs, data_cell, other_setting, prefix, logger, data_sg,
/*probe_indexing_ambiguity=*/false, 1, 1.0);
REQUIRE(moved.ok);
logger.Info("Change of basis test: took {} out of {}, R-work {:.4f}",
moved.change_of_basis_op.triplet(), moved.setting_as_read, moved.r_work);
CHECK_FALSE(moved.change_of_basis_op == gemmi::Op::identity());
CHECK(moved.r_work < 0.15);
// The isomorphous case, which is what nearly every run is: the model already describes the
// lattice the way the data do, the only operators that map its cell are its own symmetry, and
// nothing is scored or moved.
const auto same = ValidateAgainstModel(obs, data_cell, data_model, prefix, logger, data_sg,
false, 1, 1.0);
REQUIRE(same.ok);
CHECK(same.change_of_basis_op == gemmi::Op::identity());
CHECK(same.setting_as_read.empty());
CHECK(same.r_work < 0.15);
for (const char *suffix : {"_2fofc.ccp4", "_fofc.ccp4", "_anom.ccp4", "_maps.mtz"})
std::filesystem::remove(prefix + suffix);
std::filesystem::remove(data_model);
std::filesystem::remove(other_setting);
}
// An alternative indexing is settled by relabelling the DATA into the model's indexing, not by moving
// the model into the data's: that is what puts every dataset of one crystal form in one convention.
// Only where a reference has already fixed the data's indexing is the model moved instead.
TEST_CASE("ModelValidation_ReindexesTheDataIntoTheModelsIndexing", "[ModelValidation]") {
Logger logger("ModelValidation_ReindexesTheDataIntoTheModelsIndexing");
// Point group 4 on a tetragonal lattice (4/mmm): one alternative indexing.
const auto model = WriteTemp("reidx_model_test.pdb",
ClusterPdb("CRYST1 34.000 34.000 38.000 90.00 90.00 90.00 P 4 4\n").c_str());
auto own = ModelReferenceIntensities(model, {}, {}, 2.5, logger);
REQUIRE(own.size() > 500);
for (size_t i = 0; i < own.size(); i++) {
own[i].F = std::sqrt(std::max(0.0f, own[i].I));
own[i].rfree_flag = (i % 20) == 0;
}
const UnitCell cell{.a = 34, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90};
const gemmi::SpaceGroup *p4 = gemmi::find_spacegroup_by_name("P 4");
const auto laws = ReindexAmbiguityOperators(cell, *p4);
REQUIRE(laws.size() == 1);
const std::string prefix = (std::filesystem::temp_directory_path() / "reidx_test").string();
// The data as a run that picked the other indexing would have merged them.
const auto obs = ReindexReflections(own, laws.front());
const auto to_model = ValidateAgainstModel(obs, cell, model, prefix, logger, p4,
/*probe_indexing_ambiguity=*/true, 1, 1.0);
REQUIRE(to_model.ok);
CHECK(to_model.change_of_basis_op == gemmi::Op::identity());
CHECK(to_model.indexing_decided);
CHECK(to_model.indexing_op == laws.front());
CHECK(to_model.r_work < 0.15);
// A reference already fixed the data's indexing: the data stay, and the model is moved.
const auto fixed = ValidateAgainstModel(obs, cell, model, prefix, logger, p4, false, 1, 1.0);
REQUIRE(fixed.ok);
CHECK_FALSE(fixed.change_of_basis_op == gemmi::Op::identity());
CHECK(fixed.indexing_op == gemmi::Op::identity());
CHECK(fixed.r_work < 0.15);
// A near-perfect twin of the mis-indexed data, 45 % of it in the model's indexing: the model
// prefers that indexing by less than a model in a random orientation prefers one, so it has decided
// nothing and the data must keep the indexing they were merged in.
std::vector<MergedReflection> twinned = obs;
{
const gemmi::GroupOps gops = p4->operations();
const gemmi::ReciprocalAsu asu(p4);
auto key = [&](const gemmi::Miller &h) { return asu.to_asu(h, gops).first; };
std::map<gemmi::Miller, float> by_hkl;
for (const auto &r : obs)
by_hkl[key({{r.h, r.k, r.l}})] = r.I;
for (auto &r : twinned) {
const auto mate = by_hkl.find(key(laws.front().apply_to_hkl({{r.h, r.k, r.l}})));
REQUIRE(mate != by_hkl.end());
r.I = 0.55f * r.I + 0.45f * mate->second;
r.F = std::sqrt(std::max(0.0f, r.I));
}
}
const auto twin = ValidateAgainstModel(twinned, cell, model, prefix, logger, p4, true, 1, 1.0);
REQUIRE(twin.ok);
CHECK_FALSE(twin.indexing_decided);
CHECK(twin.indexing_op == gemmi::Op::identity());
CHECK(twin.change_of_basis_op == gemmi::Op::identity());
// Data already in the model's indexing: nothing moves, and there is nothing to arbitrate.
const auto same = ValidateAgainstModel(own, cell, model, prefix, logger, p4, true, 1, 1.0);
REQUIRE(same.ok);
CHECK(same.change_of_basis_op == gemmi::Op::identity());
CHECK(same.indexing_op == gemmi::Op::identity());
CHECK_FALSE(same.fit_tested);
CHECK(same.r_work < 0.15);
for (const char *suffix : {"_2fofc.ccp4", "_fofc.ccp4", "_anom.ccp4", "_maps.mtz"})
std::filesystem::remove(prefix + suffix);
std::filesystem::remove(model);
}
TEST_CASE("WriteModel_KeepsTheContentAndTakesTheGivenFrame", "[ModelValidation]") {
Logger logger("WriteModel_KeepsTheContentAndTakesTheGivenFrame");
const auto input = WriteTemp("write_model_test_input.pdb", kPdbRich);
gemmi::Structure st = gemmi::read_structure_gz(input, gemmi::CoorFormat::Detect);
// The model in the data's own group, which is the ordinary case: the frame it is written in is
// then wholly the caller's. The cell it arrives in (40/50/60) is not the data's and is replaced.
st.spacegroup_hm = "P 43 21 2";
// The tetragonal lysozyme cell and one enantiomorph of its group, standing in for what
// AdoptModelFrame settled.
const UnitCell data_cell{.a = 79, .b = 79, .c = 38, .alpha = 90, .beta = 90, .gamma = 90};
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 43 21 2");
REQUIRE(sg != nullptr);
WritePlacedModel(st, data_cell, *sg, "write_model_test", logger);
const std::string written = "write_model_test_model.cif";
REQUIRE(std::filesystem::exists(written));
const gemmi::Structure back = gemmi::read_structure_gz(written, gemmi::CoorFormat::Detect);
CHECK(back.cell.a == Catch::Approx(79.0));
CHECK(back.cell.c == Catch::Approx(38.0));
REQUIRE(back.find_spacegroup() != nullptr);
CHECK(back.find_spacegroup()->number == 96);
REQUIRE(back.models.size() == 1);
REQUIRE(back.models[0].chains.size() == 1);
const gemmi::Chain &chain = back.models[0].chains[0];
REQUIRE(chain.residues.size() == 2); // GLY and the water
REQUIRE(chain.residues[0].atoms.size() == 2);
REQUIRE(chain.residues[1].atoms.size() == 1);
// The coordinates are written as they stand - WritePlacedModel changes the label, not the atoms.
const gemmi::Atom &ca = chain.residues[0].atoms[0];
const gemmi::Atom &cb = chain.residues[0].atoms[1];
CHECK(ca.pos.x == Catch::Approx(10.0));
CHECK(ca.pos.z == Catch::Approx(14.0));
CHECK(cb.pos.y == Catch::Approx(14.0));
// B-factors, occupancies and the anisotropic U survive the round trip.
CHECK(ca.b_iso == Catch::Approx(20.0));
CHECK(cb.b_iso == Catch::Approx(25.0));
CHECK(cb.occ == Catch::Approx(0.60));
CHECK(ca.aniso.nonzero());
CHECK(ca.aniso.u11 == Catch::Approx(0.1000).margin(1e-4));
CHECK(ca.aniso.u12 == Catch::Approx(0.0100).margin(1e-4));
CHECK(chain.residues[1].name == "HOH");
std::filesystem::remove(input);
std::filesystem::remove(written);
}
// A model cannot be labelled with a group its own coordinates do not obey. Where the data were merged
// in a supergroup of the model's - an over-merge across a pseudo-symmetry operation, or the hand the
// model's fit did not earn - that operation would generate atoms the model does not contain, and a
// refinement program acts on it silently. The model keeps its own group and the caller is warned.
TEST_CASE("WriteModel_DoesNotStampAGroupTheCoordinatesDoNotObey", "[ModelValidation]") {
Logger logger("WriteModel_DoesNotStampAGroupTheCoordinatesDoNotObey");
const auto input = WriteTemp("write_model_sg_test_input.pdb", kPdbRich); // P 1
const gemmi::Structure st = gemmi::read_structure_gz(input, gemmi::CoorFormat::Detect);
const UnitCell data_cell{.a = 79, .b = 79, .c = 38, .alpha = 90, .beta = 90, .gamma = 90};
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 43 21 2");
REQUIRE(sg != nullptr);
WritePlacedModel(st, data_cell, *sg, "write_model_sg_test", logger);
const std::string written = "write_model_sg_test_model.cif";
REQUIRE(std::filesystem::exists(written));
const gemmi::Structure back = gemmi::read_structure_gz(written, gemmi::CoorFormat::Detect);
// The cell is still the data's - the coordinates do sit in it - but the group is the model's own.
CHECK(back.cell.a == Catch::Approx(79.0));
CHECK(back.cell.c == Catch::Approx(38.0));
REQUIRE(back.find_spacegroup() != nullptr);
CHECK(back.find_spacegroup()->number == 1);
std::filesystem::remove(input);
std::filesystem::remove(written);
std::filesystem::remove("write_model_sg_test_model.pdb");
}
// The other direction is not a mismatch at all: a model whose coordinates obey more symmetry than the
// data were merged in obeys the data's group too - the file simply holds more than one asymmetric unit
// of it, which is what a model in P1 always does. The reflection files' group is kept, as before.
TEST_CASE("WriteModel_KeepsTheDataGroupWhenItIsASubgroup", "[ModelValidation]") {
Logger logger("WriteModel_KeepsTheDataGroupWhenItIsASubgroup");
const auto input = WriteTemp("write_model_sub_test_input.pdb", kPdbRich);
gemmi::Structure st = gemmi::read_structure_gz(input, gemmi::CoorFormat::Detect);
st.spacegroup_hm = "P 43 21 2";
const UnitCell data_cell{.a = 79, .b = 79, .c = 38, .alpha = 90, .beta = 90, .gamma = 90};
const gemmi::SpaceGroup *p1 = gemmi::find_spacegroup_by_name("P 1");
REQUIRE(p1 != nullptr);
WritePlacedModel(st, data_cell, *p1, "write_model_sub_test", logger);
const std::string written = "write_model_sub_test_model.cif";
REQUIRE(std::filesystem::exists(written));
const gemmi::Structure back = gemmi::read_structure_gz(written, gemmi::CoorFormat::Detect);
REQUIRE(back.find_spacegroup() != nullptr);
CHECK(back.find_spacegroup()->number == 1);
std::filesystem::remove(input);
std::filesystem::remove(written);
std::filesystem::remove("write_model_sub_test_model.pdb");
}
// CC(model, data) has to follow where the signal actually is, or it cannot support the one-sided
// claim it exists for. The check is closed: the "observed" intensities are the model's own with
// Gaussian noise added, and how much noise is chosen per shell - almost none in the first, some in
// the second, enough to bury the signal in the third - so the answer is known before the run.
TEST_CASE("ModelValidation_CCModelFollowsTheSignalByShell", "[ModelValidation]") {
Logger logger("ModelValidation_CCModelFollowsTheSignalByShell");
const auto path = WriteTemp("cc_model_test.pdb", ClusterPdb().c_str());
auto obs = ModelReferenceIntensities(path, {}, {}, 2.5, logger);
REQUIRE(obs.size() > 1000);
// The shells the correlation is reported on, coarse to fine, and the noise each one gets as a
// multiple of the r.m.s. intensity of that shell. Nothing coarser than the first shell is kept:
// the reference intensities carry the bulk solvent at fixed constants while the validation fits
// its own, and below about 6 A that difference is a large part of |F| and would decorrelate a
// shell this test needs to be clean.
const std::vector<float> shells{5.0f, 3.2f, 2.5f};
const double noise[3] = {0.02, 1.0, 30.0};
auto shell_of = [&](float d) {
for (size_t s = 0; s < shells.size(); s++)
if (d > shells[s]) return s;
return shells.size();
};
std::erase_if(obs, [&](const MergedReflection &r) { return r.d > 6.0f || shell_of(r.d) >= shells.size(); });
REQUIRE(obs.size() > 500);
std::vector<double> sum_i2(shells.size(), 0.0);
std::vector<int> count(shells.size(), 0);
for (const auto &r : obs) {
sum_i2[shell_of(r.d)] += static_cast<double>(r.I) * r.I;
++count[shell_of(r.d)];
}
std::mt19937 rng(20260907);
std::normal_distribution<double> normal(0.0, 1.0);
for (size_t i = 0; i < obs.size(); i++) {
const size_t bin = shell_of(obs[i].d);
const double sd = noise[bin] * std::sqrt(sum_i2[bin] / count[bin]);
obs[i].I = static_cast<float>(obs[i].I + sd * normal(rng));
obs[i].sigma = static_cast<float>(std::max(1.0, sd));
obs[i].F = std::sqrt(std::max(0.0f, obs[i].I));
obs[i].rfree_flag = (i % 20) == 0;
}
const std::string prefix = (std::filesystem::temp_directory_path() / "cc_model_test").string();
const auto result =
ValidateAgainstModel(obs, UnitCell{.a = 30, .b = 34, .c = 38,
.alpha = 90, .beta = 90, .gamma = 90},
path, prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21"),
/*probe_indexing_ambiguity=*/false, 1, 1.0, shells);
REQUIRE(result.ok);
REQUIRE(result.cc_model_shells.size() == shells.size());
int n_total = 0;
for (size_t s = 0; s < shells.size(); s++) {
const auto &sh = result.cc_model_shells[s];
logger.Info("CC(model,data) {:.2f} A: {:.3f} on {} refl, {:+.1f} sigma",
sh.d_min, sh.cc, sh.n, sh.sigma);
CHECK(sh.d_min == shells[s]);
CHECK(sh.n > 20);
n_total += sh.n;
}
CHECK(n_total == result.cc_model_n);
// Essentially noiseless: the model is the data, so the correlation is high and hugely significant.
CHECK(result.cc_model_shells[0].cc > 0.9);
CHECK(result.cc_model_shells[0].sigma > 10.0);
// Noise at the shell's own r.m.s. still leaves plenty to see.
CHECK(result.cc_model_shells[1].cc > 0.25);
CHECK(result.cc_model_shells[1].sigma > 5.0);
// Buried: the shell must NOT come out significant, or the one-sided test would fire on noise.
CHECK(std::fabs(result.cc_model_shells[2].cc) < 0.15);
CHECK(std::fabs(result.cc_model_shells[2].sigma) < 4.0);
// No shells asked for, none reported: a run that did not measure it writes no key.
const auto no_shells =
ValidateAgainstModel(obs, UnitCell{.a = 30, .b = 34, .c = 38,
.alpha = 90, .beta = 90, .gamma = 90},
path, prefix, logger, gemmi::find_spacegroup_by_name("P 21 21 21"),
false, 1, 1.0);
CHECK(no_shells.ok);
CHECK(no_shells.cc_model_shells.empty());
for (const char *suffix : {"_2fofc.ccp4", "_fofc.ccp4", "_anom.ccp4", "_maps.mtz"})
std::filesystem::remove(prefix + suffix);
std::filesystem::remove(path);
}
// The model path's structure factors come from FFTW; they must be gemmi's own transform to float
// precision, in the same layout, so prepare_asu_data() reads the same reflections from either.
TEST_CASE("ModelValidation_MapToFPhiMatchesGemmi", "[ModelValidation]") {
// An even grid, and an odd one on every axis (FFTW and pocketfft split odd lengths differently).
const auto size = GENERATE(std::array<int, 3>{20, 24, 30}, std::array<int, 3>{15, 21, 27});
gemmi::Grid<float> map;
map.unit_cell.set(40.0, 50.0, 60.0, 90.0, 95.0, 90.0);
map.spacegroup = gemmi::find_spacegroup_by_name("P 1");
map.set_size(size[0], size[1], size[2]);
std::mt19937 rng(7);
std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
for (auto &x : map.data)
x = dist(rng);
gemmi::FPhiGrid<float> ref = gemmi::transform_map_to_f_phi(map, true);
gemmi::FPhiGrid<float> ours = MapToFPhi(map);
REQUIRE(ours.nu == ref.nu);
REQUIRE(ours.nv == ref.nv);
REQUIRE(ours.nw == ref.nw);
REQUIRE(ours.half_l == ref.half_l);
REQUIRE(ours.data.size() == ref.data.size());
double largest = 0, worst = 0;
for (size_t i = 0; i < ref.data.size(); i++) {
largest = std::max(largest, static_cast<double>(std::abs(ref.data[i])));
worst = std::max(worst, static_cast<double>(std::abs(ours.data[i] - ref.data[i])));
}
CHECK(worst <= 1e-5 * largest);
const auto a = ref.prepare_asu_data(4.0, 0, false, false, false);
const auto b = ours.prepare_asu_data(4.0, 0, false, false, false);
REQUIRE(a.v.size() == b.v.size());
for (size_t i = 0; i < a.v.size(); i++)
CHECK(a.v[i].hkl == b.v[i].hkl);
}
// The output maps come from FFTW too; each must be gemmi's own map to float precision, point for point
// in the same layout. Coefficients with arbitrary phases, negative indices, expanded by symmetry and
// Friedel, on grids that are odd along u and v (w is even by construction of a half-l grid).
TEST_CASE("ModelValidation_MapFromFPhiMatchesGemmi", "[ModelValidation]") {
const char *sg_name = GENERATE("P 1", "P 1 21 1", "P 21 21 21");
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(sg_name);
REQUIRE(sg != nullptr);
gemmi::AsuData<std::complex<float>> coef;
coef.unit_cell_.set(31.0, 43.0, 57.0, 90.0, sg->number == 4 ? 104.0 : 90.0, 90.0);
coef.spacegroup_ = sg;
const gemmi::ReciprocalAsu asu(sg);
const gemmi::GroupOps gops = sg->operations();
std::mt19937 rng(11);
std::uniform_real_distribution<float> amp(0.1f, 10.0f), phase(-3.14159f, 3.14159f);
for (int h = -7; h <= 7; h++)
for (int k = -9; k <= 9; k++)
for (int l = -11; l <= 11; l++) {
const gemmi::Op::Miller hkl{{h, k, l}};
if ((h == 0 && k == 0 && l == 0) || !asu.is_in(hkl) || gops.is_systematically_absent(hkl))
continue;
coef.v.push_back({hkl, std::polar(amp(rng), phase(rng))});
}
REQUIRE(coef.v.size() > 500);
// P 1 takes an odd grid on u and v; the screw axes need even factors.
const std::array<int, 3> size = sg->number == 1 ? std::array<int, 3>{17, 21, 26}
: std::array<int, 3>{18, 24, 26};
gemmi::FPhiGrid<float> grid = gemmi::get_f_phi_on_grid<float>(coef, size, true);
grid.data[grid.index_n(2, -3, 4)] = std::complex<float>(1.0f, NAN); // a missing coefficient
const gemmi::Grid<float> ours = MapFromFPhi(grid);
const gemmi::Grid<float> ref = gemmi::transform_f_phi_grid_to_map(gemmi::FPhiGrid<float>(grid));
REQUIRE(ours.nu == ref.nu);
REQUIRE(ours.nv == ref.nv);
REQUIRE(ours.nw == ref.nw);
REQUIRE(ours.axis_order == ref.axis_order);
REQUIRE(ours.spacegroup == ref.spacegroup);
REQUIRE(ours.data.size() == ref.data.size());
double largest = 0, worst = 0;
for (size_t i = 0; i < ref.data.size(); i++) {
REQUIRE(std::isfinite(ours.data[i]));
largest = std::max(largest, static_cast<double>(std::abs(ref.data[i])));
worst = std::max(worst, static_cast<double>(std::abs(ours.data[i] - ref.data[i])));
}
CHECK(largest > 0);
CHECK(worst <= 1e-5 * largest);
}
// Map -> coefficients -> map is the identity (the V/N and 1/V scales cancel the unnormalised
// transforms), on a grid odd along u and v.
TEST_CASE("ModelValidation_ModelFFTRoundTrip", "[ModelValidation]") {
gemmi::Grid<float> map;
map.unit_cell.set(35.0, 45.0, 55.0, 80.0, 95.0, 105.0);
map.spacegroup = gemmi::find_spacegroup_by_name("P 1");
map.set_size(15, 21, 28);
std::mt19937 rng(3);
std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
for (auto &x : map.data)
x = dist(rng);
const gemmi::Grid<float> back = MapFromFPhi(MapToFPhi(map));
REQUIRE(back.nu == map.nu);
REQUIRE(back.nv == map.nv);
REQUIRE(back.nw == map.nw);
double worst = 0;
for (size_t i = 0; i < map.data.size(); i++)
worst = std::max(worst, static_cast<double>(std::abs(back.data[i] - map.data[i])));
CHECK(worst <= 1e-5);
}
// The Jacobian's six columns are evaluated in parallel; the placement must be the serial one, bit for
// bit.
TEST_CASE("ModelValidation_RigidBodySameOnAnyNumberOfThreads", "[ModelValidation]") {
Logger logger("ModelValidation_RigidBodySameOnAnyNumberOfThreads");
const auto path = WriteTemp("rigid_body_threads_test.pdb", ClusterPdb().c_str());
gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
const gemmi::SpaceGroup *sg = st.find_spacegroup();
REQUIRE(sg != nullptr);
st.setup_cell_images();
const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger);
REQUIRE_FALSE(ref.empty());
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
fobs.unit_cell_ = st.cell;
fobs.spacegroup_ = sg;
for (const auto &r : ref)
fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}});
fobs.ensure_sorted();
std::vector<gemmi::Position> displaced;
for (const gemmi::Position &p : ModelPositions(st.models[0]))
displaced.emplace_back(p.x + 0.40, p.y - 0.30, p.z + 0.20);
gemmi::Model serial = st.models[0], parallel = st.models[0];
SetModelPositions(serial, displaced);
SetModelPositions(parallel, displaced);
const RigidBodyRefineResult r1 = RefineRigidBody(serial, st.cell, *sg, fobs, 3.0, logger, 1);
const RigidBodyRefineResult r6 = RefineRigidBody(parallel, st.cell, *sg, fobs, 3.0, logger, 6);
CHECK(r1.evaluations == r6.evaluations);
CHECK(r1.angle_deg == r6.angle_deg);
CHECK(r1.shift_A == r6.shift_A);
const auto p1 = ModelPositions(serial), p6 = ModelPositions(parallel);
REQUIRE(p1.size() == p6.size());
for (size_t i = 0; i < p1.size(); i++)
CHECK((p1[i].x == p6[i].x && p1[i].y == p6[i].y && p1[i].z == p6[i].z));
std::filesystem::remove(path);
}
// The rigid body puts each probe placement on the grid with its own parallel copy of gemmi's gridding
// (density, solvent mask and their symmetrization). It must give gemmi's grids bit for bit, on any
// number of threads, for isotropic and anisotropic atoms and for groups with and without centring.
TEST_CASE("ModelValidation_ParallelGriddingMatchesGemmi", "[ModelValidation]") {
const char *crysts[] = {
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n",
"CRYST1 40.000 50.000 60.000 90.00 100.00 90.00 C 1 2 1 4\n",
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 21 21 21 4\n",
"CRYST1 60.000 60.000 60.000 90.00 90.00 90.00 I 2 3 24\n",
"CRYST1 80.000 80.000 80.000 90.00 90.00 90.00 F 41 3 2 96\n",
};
for (const char *cryst : crysts) {
const auto path = WriteTemp("parallel_gridding_test.pdb", ClusterPdb(cryst).c_str());
gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
std::filesystem::remove(path);
const gemmi::SpaceGroup *sg = st.find_spacegroup();
REQUIRE(sg != nullptr);
int i = 0;
for (gemmi::Chain &ch : st.models[0].chains)
for (gemmi::Residue &r : ch.residues)
for (gemmi::Atom &a : r.atoms)
if (i++ % 3 == 0)
a.aniso = {0.30f, 0.25f, 0.20f, 0.02f, -0.01f, 0.03f};
const auto same = [](const gemmi::Grid<float> &a, const gemmi::Grid<float> &b) {
return a.data.size() == b.data.size() &&
std::memcmp(a.data.data(), b.data.data(), a.data.size() * sizeof(float)) == 0;
};
for (double d_min : {6.0, 3.5}) {
gemmi::DensityCalculator<gemmi::IT92<float>, float> gemmi_dc;
gemmi_dc.d_min = d_min;
gemmi_dc.rate = 1.5;
gemmi_dc.grid.unit_cell = st.cell;
gemmi_dc.grid.spacegroup = sg;
gemmi_dc.set_refmac_compatible_blur(st.models[0]);
gemmi_dc.put_model_density_on_grid(st.models[0]);
gemmi::Grid<float> gemmi_mask;
gemmi_mask.unit_cell = st.cell;
gemmi_mask.spacegroup = sg;
gemmi_mask.set_size_from_spacing(gemmi_dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
gemmi::SolventMasker(gemmi::AtomicRadiiSet::Refmac).put_mask_on_grid(gemmi_mask, st.models[0]);
for (size_t nthreads : {1, 4}) {
const std::vector<size_t> leaders = OrbitLeaders(gemmi_mask, nthreads);
gemmi::DensityCalculator<gemmi::IT92<float>, float> dc;
dc.d_min = d_min;
dc.rate = 1.5;
dc.grid.unit_cell = st.cell;
dc.grid.spacegroup = sg;
dc.set_refmac_compatible_blur(st.models[0]);
PutModelDensityOnGrid(dc, st.models[0], leaders, nthreads);
CHECK(same(dc.grid, gemmi_dc.grid));
gemmi::Grid<float> mask;
mask.unit_cell = st.cell;
mask.spacegroup = sg;
mask.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
PutMaskOnGrid(mask, st.models[0], leaders, nthreads);
CHECK(same(mask, gemmi_mask));
}
}
}
}
namespace {
// The five groups of ModelValidation_ParallelGriddingMatchesGemmi: none, a centring, screws, a
// cubic body centring and a cubic face centring with a diamond glide.
const char *kRigidBodyCrysts[] = {
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n",
"CRYST1 40.000 50.000 60.000 90.00 100.00 90.00 C 1 2 1 4\n",
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 21 21 21 4\n",
"CRYST1 60.000 60.000 60.000 90.00 90.00 90.00 I 2 3 24\n",
"CRYST1 80.000 80.000 80.000 90.00 90.00 90.00 F 41 3 2 96\n",
};
// ClusterPdb() in `cryst`, every third atom anisotropic.
gemmi::Structure AnisoCluster(const char *cryst) {
const auto path = WriteTemp("rigid_body_composition_test.pdb", ClusterPdb(cryst).c_str());
gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
std::filesystem::remove(path);
st.setup_cell_images();
int i = 0;
for (gemmi::Chain &ch : st.models[0].chains)
for (gemmi::Residue &r : ch.residues)
for (gemmi::Atom &a : r.atoms)
if (i++ % 3 == 0)
a.aniso = {0.30f, 0.25f, 0.20f, 0.02f, -0.01f, 0.03f};
return st;
}
// The zone's density calculator, as the rigid body sets it up for `model`.
gemmi::DensityCalculator<gemmi::IT92<float>, float> ZoneDensity(const gemmi::Structure &st,
const gemmi::Model &model, double d_min) {
gemmi::DensityCalculator<gemmi::IT92<float>, float> dc;
dc.d_min = d_min;
dc.rate = 1.5;
dc.grid.unit_cell = st.cell;
dc.grid.spacegroup = st.find_spacegroup();
dc.set_refmac_compatible_blur(model);
return dc;
}
}
// The rigid body's Fcalc is composed from the transform of one copy of the model instead of being
// taken from the symmetrized grid. The two are the same sum rearranged, so they must agree to
// rounding, and must give a value to exactly the reflections prepare_asu_data() does - the
// systematic absences of the group and a reflection sitting exactly on d_min (the (10 0 0) of the
// 60 A cubic cell at 6 A) are left out by both.
TEST_CASE("ModelValidation_RigidBodyP1FcalcMatchesSymmetrized", "[ModelValidation]") {
for (const char *cryst : kRigidBodyCrysts) {
const gemmi::Structure st = AnisoCluster(cryst);
const gemmi::Model &model = st.models[0];
for (double d_min : {6.0, 3.5}) {
auto sym = ZoneDensity(st, model, d_min);
sym.put_model_density_on_grid(model);
gemmi::FPhiGrid<float> sym_f = MapToFPhi(sym.grid);
const auto ref = sym_f.prepare_asu_data(d_min, sym.blur, false, false, false);
// Every index of the ASU in the grid's box, whatever its resolution, absences and 000 too.
std::vector<gemmi::Miller> candidates;
for (const auto &hv : sym_f.prepare_asu_data(0, 0, true, true, false).v)
candidates.push_back(hv.hkl);
auto copy = ZoneDensity(st, model, d_min);
PutModelDensityOnGrid(copy, model, {}, 4);
const SymmetryComposition composition(copy.grid, d_min, candidates);
std::vector<std::complex<double>> f;
composition.Compose(MapToFPhi(copy.grid), copy.blur, f, nullptr, 4);
REQUIRE(composition.Hkl().size() == ref.v.size());
double mean = 0, worst = 0;
for (size_t i = 0; i < ref.v.size(); i++) {
REQUIRE(composition.Hkl()[i] == ref.v[i].hkl);
mean += std::abs(ref.v[i].value) / static_cast<double>(ref.v.size());
worst = std::max(worst, std::abs(f[i] - std::complex<double>(ref.v[i].value)));
}
INFO(cryst << " at " << d_min << " A: worst " << worst / mean << " of the mean |F|");
CHECK(worst <= 1e-4 * mean);
}
}
}
// The derivative of the composed Fcalc with respect to a translation of the body is a phase factor
// per term, exact - against a central difference of the regridded, recomposed Fcalc. Along a
// direction the origin is free in (all three in P1, b in C2) the amplitude does not change at all.
TEST_CASE("ModelValidation_RigidBodyTranslationDerivativeIsExact", "[ModelValidation]") {
for (const char *cryst : kRigidBodyCrysts) {
const gemmi::Structure st = AnisoCluster(cryst);
const double d_min = 3.5;
const gemmi::DensityCalculator<gemmi::IT92<float>, float> zone = ZoneDensity(st, st.models[0], d_min);
gemmi::Grid<float> grid;
grid.unit_cell = zone.grid.unit_cell;
grid.spacegroup = zone.grid.spacegroup;
grid.set_size_from_spacing(zone.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
std::vector<gemmi::Miller> candidates;
const gemmi::ReciprocalAsu asu(grid.spacegroup);
for (int h = -25; h <= 25; h++)
for (int k = -25; k <= 25; k++)
for (int l = -25; l <= 25; l++)
if (asu.is_in({{h, k, l}}))
candidates.push_back({{h, k, l}});
const SymmetryComposition composition(grid, d_min, candidates);
REQUIRE(composition.Hkl().size() > 200);
auto fcalc = [&](const gemmi::Vec3 &shift, std::vector<std::array<std::complex<double>, 3>> *df_dt) {
gemmi::Model model = st.models[0];
for (gemmi::Chain &ch : model.chains)
for (gemmi::Residue &r : ch.residues)
for (gemmi::Atom &a : r.atoms)
a.pos += gemmi::Position(shift);
auto dc = ZoneDensity(st, model, d_min);
PutModelDensityOnGrid(dc, model, {}, 4);
std::vector<std::complex<double>> f;
composition.Compose(MapToFPhi(dc.grid), dc.blur, f, df_dt, 4);
return f;
};
std::vector<std::array<std::complex<double>, 3>> df_dt;
const std::vector<std::complex<double>> f = fcalc({}, &df_dt);
// A five-point difference, whose truncation error goes as eps^4: the grid is float, and at the
// 1e-3 A of a plain central difference its rounding alone reads as 3e-3 of the derivative. What
// is left, about 5e-4, is the sampling - an atom gridded after a move is the moved atom only to
// the accuracy of the grid - which the exact derivative does not have and the difference does.
const double eps = 0.05;
const bool p1 = st.find_spacegroup()->number == 1, c2 = st.find_spacegroup()->number == 5;
for (int k = 0; k < 3; k++) {
gemmi::Vec3 e;
e.at(k) = eps;
const auto p_1 = fcalc(e, nullptr), m_1 = fcalc(-e, nullptr);
const auto p_2 = fcalc(2 * e, nullptr), m_2 = fcalc(-2 * e, nullptr);
double diff = 0, norm = 0, amplitude = 0;
for (size_t m = 0; m < f.size(); m++) {
const std::complex<double> numeric = (8.0 * (p_1[m] - m_1[m]) - (p_2[m] - m_2[m])) / (12 * eps);
diff += std::norm(numeric - df_dt[m][k]);
norm += std::norm(df_dt[m][k]);
amplitude += gemmi::sq(std::real(std::conj(f[m]) * df_dt[m][k]) / std::abs(f[m]));
}
INFO(cryst << " axis " << k << ": relative error " << std::sqrt(diff / norm)
<< ", amplitude part " << std::sqrt(amplitude / norm));
CHECK(std::sqrt(diff / norm) <= 1e-3);
if (p1 || (c2 && k == 1))
CHECK(std::sqrt(amplitude / norm) <= 1e-5);
}
}
}
#ifdef JFJOCH_USE_CUDA
namespace {
// A pool of `engines`, or a SKIP where there is none because the card is busy - over half its memory
// taken by something else, where the pool falls back to the CPU by design. Anywhere else no pool fails.
std::unique_ptr<RigidBodyGPUPool> TestPool(const gemmi::Model &model, const gemmi::UnitCell &cell,
const gemmi::SpaceGroup &sg, double d_min, size_t observations,
size_t engines, Logger &logger) {
auto pool = RigidBodyGPUPool::Create(model, cell, sg, d_min, observations, engines, logger);
if (!pool) {
size_t free = 0, total = 0;
RigidBodyGPUEngine::MemoryInfo(free, total);
if (free < total / 2)
SKIP("No GPU engine: the card is busy (" << free / 1000000 << " of " << total / 1000000 << " MB free)");
}
REQUIRE(pool);
return pool;
}
}
#endif
namespace {
// The rigid body's Jacobian at q against a central difference of its own residuals - which
// re-fit the scale at every evaluation - with the bulk-solvent mask held, as the Jacobian holds
// it. Per column (those in `columns`): the cosine between the two and the ratio of their norms.
// Also logged, not checked: how far the difference moves once the mask is let move with the body,
// which is what holding it costs.
// The model's own amplitudes to `d_min`: "observed" data whose minimum is where the model is.
gemmi::AsuData<gemmi::ValueSigma<float>> OwnAmplitudes(const char *cryst, const gemmi::Structure &st,
double d_min, Logger &logger) {
const auto path = WriteTemp("rigid_body_own_amplitudes.pdb", ClusterPdb(cryst).c_str());
const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger);
std::filesystem::remove(path);
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
fobs.unit_cell_ = st.cell;
fobs.spacegroup_ = st.find_spacegroup();
for (const auto &r : ref)
if (st.cell.calculate_d({{r.h, r.k, r.l}}) >= d_min)
fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}});
fobs.ensure_sorted();
return fobs;
}
// The rigid body's target on the CPU, or on an engine of `gpu` where it is given.
std::unique_ptr<RigidBodyTargetBase> MakeTarget(gemmi::Model &model, const gemmi::Structure &st,
RigidBodyGPUPool *gpu) {
#ifdef JFJOCH_USE_CUDA
if (gpu != nullptr)
return std::make_unique<RigidBodyTargetGPU>(*gpu, model, st.cell, *st.find_spacegroup(), 4);
#endif
return std::make_unique<RigidBodyTarget>(model, st.cell, *st.find_spacegroup(), 4);
}
void CheckRigidBodyJacobian(const char *cryst, const double q0[6], const std::vector<int> &columns,
double min_cosine, double max_norm_error, bool gpu = false) {
Logger logger("CheckRigidBodyJacobian");
const gemmi::Structure st = AnisoCluster(cryst);
const gemmi::SpaceGroup *sg = st.find_spacegroup();
// "Observed" amplitudes: the model's own, as placed at q = 0.
const auto path = WriteTemp("rigid_body_jacobian_test.pdb", ClusterPdb(cryst).c_str());
const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger);
std::filesystem::remove(path);
REQUIRE_FALSE(ref.empty());
for (double zone : {6.0, 3.5}) {
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
fobs.unit_cell_ = st.cell;
fobs.spacegroup_ = sg;
for (const auto &r : ref)
if (st.cell.calculate_d({{r.h, r.k, r.l}}) >= zone)
fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}});
fobs.ensure_sorted();
gemmi::Model model = st.models[0];
RigidBodyGPUPool *pool = nullptr;
#ifdef JFJOCH_USE_CUDA
std::unique_ptr<RigidBodyGPUPool> engines;
if (gpu) {
engines = TestPool(model, st.cell, *sg, zone, fobs.v.size(), 1, logger);
pool = engines.get();
}
#else
(void) gpu;
#endif
const std::unique_ptr<RigidBodyTargetBase> target_backend = MakeTarget(model, st, pool);
RigidBodyTargetBase &target = *target_backend;
target.SetZone(fobs, zone);
const size_t n = target.NumObservations();
std::vector<double> r(n), jacobian(n * 6);
REQUIRE(target.Residuals(q0, r.data()));
REQUIRE(target.Jacobian(q0, jacobian.data()));
auto difference = [&](int j, double h) {
double qp[6], qm[6];
std::copy(q0, q0 + 6, qp);
std::copy(q0, q0 + 6, qm);
qp[j] += h;
qm[j] -= h;
std::vector<double> rp(n), rm(n), d(n);
REQUIRE(target.Residuals(qp, rp.data()));
REQUIRE(target.Residuals(qm, rm.data()));
for (size_t i = 0; i < n; i++)
d[i] = (rp[i] - rm[i]) / (2 * h);
return d;
};
std::vector<std::vector<double>> held(6), moving(6);
target.hold_mask = true;
REQUIRE(target.Residuals(q0, r.data())); // the mask the difference holds is q0's
for (int j : columns)
held[j] = difference(j, 0.005);
// The mask is binary on the grid, so with it moving the difference is taken at the step the
// Jacobian had when every column moved it (0.01 of the zone's resolution): at 0.005 A it
// mostly measures which grid points a few atoms happened to cross.
target.hold_mask = false;
for (int j : columns)
moving[j] = difference(j, 0.01 * zone);
double largest = 0;
for (int j : columns) {
double s = 0;
for (size_t i = 0; i < n; i++)
s += gemmi::sq(held[j][i]);
largest = std::max(largest, std::sqrt(s));
}
for (int j : columns) {
double dot = 0, ours = 0, theirs = 0, dot_moving = 0, norm_moving = 0;
for (size_t i = 0; i < n; i++) {
dot += jacobian[i * 6 + j] * held[j][i];
ours += gemmi::sq(jacobian[i * 6 + j]);
theirs += gemmi::sq(held[j][i]);
dot_moving += moving[j][i] * held[j][i];
norm_moving += gemmi::sq(moving[j][i]);
}
ours = std::sqrt(ours);
theirs = std::sqrt(theirs);
norm_moving = std::sqrt(norm_moving);
const double cosine = dot / (ours * theirs);
logger.Info("{} at {:.1f} A, column {}: cosine {:.5f}, norm {:.4f} of the difference's; "
"with the mask moving, cosine {:.4f} and norm {:.4f}",
std::string(cryst).substr(55, 11), zone, j, cosine, ours / theirs,
dot_moving / (norm_moving * theirs), norm_moving / theirs);
INFO(cryst << " at " << zone << " A, column " << j << ": cosine " << cosine
<< ", norm ratio " << ours / theirs);
if (theirs < 1e-3 * largest) {
// A direction the origin is free in: |F| does not change, so neither may the column.
CHECK(ours < 1e-2 * largest);
continue;
}
CHECK(cosine >= min_cosine);
CHECK(std::fabs(ours / theirs - 1) <= max_norm_error);
}
}
}
}
// The rotation columns at a placement already rotated (about 1.4 deg, so the step is taken from a
// rotated body and not from the model as read) with anisotropic atoms, whose U the placement does not
// turn with the body - the columns must differentiate exactly the function the residuals evaluate.
TEST_CASE("ModelValidation_RigidBodyRotationDerivativeMatchesDifferences", "[ModelValidation]") {
const double q0[6] = {0.12, -0.09, 0.07, 0.15, -0.10, 0.05};
for (const char *cryst : kRigidBodyCrysts)
CheckRigidBodyJacobian(cryst, q0, {0, 1, 2}, 0.99, 0.03);
}
// The whole Jacobian - analytic translation, forward-difference rotation, the scale re-fit folded in by
// projection - against the central difference of the residuals it describes. The projection is
// Kaufman's, which leaves out a term that grows with the residuals (the scale's derivatives depend on
// the placement too, weighted by how badly the model fits): with the body 0.07 A off it agrees with the
// difference to about 2e-3 in cosine, with it 0.35 A off to about 2e-2, and the bounds follow.
TEST_CASE("ModelValidation_RigidBodyJacobianMatchesNumericScaleRefit", "[ModelValidation]") {
const double close[6] = {0.0, 0.0, 0.0, 0.05, -0.04, 0.03};
const double far[6] = {0.0, 0.0, 0.0, 0.25, -0.20, 0.15};
for (const char *cryst : kRigidBodyCrysts) {
CheckRigidBodyJacobian(cryst, close, {0, 1, 2, 3, 4, 5}, 0.99, 0.02);
CheckRigidBodyJacobian(cryst, far, {0, 1, 2, 3, 4, 5}, 0.975, 0.03);
}
}
// The null checks where a replicate ENDED against the orientations equivalent to the model's, from the
// rotation the rigid body reports - which must be the rotation it applied: every atom's offset from the
// centroid after the refinement is that rotation of its offset before.
TEST_CASE("ModelValidation_RigidBodyReportsTheRotationItApplied", "[ModelValidation]") {
Logger logger("ModelValidation_RigidBodyReportsTheRotationItApplied");
const auto path = WriteTemp("rigid_body_rotation_test.pdb", ClusterPdb().c_str());
gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
const gemmi::SpaceGroup *sg = st.find_spacegroup();
REQUIRE(sg != nullptr);
st.setup_cell_images();
const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger);
std::filesystem::remove(path);
REQUIRE_FALSE(ref.empty());
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
fobs.unit_cell_ = st.cell;
fobs.spacegroup_ = sg;
for (const auto &r : ref)
fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}});
fobs.ensure_sorted();
// Turned 3 deg about z through the centroid, so there is a rotation to take back.
std::vector<gemmi::Position> turned = ModelPositions(st.models[0]);
gemmi::Vec3 centre;
for (const gemmi::Position &p : turned)
centre += p;
centre *= 1.0 / static_cast<double>(turned.size());
const double a = 3.0 * PI / 180.0;
const gemmi::Mat33 rz(std::cos(a), -std::sin(a), 0, std::sin(a), std::cos(a), 0, 0, 0, 1);
for (gemmi::Position &p : turned)
p = gemmi::Position(rz.multiply(gemmi::Vec3(p) - centre) + centre);
SetModelPositions(st.models[0], turned);
const RigidBodyRefineResult result = RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger);
REQUIRE(result.angle_deg > 1.0);
const std::vector<gemmi::Position> refined = ModelPositions(st.models[0]);
gemmi::Vec3 moved_centre;
for (const gemmi::Position &p : refined)
moved_centre += p;
moved_centre *= 1.0 / static_cast<double>(refined.size());
double worst = 0;
for (size_t i = 0; i < refined.size(); i++) {
const gemmi::Vec3 expected = result.rotation.multiply(gemmi::Vec3(turned[i]) - centre);
worst = std::max(worst, (gemmi::Vec3(refined[i]) - moved_centre - expected).length());
}
CHECK(worst < 1e-6);
const double trace = result.rotation[0][0] + result.rotation[1][1] + result.rotation[2][2];
CHECK(std::acos(std::clamp((trace - 1) / 2, -1.0, 1.0)) * 180 / PI ==
Catch::Approx(result.angle_deg).margin(1e-6));
}
#ifdef JFJOCH_USE_CUDA
// The gather's bound on the bricks a box falls in along one axis, against every box on small grids -
// wrapped boxes on grids that are not a multiple of the brick included (n = 17: 15, 16, 0 fall in bricks
// 1, 2 and 0). Host arithmetic, no device needed.
TEST_CASE("RigidBodyGPU_AxisBrickBoundCoversWrappedBoxes", "[ModelValidation][gpu]") {
CHECK(RigidBodyGPUEngine::AxisBrickBound(1, 17) == 3);
const int brick = 8;
for (int n = 1; n <= 64; n++)
for (int d = 0; 2 * d + 1 <= n; d++)
for (int c = 0; c < n; c++) {
std::set<int> bricks;
for (int p = c - d; p <= c + d; p++)
bricks.insert(((p % n) + n) % n / brick);
INFO("n " << n << ", d " << d << ", c " << c);
CHECK(bricks.size() <= RigidBodyGPUEngine::AxisBrickBound(d, n));
}
}
// The GPU target is the CPU target's function, computed on the device: at the same placement the two
// give the same residuals and the same Jacobian to rounding (float distances on the device, cuFFT for
// FFTW), on every group of the composition tests, both zones, anisotropic atoms included.
TEST_CASE("RigidBodyGPU_MatchesCPU", "[ModelValidation][gpu]") {
if (get_gpu_count() == 0)
SKIP("No GPU");
Logger logger("RigidBodyGPU_MatchesCPU");
const double q0[6] = {0.12, -0.09, 0.07, 0.15, -0.10, 0.05};
for (const char *cryst : kRigidBodyCrysts) {
const gemmi::Structure st = AnisoCluster(cryst);
const gemmi::SpaceGroup *sg = st.find_spacegroup();
for (double zone : {6.0, 3.5}) {
const auto fobs = OwnAmplitudes(cryst, st, zone, logger);
gemmi::Model cpu_model = st.models[0], gpu_model = st.models[0];
auto pool = TestPool(gpu_model, st.cell, *sg, zone, fobs.v.size(), 1, logger);
RigidBodyTarget cpu(cpu_model, st.cell, *sg, 4);
RigidBodyTargetGPU gpu(*pool, gpu_model, st.cell, *sg, 4);
cpu.SetZone(fobs, zone);
gpu.SetZone(fobs, zone);
const size_t n = cpu.NumObservations();
REQUIRE(gpu.NumObservations() == n);
std::vector<double> rc(n), rg(n), jc(6 * n), jg(6 * n);
REQUIRE(cpu.Residuals(q0, rc.data()));
REQUIRE(gpu.Residuals(q0, rg.data()));
REQUIRE(cpu.Jacobian(q0, jc.data()));
REQUIRE(gpu.Jacobian(q0, jg.data()));
CHECK(gpu.unmatched == cpu.unmatched);
CHECK(gpu.k_sol == cpu.k_sol);
CHECK(gpu.b_sol == cpu.b_sol);
double worst = 0, rms = 0;
for (size_t i = 0; i < n; i++) {
worst = std::max(worst, std::fabs(rc[i] - rg[i]));
rms += rc[i] * rc[i] / static_cast<double>(n);
}
INFO(cryst << " at " << zone << " A: residuals differ by at most " << worst << ", rms residual "
<< std::sqrt(rms));
// Most cases agree to a few 1e-6 of <Fobs>. The bound is set by gemmi's scale fit instead: its
// Levenberg-Marquardt stops at a relative change of 1e-5, so a near-tie in its accept or stop
// decision - which a 1e-5 change of Fcalc can flip, on the CPU alone - moves the scale by up
// to about 1e-4 of |F|.
CHECK(worst <= 2e-4);
for (int j = 0; j < 6; j++) {
double diff = 0, norm = 0;
for (size_t i = 0; i < n; i++) {
diff += gemmi::sq(jc[6 * i + j] - jg[6 * i + j]);
norm += gemmi::sq(jc[6 * i + j]);
}
INFO(cryst << " at " << zone << " A, column " << j << ": relative difference "
<< std::sqrt(diff / norm));
// Measured at up to 5e-5 (the rotation columns, a forward difference of float-gridded
// density; the translation columns agree to about 1e-6). Every column is proportional to
// the overall scale, which the near-tie above can move by 1e-4, so that is the floor of
// an honest bound; 5e-4 leaves room for the difference's own rounding on another card.
CHECK(std::sqrt(diff) <= 5e-4 * std::sqrt(norm) + 1e-9);
}
}
}
}
// The whole Jacobian on the GPU against the central difference of the GPU's own residuals, with the
// bounds of ModelValidation_RigidBodyJacobianMatchesNumericScaleRefit.
TEST_CASE("RigidBodyGPU_JacobianMatchesNumericScaleRefit", "[ModelValidation][gpu]") {
if (get_gpu_count() == 0)
SKIP("No GPU");
const double close[6] = {0.0, 0.0, 0.0, 0.05, -0.04, 0.03};
const double rotated[6] = {0.12, -0.09, 0.07, 0.15, -0.10, 0.05};
for (const char *cryst : kRigidBodyCrysts) {
CheckRigidBodyJacobian(cryst, close, {0, 1, 2, 3, 4, 5}, 0.99, 0.02, true);
CheckRigidBodyJacobian(cryst, rotated, {0, 1, 2}, 0.99, 0.03, true);
}
}
namespace {
// A whole placement on the GPU from the model displaced by 0.54 A and turned 2 deg, against the
// model's own amplitudes; the model is left where the fit put it.
RigidBodyRefineResult DisplacedFit(const char *cryst, gemmi::Structure &st, RigidBodyGPUPool *pool,
Logger &logger) {
const gemmi::SpaceGroup *sg = st.find_spacegroup();
const auto fobs = OwnAmplitudes(cryst, st, 3.0, logger);
std::vector<gemmi::Position> moved = ModelPositions(st.models[0]);
gemmi::Vec3 centre;
for (const gemmi::Position &p : moved)
centre += p;
centre *= 1.0 / static_cast<double>(moved.size());
const double a = 2.0 * PI / 180.0;
const gemmi::Mat33 rz(std::cos(a), -std::sin(a), 0, std::sin(a), std::cos(a), 0, 0, 0, 1);
for (gemmi::Position &p : moved)
p = gemmi::Position(rz.multiply(gemmi::Vec3(p) - centre) + centre + gemmi::Vec3(0.40, -0.30, 0.20));
SetModelPositions(st.models[0], moved);
return RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger, 4, pool);
}
}
// The GPU fit walks where the CPU fit walks: from the same displaced start, to the same placement.
TEST_CASE("RigidBodyGPU_FitAgreesWithCPU", "[ModelValidation][gpu]") {
if (get_gpu_count() == 0)
SKIP("No GPU");
Logger logger("RigidBodyGPU_FitAgreesWithCPU");
for (const char *cryst : {kCryst, kPolarCryst, kRigidBodyCrysts[3]}) {
gemmi::Structure cpu_st = AnisoCluster(cryst), gpu_st = AnisoCluster(cryst);
auto pool = TestPool(gpu_st.models[0], gpu_st.cell, *gpu_st.find_spacegroup(), 3.0, 100000, 1, logger);
const RigidBodyRefineResult cpu = DisplacedFit(cryst, cpu_st, nullptr, logger);
const RigidBodyRefineResult gpu = DisplacedFit(cryst, gpu_st, pool.get(), logger);
CHECK(gpu.converged == cpu.converged);
const std::vector<gemmi::Position> pc = ModelPositions(cpu_st.models[0]), pg = ModelPositions(gpu_st.models[0]);
double rmsd = 0;
for (size_t i = 0; i < pc.size(); i++)
rmsd += pc[i].dist_sq(pg[i]) / static_cast<double>(pc.size());
INFO(cryst << ": CPU " << cpu.angle_deg << " deg " << cpu.shift_A << " A, GPU " << gpu.angle_deg << " deg "
<< gpu.shift_A << " A, " << std::sqrt(rmsd) << " A apart");
// Measured at up to 5e-7 A: the endpoint moves only by what the residuals differ by (a few
// 1e-6 of <Fobs>) over the target's curvature. The bound leaves a factor 20 for another card's
// cuFFT and float rounding, and is still four orders under anything a placement is judged at.
CHECK(std::sqrt(rmsd) < 1e-5);
}
}
// Deterministic: the same fit twice, and on a pool of one engine and of four, gives the same placement
// bit for bit.
TEST_CASE("RigidBodyGPU_Deterministic", "[ModelValidation][gpu]") {
if (get_gpu_count() == 0)
SKIP("No GPU");
Logger logger("RigidBodyGPU_Deterministic");
const char *cryst = kRigidBodyCrysts[4];
std::vector<std::vector<gemmi::Position>> placed;
for (size_t engines : {1, 1, 4}) {
gemmi::Structure st = AnisoCluster(cryst);
auto pool = TestPool(st.models[0], st.cell, *st.find_spacegroup(), 3.0, 100000, engines, logger);
const RigidBodyRefineResult r = DisplacedFit(cryst, st, pool.get(), logger);
CHECK(r.evaluations > 0);
placed.push_back(ModelPositions(st.models[0]));
}
for (size_t k = 1; k < placed.size(); k++)
for (size_t i = 0; i < placed[0].size(); i++) {
CHECK(placed[k][i].x == placed[0][i].x);
CHECK(placed[k][i].y == placed[0][i].y);
CHECK(placed[k][i].z == placed[0][i].z);
}
}
// The model's structure factors on the GPU against the CPU path model validation takes - gemmi's density and
// mask on the symmetrized grid, transformed by FFTW and read off by prepare_asu_data() - in the five groups,
// at a resolution where the mask's shrink step does nothing and at one where it does. The same reflections
// in the same order; F_calc to rounding; F_mask to the few grid points at an atom's mask radius that float
// and double distances put on different sides. The same atoms twice give the same bits.
TEST_CASE("ModelStructureFactorsGPU_MatchesCPU", "[ModelValidation][gpu]") {
if (get_gpu_count() == 0)
SKIP("No GPU");
for (const char *cryst : kRigidBodyCrysts) {
const gemmi::Structure st = AnisoCluster(cryst);
const gemmi::Model &model = st.models[0];
const gemmi::SpaceGroup *sg = st.find_spacegroup();
for (double d_min : {3.5, 1.5}) {
auto dc = ZoneDensity(st, model, d_min);
dc.put_model_density_on_grid(model);
const auto ref_fc = MapToFPhi(dc.grid).prepare_asu_data(d_min, dc.blur, false, false, false);
gemmi::Grid<float> mask;
mask.unit_cell = st.cell;
mask.spacegroup = sg;
mask.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
gemmi::SolventMasker(gemmi::AtomicRadiiSet::Refmac).put_mask_on_grid(mask, model);
const auto ref_fm = MapToFPhi(mask).prepare_asu_data(d_min, 0);
ModelStructureFactorsGPU gpu(0, model, st.cell, *sg, d_min);
REQUIRE(gpu.Supported());
gpu.Reserve();
gemmi::AsuData<std::complex<float>> fc, fm, fc2, fm2;
gpu.Compute(model, fc, fm);
gpu.Compute(model, fc2, fm2);
REQUIRE(fc.v.size() == ref_fc.v.size());
REQUIRE(fm.v.size() == ref_fm.v.size());
double fc_mean = 0, fc_worst = 0, fm_norm = 0, fm_diff = 0;
bool repeat_same = true;
for (size_t i = 0; i < fc.v.size(); i++) {
REQUIRE(fc.v[i].hkl == ref_fc.v[i].hkl);
REQUIRE(fm.v[i].hkl == ref_fm.v[i].hkl);
fc_mean += std::abs(ref_fc.v[i].value) / static_cast<double>(fc.v.size());
fc_worst = std::max(fc_worst, static_cast<double>(std::abs(fc.v[i].value - ref_fc.v[i].value)));
fm_norm += std::norm(ref_fm.v[i].value);
fm_diff += std::norm(fm.v[i].value - ref_fm.v[i].value);
repeat_same = repeat_same && fc2.v[i].value == fc.v[i].value && fm2.v[i].value == fm.v[i].value;
}
INFO(cryst << " at " << d_min << " A: F_calc worst " << fc_worst / fc_mean << " of the mean |F|, F_mask "
<< std::sqrt(fm_diff / fm_norm) << " relative rms");
CHECK(fc_worst <= 1e-4 * fc_mean);
CHECK(std::sqrt(fm_diff) <= 1e-3 * std::sqrt(fm_norm));
CHECK(repeat_same);
}
}
}
// A map from coefficients on the GPU against MapFromFPhi() of gemmi's get_f_phi_on_grid(): the same grid,
// the same values to rounding, the same bits twice.
TEST_CASE("ModelStructureFactorsGPU_MapMatchesCPU", "[ModelValidation][gpu]") {
if (get_gpu_count() == 0)
SKIP("No GPU");
for (const char *cryst : kRigidBodyCrysts) {
const gemmi::Structure st = AnisoCluster(cryst);
const gemmi::Model &model = st.models[0];
const gemmi::SpaceGroup *sg = st.find_spacegroup();
const double d_min = 2.0;
ModelStructureFactorsGPU gpu(0, model, st.cell, *sg, d_min);
REQUIRE(gpu.Supported());
gpu.Reserve();
gemmi::AsuData<std::complex<float>> coef, fm;
gpu.Compute(model, coef, fm);
// Coefficients on every second reflection only, as a map's are on the observed ones.
gemmi::AsuData<std::complex<float>> some = coef;
some.v.clear();
for (size_t i = 0; i < coef.v.size(); i += 2)
some.v.push_back(coef.v[i]);
gemmi::AsuData<std::complex<float>> cpu_coef = some;
cpu_coef.ensure_sorted();
const auto size = gemmi::get_size_for_hkl(cpu_coef, {{0, 0, 0}}, 3.0);
const gemmi::Grid<float> cpu = MapFromFPhi(gemmi::get_f_phi_on_grid<float>(cpu_coef, size, true));
gemmi::AsuData<std::complex<float>> gpu_coef = some;
const gemmi::Grid<float> map = gpu.Map(gpu_coef);
const gemmi::Grid<float> map2 = gpu.Map(gpu_coef);
REQUIRE(map.nu == cpu.nu);
REQUIRE(map.nv == cpu.nv);
REQUIRE(map.nw == cpu.nw);
double rms = 0, worst = 0;
for (size_t i = 0; i < cpu.data.size(); i++) {
rms += gemmi::sq(cpu.data[i]) / static_cast<double>(cpu.data.size());
worst = std::max(worst, static_cast<double>(std::fabs(map.data[i] - cpu.data[i])));
}
INFO(cryst << ": map differs by at most " << worst / std::sqrt(rms) << " of its rms");
CHECK(worst <= 1e-4 * std::sqrt(rms));
CHECK(std::memcmp(map.data.data(), map2.data.data(), map.data.size() * sizeof(float)) == 0);
}
}
#endif
Reference in New Issue View Git Blame Copy Permalink