Jungfraujoch/tests/MolecularReplacementTest.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 <random>

#include <gemmi/dencalc.hpp>
#include <gemmi/it92.hpp>

#include "../image_analysis/structure_refinement/ModelFFT.h"
#include "../image_analysis/structure_refinement/MolecularReplacement.h"
#ifdef JFJOCH_USE_CUDA
#include "../common/CUDAWrapper.h"
#include "../image_analysis/structure_refinement/MolecularReplacementGPU.h"
#endif

namespace {
    // A synthetic body: carbon atoms scattered in a 12 A sphere, centroid at the origin.
    SearchModel synthetic_model() {
        SearchModel m;
        std::mt19937 rng(7);
        std::uniform_real_distribution<double> u(-12.0, 12.0);
        gemmi::Model model(1);
        gemmi::Chain chain("A");
        gemmi::Residue res;
        res.name = "UNK";
        res.seqid = gemmi::SeqId(1, ' ');
        std::vector<gemmi::Vec3> pos;
        while (pos.size() < 300) {
            const gemmi::Vec3 p(u(rng), u(rng), u(rng));
            if (p.length() < 12.0)
                pos.push_back(p);
        }
        gemmi::Vec3 c;
        for (const auto &p : pos)
            c += p;
        c *= 1.0 / pos.size();
        for (const auto &p : pos) {
            gemmi::Atom a;
            a.name = "C";
            a.element = gemmi::El::C;
            a.pos = gemmi::Position(p - c);
            a.b_iso = 15.0f;
            a.occ = 1.0f;
            res.atoms.push_back(a);
            m.radius_A = std::max(m.radius_A, (p - c).length());
        }
        chain.residues.push_back(res);
        model.chains.push_back(chain);
        m.st.models.push_back(model);
        m.residues_read = m.residues_kept = 1;
        m.mass_da = 3600;
        return m;
    }

    gemmi::Mat33 axis_angle(gemmi::Vec3 axis, double deg) {
        axis = axis.normalized();
        const double k = deg * M_PI / 180.0, s = std::sin(k), c1 = 1.0 - std::cos(k);
        const double x = axis.x, y = axis.y, z = axis.z;
        return {1 + c1 * (x * x - 1), -s * z + c1 * x * y, s * y + c1 * x * z,
                s * z + c1 * x * y, 1 + c1 * (y * y - 1), -s * x + c1 * y * z,
                -s * y + c1 * x * z, s * x + c1 * y * z, 1 + c1 * (z * z - 1)};
    }

    // |F|^2 of the placed crystal to 4 A, as merged reflections.
    std::vector<MergedReflection> observed(const gemmi::Structure &crystal) {
        gemmi::DensityCalculator<gemmi::IT92<float>, float> dc;
        dc.d_min = 3.5;
        dc.rate = 1.5;
        dc.set_grid_cell_and_spacegroup(crystal);
        dc.put_model_density_on_grid(crystal.models[0]);
        const auto f = MapToFPhi(dc.grid).prepare_asu_data(3.5, 0, false, false, false);
        std::vector<MergedReflection> out;
        for (const auto &hv : f.v) {
            MergedReflection r;
            r.h = hv.hkl[0];
            r.k = hv.hkl[1];
            r.l = hv.hkl[2];
            r.I = std::norm(hv.value);
            r.sigma = 1;
            r.F = std::abs(hv.value);
            r.d = 1.0 / std::sqrt(crystal.cell.calculate_1_d2(hv.hkl));
            out.push_back(r);
        }
        return out;
    }
}

TEST_CASE("MolecularReplacement_UniformRotationsCoverSO3", "[MolecularReplacement]") {
    const auto rotations = UniformRotations(6.0);
    const double step = 6.0 * M_PI / 180.0;
    // 8 pi^2 is the volume of SO(3) in the measure a 6 deg cube of Euler-like angles samples.
    const double expected = 8 * M_PI * M_PI / (step * step * step);
    CHECK(rotations.size() > 0.8 * expected);
    CHECK(rotations.size() < 1.2 * expected);
    // Every rotation has a sampled one nearby.
    std::mt19937 rng(11);
    std::normal_distribution<double> g;
    const std::vector<gemmi::Mat33> identity{gemmi::Mat33()};
    double worst = 0;
    for (int i = 0; i < 100; i++) {
        const gemmi::Mat33 r = axis_angle({g(rng), g(rng), g(rng)}, std::uniform_real_distribution<double>(0, 180)(rng));
        double best = 180;
        for (const auto &s : rotations)
            best = std::min(best, RotationDistanceDeg(r, s, identity));
        worst = std::max(worst, best);
    }
    CHECK(worst < 6.0);
}

TEST_CASE("MolecularReplacement_RotationAsu", "[MolecularReplacement]") {
    const gemmi::UnitCell cell(40, 40, 60, 90, 90, 90);
    const auto pg = CartesianPointGroup(gemmi::get_spacegroup_by_name("P 41 2 2"), cell);
    REQUIRE(pg.size() == 8);
    const auto all = UniformRotations(10.0);
    const auto asu = RotationAsu(all, pg);
    // One class in eight, give or take the boundaries.
    CHECK(asu.size() > all.size() / 8 * 0.9);
    CHECK(asu.size() < all.size() / 8 * 1.3);
    // Every sampled rotation has a class member in the ASU close by.
    for (size_t i = 0; i < all.size(); i += 37) {
        double best = 180;
        for (const auto &s : asu)
            best = std::min(best, RotationDistanceDeg(s, all[i], pg));
        CHECK(best < 10.0);
    }
}

