Jungfraujoch/tests/ModelValidationTest.cpp

// 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 <sstream>

#include <gemmi/mmread_gz.hpp>
#include <gemmi/fourier.hpp>
#include <gemmi/solmask.hpp>

#include "../common/Logger.h"
#include "../rugnux/ModelFFT.h"
#include "../rugnux/ModelGrid.h"
#include "../rugnux/ModelValidation.h"
#include "../rugnux/RigidBodyRefine.h"
#include "../rugnux/SigmaA.h"
#include "../rugnux/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 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);
        }
    }
}

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.
    void CheckRigidBodyJacobian(const char *cryst, const double q0[6], const std::vector<int> &columns,
                                double min_cosine, double max_norm_error) {
        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];
            RigidBodyTarget target(model, st.cell, *sg, 4);
            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);
    }
}