TEST_CASE("MolecularReplacement_SameRotationGroups", "[MolecularReplacement]") {
    const auto groups = SameRotationGroups(gemmi::get_spacegroup_by_name("P 41"));
    REQUIRE(groups.size() == 4);
    CHECK(groups[0].xhm() == "P 41");
    const auto i = SameRotationGroups(gemmi::get_spacegroup_by_name("I 41 2 2"));
    CHECK(i.size() == 2);   // I 4 2 2, I 41 2 2
}

TEST_CASE("MolecularReplacement_SyntheticPlacement", "[MolecularReplacement]") {
    const SearchModel model = synthetic_model();
    const gemmi::UnitCell cell(42, 50, 58, 90, 90, 90);
    const gemmi::SpaceGroup &sg = gemmi::get_spacegroup_by_name("P 21 21 21");
    const gemmi::Mat33 truth = axis_angle({1, 2, 3}, 50);
    const gemmi::Fractional t(0.17, 0.31, 0.09);
    const gemmi::Structure crystal = PlaceSearchModel(model, truth, t, cell, sg);
    const auto merged = observed(crystal);
    const SearchReflections refl = MakeSearchReflections(merged, cell, sg, 4.0, 12.0);
    REQUIRE(refl.hkl.size() > 300);
    const MolecularTransform ft(model, 4.0, 4);

    // The rotation function's best orientation is the true one, to the sampling.
    const auto pg = CartesianPointGroup(sg, cell);
    const auto rotations = RotationAsu(UniformRotations(6.0), pg);
    const auto rf = RotationFunction(ft, refl, cell, sg, rotations, 4);
    size_t best = 0;
    for (size_t i = 1; i < rf.size(); i++)
        if (rf[i] > rf[best])
            best = i;
    CHECK(RotationDistanceDeg(rotations[best], truth, pg) < 6.0);

    // At the true orientation the translation function peaks at the true position, up to the
    // origins P 21 21 21 allows (a half in each coordinate), and stands far out of its own noise.
    const TranslationPeak peak = TranslationFunction(ft, MakeTranslationSetup(refl, cell, sg), truth);
    CHECK(peak.z > 10.0);
    for (double d : {peak.t.x - t.x, peak.t.y - t.y, peak.t.z - t.z}) {
        const double half = 2 * d - std::round(2 * d);
        CHECK(std::fabs(half) < 0.06);
    }
}


#ifdef JFJOCH_USE_CUDA
TEST_CASE("MolecularReplacement_GPUMatchesCPU", "[MolecularReplacement]") {
    if (get_gpu_count() == 0)
        SKIP("no GPU");
    const SearchModel model = synthetic_model();
    const gemmi::UnitCell cell(42, 50, 58, 90, 90, 90);
    const gemmi::SpaceGroup &sg = gemmi::get_spacegroup_by_name("P 21 21 21");
    const gemmi::Mat33 truth = axis_angle({1, 2, 3}, 50);
    const gemmi::Fractional t(0.17, 0.31, 0.09);
    const auto merged = observed(PlaceSearchModel(model, truth, t, cell, sg));
    const SearchReflections refl = MakeSearchReflections(merged, cell, sg, 4.0, 12.0);
    const MolecularTransform ft(model, 4.0, 4);
    MolecularReplacementGPU gpu(ft);

    const auto rotations = RotationAsu(UniformRotations(15.0), CartesianPointGroup(sg, cell));
    const auto rf_cpu = RotationFunction(ft, refl, cell, sg, rotations, 4);
    const auto rf_gpu = gpu.RotationFunction(refl, cell, sg, rotations);
    REQUIRE(rf_gpu.size() == rf_cpu.size());
    for (size_t i = 0; i < rf_cpu.size(); i++)
        CHECK(std::fabs(rf_gpu[i] - rf_cpu[i]) < 1e-4);

    const TranslationSetup setup = MakeTranslationSetup(refl, cell, sg);
    std::vector<gemmi::Mat33> orientations{truth, rotations[3], rotations[17]};
    const auto tf_gpu = gpu.TranslationFunction(setup, orientations);
    for (size_t i = 0; i < orientations.size(); i++) {
        const TranslationPeak cpu = TranslationFunction(ft, setup, orientations[i]);
        CHECK(std::fabs(tf_gpu[i].z - cpu.z) < 1e-3);
        if (i == 0) {   // the solution's peak is unique; noise peaks can tie between grid points
            CHECK(tf_gpu[i].t.x == Catch::Approx(cpu.t.x));
            CHECK(tf_gpu[i].t.y == Catch::Approx(cpu.t.y));
            CHECK(tf_gpu[i].t.z == Catch::Approx(cpu.t.z));
        }
    }
}
#endif