mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions 4 Packages Projects Releases 81 Wiki Activity
Files
babff6d8bff8f0f039e510f6b12b09ca381df4ea
Jungfraujoch/image_analysis/structure_refinement/MolecularReplacement.cpp
T
leonarski_fandClaude Opus 5.5 babff6d8bf rugnux --af-model: the null keeps the copies' arrangement 
The null a placement must beat turned the whole assembly about its common
centroid, so with several copies it also moved them off the arrangement the
search found. Under tNCS, two parallel copies a tNCS vector apart fit the data
whatever the model is, and the null threw that fit away: wrong deposited models
reached TF Z 14-47 and 8-34 sd over the null on tNCS crystals, and two of them
decided the space group (SPACE_GROUP, a reflection file in the wrong group).

Each copy is now turned about its own centroid, all by the same rotation, so
only the model's own orientation is random. On a control of 38 open-arm sets
(deposited chain randomly re-oriented as search model, 3 unrelated deposited
models and a fold-scrambled decoy per set) the wrong models' null falls to
1.4-4.2 sd on the tNCS sets (one at 7.2, left UNDECIDED), both wrong decisions
become NOT_PLACED, and no right-model or AlphaFold-model decision changes
(nulls stay 15-107 sd). The 193-set AF survey replay loses no GREEN (two
UNDECIDED become NOT_PLACED); the four twin/tNCS targets stay correct.

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

2708 lines
135 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 "MolecularReplacement.h"
#include <algorithm>
#include <array>
#include <atomic>
#include <chrono>
#include <cmath>
#include <future>
#include <functional>
#include <numeric>
#include <optional>
#include <map>
#include <mutex>
#include <tuple>
#include <random>
#include <unordered_map>
#include <fstream>
#include <nlohmann/json.hpp>
#include <gemmi/mmread_gz.hpp> // read_structure_gz
#include <gemmi/modify.hpp> // remove_hydrogens
#include <gemmi/dencalc.hpp> // DensityCalculator
#include <gemmi/it92.hpp> // IT92 x-ray form factors
#include <gemmi/scaling.hpp> // Scaling (bulk solvent + anisotropic B)
#include "../../common/JFJochMath.h" // PI
#include "../../common/Logger.h"
#include "../../common/ParallelFor.h"
#include "ModelFFT.h"
#include "ModelGrid.h"
#include "ModelScaling.h"
#include "ModelValidation.h"
#include "RigidBodyRefine.h"
#ifdef JFJOCH_USE_CUDA
#include "MolecularReplacementGPU.h"
#include "RigidBodyGPU.h"
#include "RigidBodyGPUEngine.h"
#include "ModelStructureFactorsGPU.h"
#include "../../common/CUDAWrapper.h"
#include "../../common/JFJochException.h"
#endif
namespace {
using Table = gemmi::IT92<float>;
// The search's resolution range: the low-resolution end is where the bulk solvent, which the search
// does not model, stops dominating; the high-resolution end is where a predicted model still agrees
// with the crystal's structure in its details.
constexpr double SEARCH_D_MIN = 4.0;
constexpr double SEARCH_D_MAX = 12.0;
constexpr int SEARCH_BINS = 10;
constexpr double ROTATION_STEP_DEG = 6.0;
// Orientations the translation function is run for, per merge, and placements per space group taken
// through the rigid body.
constexpr int TRANSLATION_ROTATIONS = 50;
constexpr int RESCORED_PER_GROUP = 3;
// Hill-climbing steps per step size of a rotation-function peak's refinement.
constexpr int PEAK_REFINE_ITERATIONS = 4;
// The resolution the placements are refined and scored to: whether a model sits in the right place
// is a low-resolution question, the rigid body's ladder ends at 3.5 A anyway, and every R-factor of
// the comparison is made to the same limit.
constexpr double SCORE_D_MIN = 3.0;
// The likelihood a placement is judged by is summed over the reflections to LLG_D_MIN - finer than the
// search, where the information a small cell's few reflections hold is - or over as many as LLG_TERMS.
constexpr double LLG_D_MIN = 2.5;
constexpr size_t LLG_TERMS = 20000;
// The decision. A translation-function peak this many standard deviations above the function's mean
// is a solution, where the model packs and one copy fits the cell; the McCoy et al. (2007) J. Appl.
// Cryst. 40, 658-674 convention for a single-copy search, taken over to this Z-score's own null on
// the Crowther & Blow function (measured: the highest peak of a wrong orientation reads 4-6).
constexpr double TF_Z_PLACED = 8.0;
constexpr double CLASH_MAX = 0.05;
// Matthews (1968) J. Mol. Biol. 33, 491-497: protein crystals have Vm of about 1.7 to 3.5 A^3/Da.
// Below VM_MIN one copy does not fit the cell.
constexpr double VM_MIN = 1.5;
// Copies are searched for one after another up to the largest copy number at least this probable
// (MatthewsProbabilities); a copy is kept where its translation function stands out beside the copies
// already placed, it packs, and the R-free of all of them drops by at least COPY_R_FREE_GAIN.
constexpr double MATTHEWS_P_MIN = 0.01;
constexpr double COPY_R_FREE_GAIN = 0.01;
// The translation-function peaks of a further copy its assembly correlation is computed for.
constexpr size_t COPY_TRIALS = 10;
constexpr double DALTON_PER_RESIDUE = 110.0;
// How far below every other candidate of its merge the chosen one's R-free must be.
constexpr double R_FREE_MARGIN = 0.03;
// The null (as in model validation): the same model refitted from random orientations about its
// centroid, at a coarse limit, and how far above that the real placement's R-work has to sit.
constexpr int NULL_REPLICATES = 9;
constexpr unsigned NULL_SEED = 20261011;
constexpr double NULL_D_MIN = 3.5;
// 3 sd of the null, one-sided, allowing for its mean and sd being estimated from NULL_REPLICATES draws:
// a new draw lies t(n-1) sqrt(1 + 1/n) sample sds above their mean at that probability, 4.28 x 1.054 for
// n = 9 - the gate model validation's fit is read at (ModelValidation.cpp, MODEL_FIT_SIGMA).
constexpr double NULL_SIGMA_MIN = 4.51;
double seconds_since(std::chrono::steady_clock::time_point t0) {
return std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count();
}
// A rotation drawn uniformly from SO(3), through a uniform random unit quaternion.
// Following Shoemake (1992) Graphics Gems III, 124-132
gemmi::Mat33 random_rotation(std::mt19937 &rng) {
std::uniform_real_distribution<double> u(0.0, 1.0);
const double u1 = u(rng), t2 = 2 * PI * u(rng), t3 = 2 * PI * u(rng);
const double r1 = std::sqrt(1 - u1), r2 = std::sqrt(u1);
const double x = r1 * std::sin(t2), y = r1 * std::cos(t2), z = r2 * std::sin(t3), w = r2 * std::cos(t3);
return {1 - 2 * (y * y + z * z), 2 * (x * y - z * w), 2 * (x * z + y * w),
2 * (x * y + z * w), 1 - 2 * (x * x + z * z), 2 * (y * z - x * w),
2 * (x * z - y * w), 2 * (y * z + x * w), 1 - 2 * (x * x + y * y)};
}
// A rotation by `deg` about Cartesian axis 0, 1 or 2.
gemmi::Mat33 axis_rotation(int axis, double deg) {
const double c = std::cos(deg * PI / 180.0), s = std::sin(deg * PI / 180.0);
if (axis == 0)
return {1, 0, 0, 0, c, -s, 0, s, c};
if (axis == 1)
return {c, 0, s, 0, 1, 0, -s, 0, c};
return {c, -s, 0, s, c, 0, 0, 0, 1};
}
double trace(const gemmi::Mat33 &m) {
return m.a[0][0] + m.a[1][1] + m.a[2][2];
}
gemmi::Vec3 miller_vec(const std::array<int, 3> &h) {
return {static_cast<double>(h[0]), static_cast<double>(h[1]), static_cast<double>(h[2])};
}
// h R for a symmetry operator's rotation: the index the operator takes h to.
gemmi::Miller rotated_index(const gemmi::Op &op, const gemmi::Miller &h) {
const auto k = op.apply_to_hkl_without_division(h);
return {{k[0] / gemmi::Op::DEN, k[1] / gemmi::Op::DEN, k[2] / gemmi::Op::DEN}};
}
long hkl_key(const gemmi::Miller &h) {
return (static_cast<long>(h[0]) + 512) * 1048576L + (static_cast<long>(h[1]) + 512) * 1024L + (h[2] + 512);
}
int good_even_size(int n) {
while (n % 2 != 0 || !gemmi::has_small_factorization(n))
++n;
return n;
}
// --- placements scored against the data: R-work / R-free after the bulk solvent and the scale ---
struct PlacementScore {
double r_work = 1, r_free = 1;
double k_sol = NAN, b_sol = NAN;
// Against the data as twinned by a law (score_placement's `law`): at the fraction that gives the lowest
// R-work.
double r_work_twin = NAN, r_free_twin = NAN, twin_fraction = NAN;
};
gemmi::AsuData<gemmi::ValueSigma<float>> working_amplitudes(const std::vector<MergedReflection> &merged,
const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
double d_min) {
const gemmi::GroupOps gops = sg.operations();
const gemmi::ReciprocalAsu asu(&sg);
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
fobs.unit_cell_ = cell;
fobs.spacegroup_ = &sg;
for (const MergedReflection &r : merged) {
if (std::isnan(r.F) || r.rfree_flag || r.d < d_min)
continue;
gemmi::Miller h{{r.h, r.k, r.l}};
if (!asu.is_in(h))
h = asu.to_asu(h, gops).first;
fobs.v.push_back({h, {r.F, 1.0f}});
}
fobs.ensure_asu();
fobs.ensure_sorted();
return fobs;
}
// The model's Fcalc and bulk-solvent mask to d_min on the CPU, as model validation makes them.
void model_factors(const gemmi::Structure &st, double d_min, gemmi::AsuData<std::complex<float>> &fmodel,
gemmi::AsuData<std::complex<float>> &fmask, size_t nthreads) {
gemmi::DensityCalculator<Table, float> dc;
dc.d_min = d_min;
dc.rate = 1.5;
dc.set_grid_cell_and_spacegroup(st);
dc.set_refmac_compatible_blur(st.models[0]);
gemmi::Grid<float> mask;
mask.unit_cell = dc.grid.unit_cell;
mask.spacegroup = dc.grid.spacegroup;
mask.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
const std::vector<size_t> leaders = OrbitLeaders(mask, nthreads);
PutMaskOnGrid(mask, st.models[0], leaders, nthreads);
fmask = MapToFPhi(mask).prepare_asu_data(d_min, 0);
PutModelDensityOnGrid(dc, st.models[0], leaders, nthreads);
fmodel = MapToFPhi(dc.grid).prepare_asu_data(d_min, dc.blur, false, false, false);
}
// R-work and R-free of the scaled Fcalc over the reflections to d_min.
PlacementScore r_factors(const gemmi::AsuData<std::complex<float>> &fmodel, const gemmi::Scaling<float> &scaling,
const gemmi::SpaceGroup &sg, const std::vector<MergedReflection> &merged, double d_min) {
const gemmi::GroupOps gops = sg.operations();
const gemmi::ReciprocalAsu asu(&sg);
std::unordered_map<long, std::pair<double, bool>> obs;
for (const MergedReflection &r : merged) {
if (std::isnan(r.F) || r.d < d_min)
continue;
gemmi::Miller h{{r.h, r.k, r.l}};
if (!asu.is_in(h))
h = asu.to_asu(h, gops).first;
obs[hkl_key(h)] = {r.F, r.rfree_flag};
}
double num_w = 0, den_w = 0, num_f = 0, den_f = 0;
for (const auto &hv : fmodel.v) {
const auto it = obs.find(hkl_key(hv.hkl));
if (it == obs.end())
continue;
const double fo = it->second.first, fc = std::abs(hv.value);
if (it->second.second) {
num_f += std::fabs(fo - fc);
den_f += fo;
} else {
num_w += std::fabs(fo - fc);
den_w += fo;
}
}
PlacementScore out;
out.k_sol = scaling.k_sol;
out.b_sol = scaling.b_sol;
out.r_work = den_w > 0 ? num_w / den_w : 1;
out.r_free = den_f > 0 ? num_f / den_f : 1;
return out;
}
// Model validation's structure-factor engines on the GPU, one per group, copy count and resolution, made
// where the first placement of that kind is scored - or none, where there is no GPU, the cell is too small
// for its gridding or one would take over a quarter of the card. Decided on the card's total memory, so
// the same input always takes the same path.
#ifdef JFJOCH_USE_CUDA
class ScoreEngines {
public:
ModelStructureFactorsGPU *get(const gemmi::Structure &st, const gemmi::SpaceGroup &sg, double d_min) {
if (get_gpu_count() == 0)
return nullptr;
size_t atoms = 0;
for (const gemmi::Chain &ch : st.models[0].chains)
for (const gemmi::Residue &r : ch.residues)
atoms += r.atoms.size();
std::lock_guard lock(m_);
auto &e = engines_[{sg.xhm(), atoms, d_min}];
if (!e.tried) {
e.tried = true;
try {
auto sf = std::make_unique<ModelStructureFactorsGPU>(RigidBodyGPUEngine::CurrentDevice(), st.models[0],
st.cell, sg, d_min);
if (sf->Supported() && sf->DeviceBytes() <= sf->DeviceTotalMemory() / 4) {
sf->Reserve();
e.sf = std::move(sf);
}
} catch (const JFJochException &) {
cuda_clear_error();
}
}
return e.sf.get();
}
private:
struct Entry {
bool tried = false;
std::unique_ptr<ModelStructureFactorsGPU> sf;
};
std::map<std::tuple<std::string, size_t, double>, Entry> engines_;
std::mutex m_;
};
#else
class ScoreEngines {};
#endif
// The model's Fcalc and bulk-solvent mask to d_min - on the GPU where `engines` has one - scaled to the
// working amplitudes as model validation scales them (ModelValidation.cpp, fit_model), and the R-factors
// over the reflections to d_min. `solvent`: k_sol and b_sol held at its, not fitted.
// R-work and R-free against the data as twinned by `law`: |F_twin(h)|^2 = (1-a)|F(h)|^2 + a|F(law h)|^2 of the
// scaled Fcalc, one overall scale refitted at each fraction a, a scanned from 0 to 0.5 in steps of 0.01 and
// chosen on the working set.
// Following Yeates (1997) Methods Enzymol. 276, 344-358
void twin_r_factors(const gemmi::AsuData<std::complex<float>> &fmodel, const gemmi::SpaceGroup &sg,
const std::vector<MergedReflection> &merged, double d_min, const gemmi::Op &law,
PlacementScore &out) {
const gemmi::GroupOps gops = sg.operations();
const gemmi::ReciprocalAsu asu(&sg);
std::unordered_map<long, double> model_i;
for (const auto &hv : fmodel.v)
model_i[hkl_key(hv.hkl)] = std::norm(hv.value);
struct Term { double fo, i1, i2; bool free; };
std::vector<Term> terms;
for (const MergedReflection &r : merged) {
if (std::isnan(r.F) || r.d < d_min)
continue;
gemmi::Miller h{{r.h, r.k, r.l}};
if (!asu.is_in(h))
h = asu.to_asu(h, gops).first;
const auto a = model_i.find(hkl_key(h));
const auto b = model_i.find(hkl_key(asu.to_asu(law.apply_to_hkl(h), gops).first));
if (a != model_i.end() && b != model_i.end())
terms.push_back({r.F, a->second, b->second, r.rfree_flag});
}
double best = INFINITY;
for (int step = 0; step <= 50; step++) {
const double a = 0.01 * step;
double fo_fc = 0, fc_fc = 0;
for (const Term &t : terms)
if (!t.free) {
const double fc = std::sqrt((1 - a) * t.i1 + a * t.i2);
fo_fc += t.fo * fc;
fc_fc += fc * fc;
}
const double k = fc_fc > 0 ? fo_fc / fc_fc : 1.0;
double num_w = 0, den_w = 0, num_f = 0, den_f = 0;
for (const Term &t : terms) {
const double d = std::fabs(t.fo - k * std::sqrt((1 - a) * t.i1 + a * t.i2));
(t.free ? num_f : num_w) += d;
(t.free ? den_f : den_w) += t.fo;
}
if (den_w > 0 && den_f > 0 && num_w / den_w < best) {
best = num_w / den_w;
out.r_work_twin = best;
out.r_free_twin = num_f / den_f;
out.twin_fraction = a;
}
}
}
PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGroup &sg,
const std::vector<MergedReflection> &merged, double d_min, size_t nthreads,
ScoreEngines *engines = nullptr, const PlacementScore *solvent = nullptr,
const gemmi::Op *law = nullptr) {
gemmi::AsuData<std::complex<float>> fmask, fmodel;
#ifdef JFJOCH_USE_CUDA
if (ModelStructureFactorsGPU *gpu = engines != nullptr ? engines->get(st, sg, d_min) : nullptr)
gpu->Compute(st.models[0], fmodel, fmask);
else
#endif
model_factors(st, d_min, fmodel, fmask, nthreads);
gemmi::AsuData<gemmi::ValueSigma<float>> fobs = working_amplitudes(merged, st.cell, sg, d_min);
gemmi::Scaling<float> scaling(st.cell, &sg);
scaling.use_solvent = true;
scaling.prepare_points(fmodel, fobs, &fmask);
if (solvent != nullptr) {
// The bulk solvent of an earlier fit of the same crystal, held, as the rigid body holds it.
scaling.k_sol = solvent->k_sol;
scaling.b_sol = solvent->b_sol;
scaling.fix_k_sol = true;
scaling.fix_b_sol = true;
scaling.fit_isotropic_b_approximately();
scaling.fit_parameters();
} else {
FitModelScale(scaling, {.b_hi = SolventBMax(st.models[0])}, nthreads);
}
scaling.scale_data(fmodel, &fmask);
PlacementScore out = r_factors(fmodel, scaling, sg, merged, d_min);
if (law != nullptr)
twin_r_factors(fmodel, sg, merged, d_min, *law, out);
return out;
}
// The rigid-body engines for one model in one group, where there is a GPU.
struct RigidBodyEngines {
#ifdef JFJOCH_USE_CUDA
std::unique_ptr<RigidBodyGPUPool> pool;
#endif
RigidBodyGPUPool *get() {
#ifdef JFJOCH_USE_CUDA
return pool.get();
#else
return nullptr;
#endif
}
};
RigidBodyEngines make_engines(const gemmi::Structure &st, const gemmi::SpaceGroup &sg,
const std::vector<MergedReflection> &merged, double d_min, size_t engines,
Logger &logger) {
RigidBodyEngines out;
#ifdef JFJOCH_USE_CUDA
const double finest = RigidBodyLadder(d_min).back();
size_t observations = 0;
for (const MergedReflection &r : merged)
if (!std::isnan(r.F) && r.d >= finest)
++observations;
out.pool = RigidBodyGPUPool::Create(st.models[0], st.cell, sg, d_min, observations, engines, logger);
#else
(void) st; (void) sg; (void) merged; (void) d_min; (void) engines; (void) logger;
#endif
return out;
}
// The placement through the rigid body, and its R-factors where the rigid body left it. A search
// placement is a few degrees and a grid step from its answer, well inside the rigid body's reach, so
// the step is taken without the R-free check model validation makes before committing one.
// `moved`: where given, the rotation the rigid body turned the model by about its centroid. `fixed`: the
// other copies, held where they are.
void refine_rigid_body(gemmi::Structure &st, const gemmi::SpaceGroup &sg, const std::vector<MergedReflection> &merged,
double d_min, RigidBodyGPUPool *pool, size_t nthreads, gemmi::Mat33 *moved = nullptr,
const gemmi::Model *fixed = nullptr) {
const std::vector<gemmi::Position> start = ModelPositions(st.models[0]);
Logger quiet = Logger::Buffered();
const RigidBodyRefineResult rb = RefineRigidBody(st.models[0], st.cell, sg,
working_amplitudes(merged, st.cell, sg, d_min), d_min, quiet,
nthreads, pool, fixed);
if (!rb.converged)
SetModelPositions(st.models[0], start);
if (moved != nullptr)
*moved = rb.converged ? rb.rotation : gemmi::Mat33();
}
PlacementScore refine_placement(gemmi::Structure &st, const gemmi::SpaceGroup &sg,
const std::vector<MergedReflection> &merged, double d_min, RigidBodyGPUPool *pool,
size_t nthreads, ScoreEngines *engines, gemmi::Mat33 *moved = nullptr) {
refine_rigid_body(st, sg, merged, d_min, pool, nthreads, moved);
return score_placement(st, sg, merged, d_min, nthreads, engines);
}
// The copies of `assembly` (its chains split evenly, in order) as a crystal of group `sg`, each refined beside
// the others, and scored.
PlacementScore carry_copies(const gemmi::Structure &assembly, int copies, const gemmi::SpaceGroup &sg,
const std::vector<MergedReflection> &merged, double d_min, size_t nthreads,
ScoreEngines *scorers, Logger &logger) {
std::vector<gemmi::Structure> body = SplitCopies(assembly, copies);
for (gemmi::Structure &b : body) {
b.spacegroup_hm = sg.xhm();
b.setup_cell_images();
}
RigidBodyEngines engines = make_engines(body[0], sg, merged, d_min, 1, logger);
for (int k = 0; k < copies; k++) {
std::vector<gemmi::Structure> rest;
for (int o = 0; o < copies; o++)
if (o != k)
rest.push_back(body[o]);
const gemmi::Structure fixed = AssembleCopies(rest);
refine_rigid_body(body[k], sg, merged, d_min, engines.get(), nthreads, nullptr, &fixed.models[0]);
}
return score_placement(AssembleCopies(body), sg, merged, d_min, nthreads, scorers);
}
// The fractional position of the model's centroid.
gemmi::Fractional centroid(const gemmi::Structure &st) {
const std::vector<gemmi::Position> pos = ModelPositions(st.models[0]);
gemmi::Vec3 c;
for (const gemmi::Position &p : pos)
c += p;
return st.cell.fractionalize(gemmi::Position(c * (1.0 / static_cast<double>(pos.size()))));
}
// The placed model's R-work at NULL_D_MIN against the same model refitted from random orientations
// about its centroid, in standard deviations of the latter: the null model validation builds for a
// model that claims something (ModelValidation.cpp). An orientation the crystal cannot tell from the
// placement's - within the rigid body's reach of one of its symmetry classes - is drawn again.
// With several copies each is turned about its own centroid, all by the same rotation: the copies keep
// their arrangement, and with it the fit that tNCS gives any model placed parallel at its vector, so that
// only the model's own orientation is random.
// `translate`: where given, each replicate is placed where the translation function puts the model in its
// orientation (the centroids, fractional, of the model in each of the orientations given), as the search
// placed the real one, instead of at the real placement's centroid; the search chose that position among
// all others for the real orientation, and a null that is not given the same choice reads high.
double null_sigma(const gemmi::Structure &placed, int copies, const gemmi::Mat33 &orientation, const gemmi::SpaceGroup &sg,
const std::vector<MergedReflection> &merged, double d_min, RigidBodyGPUPool *pool,
size_t nthreads, ScoreEngines *scorers,
const std::function<std::vector<gemmi::Fractional>(const std::vector<gemmi::Mat33> &)> &translate) {
const double null_d_min = std::max(d_min, NULL_D_MIN);
std::vector<MergedReflection> coarse;
for (const MergedReflection &r : merged)
if (r.d >= null_d_min)
coarse.push_back(r);
const double real = score_placement(placed, sg, coarse, null_d_min, nthreads, scorers).r_work;
const std::vector<gemmi::Mat33> pg = CartesianPointGroup(sg, placed.cell);
const double reach = RigidBodyReachDeg(placed.models[0], null_d_min);
std::vector<gemmi::Mat33> draws;
std::mt19937 rng(NULL_SEED);
while (static_cast<int>(draws.size()) < NULL_REPLICATES) {
const gemmi::Mat33 r = random_rotation(rng);
if (RotationDistanceDeg(r.multiply(orientation), orientation, pg) >= reach)
draws.push_back(r);
}
std::vector<gemmi::Fractional> where;
if (translate) {
std::vector<gemmi::Mat33> turned;
for (const gemmi::Mat33 &r : draws)
turned.push_back(r.multiply(orientation));
where = translate(turned);
}
std::vector<double> r_work(NULL_REPLICATES);
std::atomic<int> next{0};
std::vector<std::future<void>> workers;
for (size_t t = 0; t < std::min<size_t>(std::max<size_t>(nthreads, 1), NULL_REPLICATES); t++)
workers.push_back(std::async(std::launch::async, [&] {
for (int i = next++; i < NULL_REPLICATES; i = next++) {
gemmi::Structure rep = placed;
std::vector<gemmi::Chain> &chains = rep.models[0].chains;
for (int k = 0; k < copies; k++) {
// Copy k: the chains SplitCopies gives it.
std::vector<gemmi::Atom *> atoms;
for (size_t ch = 0; ch < chains.size(); ch++)
if (static_cast<int>(ch * copies / chains.size()) == k)
for (gemmi::Residue &r : chains[ch].residues)
for (gemmi::Atom &a : r.atoms)
atoms.push_back(&a);
gemmi::Vec3 c;
for (const gemmi::Atom *a : atoms)
c += a->pos;
c *= 1.0 / static_cast<double>(atoms.size());
const gemmi::Vec3 to = translate ? gemmi::Vec3(placed.cell.orthogonalize(where[i])) : c;
for (gemmi::Atom *a : atoms)
a->pos = gemmi::Position(draws[i].multiply(gemmi::Vec3(a->pos) - c) + to);
}
r_work[i] = refine_placement(rep, sg, coarse, null_d_min, pool, nthreads, scorers).r_work;
}
}));
for (auto &w : workers)
w.get();
const double mean = std::accumulate(r_work.begin(), r_work.end(), 0.0) / NULL_REPLICATES;
double s2 = 0;
for (double x : r_work)
s2 += (x - mean) * (x - mean);
const double sd = std::sqrt(s2 / (NULL_REPLICATES - 1));
return sd > 0 ? (mean - real) / sd : 0.0;
}
bool same_rotations_and_centring(const gemmi::GroupOps &a, const gemmi::GroupOps &b) {
if (a.sym_ops.size() != b.sym_ops.size() || a.cen_ops.size() != b.cen_ops.size())
return false;
for (const gemmi::Op &op : a.sym_ops)
if (b.find_by_rotation(op.rot) == nullptr)
return false;
for (const gemmi::Op::Tran &c : a.cen_ops)
if (std::find(b.cen_ops.begin(), b.cen_ops.end(), c) == b.cen_ops.end())
return false;
return true;
}
// Whether b is a on a shifted origin, the shift on the twelfths every crystallographic one lies on.
bool same_up_to_origin(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b) {
if (a.number != b.number)
return false;
const gemmi::GroupOps bops = b.operations();
for (int i = 0; i < 12; i++)
for (int j = 0; j < 12; j++)
for (int k = 0; k < 12; k++) {
gemmi::Op shift = gemmi::Op::identity();
shift.tran = {i * gemmi::Op::DEN / 12, j * gemmi::Op::DEN / 12, k * gemmi::Op::DEN / 12};
gemmi::GroupOps aops = a.operations();
aops.change_basis_forward(shift);
if (aops.is_same_as(bops))
return true;
}
return false;
}
bool enantiomorphs(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b) {
if (!a.is_enantiomorphic())
return false;
gemmi::GroupOps e = a.operations();
e.change_basis_forward(a.change_of_hand_op());
const gemmi::SpaceGroup *other = gemmi::find_spacegroup_by_ops(e);
return other != nullptr && other->number == b.number;
}
} // namespace
// --- the model ---
namespace {
// The atoms' centroid moved to the origin, and the radius, mass and coordinate error of what is left.
void finish_search_model(SearchModel &m) {
std::vector<gemmi::Position> pos = ModelPositions(m.st.models[0]);
double error2 = 0;
for (const gemmi::Chain &ch : m.st.models[0].chains)
for (const gemmi::Residue &r : ch.residues)
for (const gemmi::Atom &a : r.atoms)
error2 += 3.0 * a.b_iso / (8.0 * PI * PI); // the B from pLDDT, back to its rmsd^2
m.error_A = std::sqrt(error2 / static_cast<double>(pos.size()));
gemmi::Vec3 c;
for (const gemmi::Position &p : pos)
c += p;
c *= 1.0 / static_cast<double>(pos.size());
m.radius_A = 0;
for (gemmi::Position &p : pos) {
p = gemmi::Position(gemmi::Vec3(p) - c);
m.radius_A = std::max(m.radius_A, p.length());
}
SetModelPositions(m.st.models[0], pos);
m.mass_da = DALTON_PER_RESIDUE * m.residues_kept;
}
} // namespace
SearchModel PrepareSearchModel(const std::string &path) {
SearchModel out;
try {
out.st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
} catch (const std::exception &e) {
out.failure_reason = std::string("cannot read the model: ") + e.what();
return out;
}
if (out.st.models.empty()) {
out.failure_reason = "the model has no atoms";
return out;
}
out.st.models.resize(1);
gemmi::remove_hydrogens(out.st);
gemmi::Model &model = out.st.models[0];
for (gemmi::Chain &ch : model.chains) {
std::vector<gemmi::Residue> kept;
for (gemmi::Residue &r : ch.residues) {
for (const gemmi::Atom &a : r.atoms)
if (a.b_iso > 100.0f) {
out.failure_reason = "the B-factor column holds values above 100, so it is not a pLDDT: "
"not a predicted model";
return out;
}
// The residues a prediction is confident of: pLDDT 70 and above, as Oeffner et al. (2022)
// Acta Cryst. D78, 1303-1314 trim an AlphaFold model for molecular replacement. Each kept with
// its position in the whole prediction, which its predicted aligned error is indexed by.
if (!r.atoms.empty() && r.atoms[0].b_iso >= 70.0f) {
kept.push_back(std::move(r));
out.residue_index.push_back(out.residues_read);
}
out.residues_read++;
}
ch.residues = std::move(kept);
out.residues_kept += static_cast<int>(ch.residues.size());
}
std::erase_if(model.chains, [](const gemmi::Chain &ch) { return ch.residues.empty(); });
if (out.residues_kept == 0) {
out.failure_reason = "no residue of the model has a pLDDT of 70 or more";
return out;
}
// B from pLDDT through the coordinate error it predicts: rmsd = 1.5 exp(4 (0.7 - pLDDT/100)) A
// and B = 8 pi^2 rmsd^2 / 3. Following Oeffner et al. (2022) Acta Cryst. D78, 1303-1314
for (gemmi::Chain &ch : model.chains)
for (gemmi::Residue &r : ch.residues)
for (gemmi::Atom &a : r.atoms) {
const double rmsd = 1.5 * std::exp(4.0 * (0.7 - a.b_iso / 100.0));
a.b_iso = static_cast<float>(8.0 * PI * PI * rmsd * rmsd / 3.0);
}
finish_search_model(out);
return out;
}
// --- domains from the predicted aligned error ---
std::string PaePathFor(const std::string &model_path) {
const size_t slash = model_path.find_last_of("/\\");
const std::string dir = slash == std::string::npos ? std::string() : model_path.substr(0, slash + 1);
const std::string name = model_path.substr(dir.size());
const size_t at = name.rfind("-model_v");
if (at == std::string::npos)
return {};
const size_t dot = name.find('.', at);
const std::string version = name.substr(at + 8, dot == std::string::npos ? std::string::npos : dot - at - 8);
return dir + name.substr(0, at) + "-predicted_aligned_error_v" + version + ".json";
}
// A prediction cut into domains by its predicted aligned error before molecular replacement, as Oeffner et
// al. (2022) Acta Cryst. D78, 1303-1314 do; here by average linkage.
std::vector<std::vector<int>> PaeDomains(const std::string &pae_path, const SearchModel &model) {
// Two domains a prediction does not place relative to each other with confidence are this far apart in
// mean PAE (within one domain it is a few angstroms); a domain is at least this many residues.
constexpr double PAE_DOMAIN_A = 10.0;
constexpr int DOMAIN_MIN_RESIDUES = 40;
// The clustering starts from runs of this many consecutive residues kept.
constexpr int PIECE = 5;
std::vector<std::vector<double>> pae;
try {
std::ifstream in(pae_path);
if (!in)
return {};
const nlohmann::json j = nlohmann::json::parse(in);
const nlohmann::json &m = j.is_array() ? j.at(0).at("predicted_aligned_error") : j.at("predicted_aligned_error");
for (const auto &row : m)
pae.push_back(row.get<std::vector<double>>());
} catch (const std::exception &) {
return {};
}
const size_t n = pae.size();
if (n != static_cast<size_t>(model.residues_read)
|| std::any_of(pae.begin(), pae.end(), [n](const std::vector<double> &r) { return r.size() != n; }))
return {};
// The pieces: runs of consecutive kept residues, at most PIECE long.
std::vector<std::vector<int>> members;
for (size_t i = 0; i < model.residue_index.size(); i++) {
const int r = model.residue_index[i];
if (members.empty() || static_cast<int>(members.back().size()) == PIECE || members.back().back() != r - 1)
members.push_back({});
members.back().push_back(r);
}
// Average linkage on the symmetrised PAE: the two clusters with the lowest mean PAE between them merged
// until that is PAE_DOMAIN_A or more.
const size_t m = members.size();
std::vector<double> sum(m * m, 0.0);
for (size_t a = 0; a < m; a++)
for (size_t b = 0; b < m; b++)
for (int i : members[a])
for (int k : members[b])
sum[a * m + b] += (pae[i][k] + pae[k][i]) / 2.0;
std::vector<char> active(m, 1);
while (true) {
double best = PAE_DOMAIN_A;
size_t ba = 0, bb = 0;
for (size_t a = 0; a < m; a++)
for (size_t b = a + 1; b < m; b++)
if (active[a] && active[b]) {
const double mean = sum[a * m + b] / static_cast<double>(members[a].size() * members[b].size());
if (mean < best) {
best = mean;
ba = a;
bb = b;
}
}
if (!(best < PAE_DOMAIN_A))
break;
for (size_t c = 0; c < m; c++) {
sum[ba * m + c] += sum[bb * m + c];
sum[c * m + ba] += sum[c * m + bb];
}
members[ba].insert(members[ba].end(), members[bb].begin(), members[bb].end());
active[bb] = 0;
}
std::vector<std::vector<int>> out;
for (size_t a = 0; a < m; a++)
if (active[a] && static_cast<int>(members[a].size()) >= DOMAIN_MIN_RESIDUES) {
std::sort(members[a].begin(), members[a].end());
out.push_back(members[a]);
}
std::stable_sort(out.begin(), out.end(), [](const std::vector<int> &a, const std::vector<int> &b) {
return a.size() > b.size();
});
return out;
}
SearchModel DomainModel(const SearchModel &model, const std::vector<int> &residues) {
SearchModel out = model;
out.residue_index.clear();
out.residues_kept = 0;
size_t i = 0;
for (gemmi::Chain &ch : out.st.models[0].chains) {
std::vector<gemmi::Residue> kept;
for (gemmi::Residue &r : ch.residues) {
const int index = model.residue_index[i++];
if (std::binary_search(residues.begin(), residues.end(), index)) {
kept.push_back(std::move(r));
out.residue_index.push_back(index);
}
}
ch.residues = std::move(kept);
out.residues_kept += static_cast<int>(ch.residues.size());
}
std::erase_if(out.st.models[0].chains, [](const gemmi::Chain &ch) { return ch.residues.empty(); });
finish_search_model(out);
// The domain's residues as ranges of the prediction's numbering, from 1.
std::string ranges;
for (size_t k = 0; k < residues.size(); k++) {
if (k > 0 && residues[k] == residues[k - 1] + 1 && k + 1 < residues.size() && residues[k + 1] == residues[k] + 1)
continue;
if (k > 0 && residues[k] == residues[k - 1] + 1)
ranges += fmt::format("-{}", residues[k] + 1);
else
ranges += fmt::format("{}{}", ranges.empty() ? "" : ",", residues[k] + 1);
}
out.description = fmt::format("PAE domain of {} residues ({})", out.residues_kept, ranges);
return out;
}
// --- rotations ---
std::vector<gemmi::Mat33> UniformRotations(double step_deg) {
const double step = step_deg * PI / 180.0;
const double golden = PI * (3.0 - std::sqrt(5.0));
const int shells = std::max(1, static_cast<int>(std::lround(PI / step)));
std::vector<gemmi::Mat33> out{gemmi::Mat33()};
for (int i = 1; i <= shells; i++) {
const double kappa = PI * i / shells;
const double reach = step / (2.0 * std::sin(kappa / 2.0)); // the axis spacing this shell needs
int n = std::max(1, static_cast<int>(std::lround(4.0 * PI / (reach * reach))));
// At 180 deg an axis and its opposite are the same rotation: one hemisphere.
const bool half = i == shells;
if (half)
n = (n + 1) / 2;
for (int k = 0; k < n; k++) {
const double z = half ? 1.0 - (k + 0.5) / n : 1.0 - (2.0 * k + 1.0) / n;
const double r = std::sqrt(std::max(0.0, 1.0 - z * z));
const double x = r * std::cos(golden * k), y = r * std::sin(golden * k);
const double s = std::sin(kappa), c1 = 1.0 - std::cos(kappa);
out.push_back({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)});
}
}
return out;
}
std::vector<gemmi::Mat33> CartesianPointGroup(const gemmi::SpaceGroup &sg, const gemmi::UnitCell &cell) {
std::vector<gemmi::Mat33> out;
for (const gemmi::Op &op : sg.operations().sym_ops) {
gemmi::Mat33 f;
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
f.a[i][j] = static_cast<double>(op.rot[i][j]) / gemmi::Op::DEN;
out.push_back(cell.orth.mat.multiply(f).multiply(cell.frac.mat));
}
return out;
}
std::vector<gemmi::Mat33> RotationAsu(const std::vector<gemmi::Mat33> &all,
const std::vector<gemmi::Mat33> &point_group) {
std::vector<gemmi::Mat33> out;
for (const gemmi::Mat33 &r : all) {
const double t = trace(r);
if (std::all_of(point_group.begin(), point_group.end(),
[&](const gemmi::Mat33 &s) { return t >= trace(s.multiply(r)) - 1e-9; }))
out.push_back(r);
}
return out;
}
double RotationDistanceDeg(const gemmi::Mat33 &a, const gemmi::Mat33 &b,
const std::vector<gemmi::Mat33> &point_group) {
double best = -1.0;
for (const gemmi::Mat33 &s : point_group)
best = std::max(best, trace(a.multiply(s.multiply(b).transpose())));
return std::acos(std::clamp((best - 1.0) / 2.0, -1.0, 1.0)) * 180.0 / PI;
}
// --- the molecular transform ---
MolecularTransform::MolecularTransform(const SearchModel &model, double d_min, size_t nthreads) {
// On a P1 box six times the model's radius: the transform of a body of radius r turns its phase
// by pi r / box per grid step, and trilinear interpolation reads it to a few percent at pi/6.
// Interpolating the molecular transform of one model in place of a Fourier summation per
// orientation, following Navaza (1994) Acta Cryst. A50, 157-163
box_ = std::max(6.0 * model.radius_A, 120.0);
gemmi::Structure st = model.st;
st.cell = gemmi::UnitCell(box_, box_, box_, 90, 90, 90);
st.spacegroup_hm = "P 1";
st.setup_cell_images();
gemmi::DensityCalculator<Table, float> dc;
dc.d_min = d_min;
dc.rate = 1.5;
dc.set_grid_cell_and_spacegroup(st);
dc.set_refmac_compatible_blur(st.models[0]);
PutModelDensityOnGrid(dc, st.models[0], {}, nthreads);
#ifdef JFJOCH_USE_CUDA
if (get_gpu_count() > 0 && MapToFPhiGPUFits(dc.grid))
grid_ = MapToFPhiGPU(dc.grid);
else
#endif
grid_ = MapToFPhi(dc.grid);
d_min_ = d_min;
blur_ = dc.blur;
}
std::complex<double> MolecularTransform::At(const gemmi::Vec3 &s) const {
const double u[3] = {s.x * box_, s.y * box_, s.z * box_};
const int i0[3] = {static_cast<int>(std::floor(u[0])), static_cast<int>(std::floor(u[1])),
static_cast<int>(std::floor(u[2]))};
const double f[3] = {u[0] - i0[0], u[1] - i0[1], u[2] - i0[2]};
const auto get = [&](int i, int j, int l) {
bool conj = false;
if (l < 0) { // the half-l grid holds l >= 0; F(-h) = conj F(h)
i = -i; j = -j; l = -l;
conj = true;
}
i %= grid_.nu; if (i < 0) i += grid_.nu;
j %= grid_.nv; if (j < 0) j += grid_.nv;
const std::complex<float> v = grid_.data[i + static_cast<size_t>(grid_.nu) * (j + static_cast<size_t>(grid_.nv) * l)];
return conj ? std::conj(std::complex<double>(v)) : std::complex<double>(v);
};
std::complex<double> out = 0;
for (int dx = 0; dx < 2; dx++)
for (int dy = 0; dy < 2; dy++)
for (int dz = 0; dz < 2; dz++) {
const double w = (dx ? f[0] : 1 - f[0]) * (dy ? f[1] : 1 - f[1]) * (dz ? f[2] : 1 - f[2]);
out += w * get(i0[0] + dx, i0[1] + dy, i0[2] + dz);
}
return out * std::exp(blur_ * s.length_sq() / 4.0);
}
// --- the observations ---
SearchReflections MakeSearchReflections(const std::vector<MergedReflection> &merged, const gemmi::UnitCell &cell,
const gemmi::SpaceGroup &sg, double d_min, double d_max) {
const gemmi::GroupOps gops = sg.operations();
const gemmi::ReciprocalAsu asu(&sg);
struct Row { gemmi::Miller h; double i_eps; double s2; double mult; };
std::vector<Row> rows;
for (const MergedReflection &r : merged) {
if (!std::isfinite(r.I) || r.d < d_min || r.d > d_max)
continue;
gemmi::Miller h{{r.h, r.k, r.l}};
if (!asu.is_in(h))
h = asu.to_asu(h, gops).first;
if (gops.is_systematically_absent(h))
continue;
// The reflections the whole sphere holds for this one: its images under the point group and
// their Friedel mates, each counted once.
std::vector<gemmi::Miller> images;
for (const gemmi::Op &op : gops.sym_ops)
for (int sign : {1, -1}) {
gemmi::Miller k = rotated_index(op, h);
for (int &x : k)
x *= sign;
if (std::find(images.begin(), images.end(), k) == images.end())
images.push_back(k);
}
rows.push_back({h, r.I / gops.epsilon_factor_without_centering(h), cell.calculate_1_d2(h),
static_cast<double>(images.size())});
}
std::stable_sort(rows.begin(), rows.end(), [](const Row &a, const Row &b) { return a.s2 < b.s2; });
SearchReflections out;
const size_t n = rows.size();
out.bins = static_cast<int>(std::min<size_t>(SEARCH_BINS, std::max<size_t>(n / 100, 1)));
std::vector<double> sum(out.bins, 0.0);
std::vector<int> count(out.bins, 0);
for (size_t i = 0; i < n; i++) {
const int b = static_cast<int>(i * out.bins / n);
out.hkl.push_back(rows[i].h);
out.multiplicity.push_back(rows[i].mult);
out.bin.push_back(b);
sum[b] += rows[i].i_eps;
count[b]++;
}
for (size_t i = 0; i < n; i++) {
const double mean = sum[out.bin[i]] / count[out.bin[i]];
out.e2.push_back(mean > 0 ? rows[i].i_eps / mean : 0.0);
}
return out;
}
// --- the rotation function ---
std::vector<gemmi::Vec3> CopyVectors(const SearchReflections &refl, const gemmi::UnitCell &cell,
const gemmi::SpaceGroup &sg) {
const gemmi::GroupOps gops = sg.operations();
const size_t nops = gops.sym_ops.size();
const gemmi::Mat33 frac_t = cell.frac.mat.transpose();
std::vector<gemmi::Vec3> s(refl.hkl.size() * nops);
for (size_t i = 0; i < refl.hkl.size(); i++)
for (size_t j = 0; j < nops; j++)
s[i * nops + j] = frac_t.multiply(miller_vec(rotated_index(gops.sym_ops[j], refl.hkl[i])));
return s;
}
std::vector<double> RotationFunction(const MolecularTransform &ft, const SearchReflections &refl,
const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
const std::vector<gemmi::Mat33> &rotations, size_t nthreads) {
const size_t n = refl.hkl.size();
const std::vector<gemmi::Vec3> s = CopyVectors(refl, cell, sg);
const size_t nops = s.size() / n;
double w_mean = 0;
for (double e2 : refl.e2)
w_mean += e2;
w_mean /= static_cast<double>(n);
double w_ss = 0;
for (double e2 : refl.e2)
w_ss += (e2 - w_mean) * (e2 - w_mean);
std::vector<double> out(rotations.size(), 0.0);
ParallelFor(static_cast<int>(rotations.size()), nthreads, [&](int r) {
// The model turned by R has F_R(s) = F_model(R^T s).
const gemmi::Mat33 rt = rotations[r].transpose();
std::vector<double> ic(n, 0.0), bin_sum(refl.bins, 0.0);
std::vector<int> bin_n(refl.bins, 0);
for (size_t i = 0; i < n; i++) {
for (size_t j = 0; j < nops; j++)
ic[i] += std::norm(ft.At(rt.multiply(s[i * nops + j])));
bin_sum[refl.bin[i]] += ic[i];
bin_n[refl.bin[i]]++;
}
double x_sum = 0, x_ss = 0, xy = 0;
for (size_t i = 0; i < n; i++) {
const double x = ic[i] / (bin_sum[refl.bin[i]] / bin_n[refl.bin[i]]);
x_sum += x;
x_ss += x * x;
xy += x * (refl.e2[i] - w_mean);
}
const double sxx = x_ss - x_sum * x_sum / static_cast<double>(n);
out[r] = sxx > 0 && w_ss > 0 ? xy / std::sqrt(sxx * w_ss) : 0.0;
});
return out;
}
// --- the translation function ---
//
// The full-symmetry translation function as one Fourier synthesis, following Crowther & Blow (1967)
// Acta Cryst. 23, 544-548: with the model at t, copy j adds G_j(h) exp(2 pi i (h R_j).t) exp(2 pi i h.t_j)
// to F(h), so |F(h)|^2 summed against E^2 - 1 over the sphere is, for each pair of copies j != k, a
// Fourier series in t with coefficient (E^2 - 1) g_j conj(g_k), g_j = G_j(h) exp(2 pi i h.t_j), at the
// index h R_j - h R_k. Every pair's terms go on one grid, and one FFT gives the function at every t.
// The pairs j = k are the rotation function's, constant in t.
TranslationSetup MakeTranslationSetup(const SearchReflections &refl, const gemmi::UnitCell &cell,
const gemmi::SpaceGroup &sg, bool known_part) {
const gemmi::GroupOps gops = sg.operations();
const int nops = static_cast<int>(gops.sym_ops.size());
const gemmi::Mat33 frac_t = cell.frac.mat.transpose();
TranslationSetup out;
out.ops = nops;
out.bins = refl.bins;
for (int j = 0; j < nops; j++)
for (int k = 0; k < nops; k++)
if (j != k) {
out.pair_j.push_back(j);
out.pair_k.push_back(k);
}
std::vector<gemmi::Miller> rotated; // terms x ops
for (size_t i = 0; i < refl.hkl.size(); i++) {
const gemmi::Miller &h = refl.hkl[i];
if (gops.is_systematically_absent(h))
continue;
for (const gemmi::Op &op : gops.sym_ops) {
rotated.push_back(rotated_index(op, h));
out.s.push_back(frac_t.multiply(miller_vec(rotated.back())));
const double ht = (h[0] * op.tran[0] + h[1] * op.tran[1] + h[2] * op.tran[2])
/ static_cast<double>(gemmi::Op::DEN);
out.phase.push_back(std::polar(1.0, 2.0 * PI * ht));
}
out.weight.push_back(refl.multiplicity[i] * (refl.e2[i] - 1.0));
out.e2.push_back(refl.e2[i]);
out.multiplicity.push_back(refl.multiplicity[i]);
out.bin.push_back(refl.bin[i]);
}
int max_m[3] = {0, 0, 0};
for (size_t t = 0; t < out.Terms(); t++)
for (size_t p = 0; p < out.pair_j.size(); p++)
for (int a = 0; a < 3; a++)
max_m[a] = std::max(max_m[a], std::abs(rotated[t * nops + out.pair_j[p]][a]
- rotated[t * nops + out.pair_k[p]][a]));
if (known_part)
for (const gemmi::Miller &k : rotated)
for (int a = 0; a < 3; a++)
max_m[a] = std::max(max_m[a], std::abs(k[a]));
out.nu = good_even_size(2 * max_m[0] + 2);
out.nv = good_even_size(2 * max_m[1] + 2);
out.nw = good_even_size(2 * max_m[2] + 2);
for (size_t t = 0; t < out.Terms(); t++)
for (size_t p = 0; p < out.pair_j.size(); p++) {
const gemmi::Miller &kj = rotated[t * nops + out.pair_j[p]], &kk = rotated[t * nops + out.pair_k[p]];
const int l = kj[2] - kk[2];
out.target_u.push_back(l < 0 ? -1 : ((kj[0] - kk[0]) % out.nu + out.nu) % out.nu);
out.target_v.push_back(((kj[1] - kk[1]) % out.nv + out.nv) % out.nv);
out.target_l.push_back(l);
}
if (known_part)
for (size_t t = 0; t < out.Terms(); t++)
for (int j = 0; j < nops; j++) {
const gemmi::Miller &k = rotated[t * nops + j];
const auto add = [&](int sign) {
out.cross_term.push_back(static_cast<int>(t));
out.cross_op.push_back(j);
out.cross_u.push_back(((sign * k[0]) % out.nu + out.nu) % out.nu);
out.cross_v.push_back(((sign * k[1]) % out.nv + out.nv) % out.nv);
out.cross_l.push_back(sign * k[2]);
out.cross_conj.push_back(sign < 0);
};
if (k[2] >= 0)
add(1);
if (k[2] <= 0)
add(-1);
}
return out;
}
std::vector<std::complex<double>> KnownStructureFactors(const MolecularTransform &ft, const TranslationSetup &setup,
const gemmi::UnitCell &cell,
const std::vector<gemmi::Mat33> &rotations,
const std::vector<gemmi::Fractional> &positions) {
const size_t terms = setup.Terms(), nops = setup.ops;
std::vector<std::complex<double>> out(terms, 0.0);
for (size_t c = 0; c < rotations.size(); c++) {
const gemmi::Mat33 rt = rotations[c].transpose();
const gemmi::Position x = cell.orthogonalize(positions[c]);
for (size_t t = 0; t < terms; t++)
for (size_t j = 0; j < nops; j++) {
const gemmi::Vec3 &sj = setup.s[t * nops + j];
// (h R_j).t = s_j.x for the Cartesian s_j and x
out[t] += ft.At(rt.multiply(sj)) * setup.phase[t * nops + j] * std::polar(1.0, 2.0 * PI * sj.dot(x));
}
}
return out;
}
TranslationPeak TranslationFunction(const MolecularTransform &ft, const TranslationSetup &setup,
const gemmi::Mat33 &rotation, const std::vector<std::complex<double>> *known) {
TranslationPeak out;
const size_t terms = setup.Terms(), nops = setup.ops, npairs = setup.pair_j.size();
if (terms == 0 || (nops < 2 && known == nullptr))
return out; // one copy in P1: every position is the same
const gemmi::Mat33 rt = rotation.transpose();
std::vector<std::complex<double>> g(terms * nops);
std::vector<double> bin_sum(setup.bins, 0.0);
std::vector<int> bin_n(setup.bins, 0);
for (size_t t = 0; t < terms; t++) {
double sum = 0;
for (size_t j = 0; j < nops; j++) {
g[t * nops + j] = ft.At(rt.multiply(setup.s[t * nops + j])) * setup.phase[t * nops + j];
sum += std::norm(g[t * nops + j]);
}
bin_sum[setup.bin[t]] += sum;
bin_n[setup.bin[t]]++;
}
const int nu = setup.nu, nv = setup.nv, nwh = setup.nw / 2 + 1;
std::vector<std::complex<double>> c(static_cast<size_t>(nu) * nv * nwh, 0.0);
for (size_t t = 0; t < terms; t++) {
const double w = setup.weight[t] / (bin_sum[setup.bin[t]] / bin_n[setup.bin[t]]);
for (size_t p = 0; p < npairs; p++) {
const size_t e = t * npairs + p;
if (setup.target_u[e] < 0)
continue;
c[setup.target_u[e] + static_cast<size_t>(nu) * (setup.target_v[e] + static_cast<size_t>(nv) * setup.target_l[e])]
+= w * g[t * nops + setup.pair_j[p]] * std::conj(g[t * nops + setup.pair_k[p]]);
}
}
if (known != nullptr)
for (size_t e = 0; e < setup.cross_term.size(); e++) {
const size_t t = setup.cross_term[e];
const double w = setup.weight[t] / (bin_sum[setup.bin[t]] / bin_n[setup.bin[t]]);
const std::complex<double> v = std::conj((*known)[t]) * g[t * nops + setup.cross_op[e]];
c[setup.cross_u[e] + static_cast<size_t>(nu) * (setup.cross_v[e] + static_cast<size_t>(nv) * setup.cross_l[e])]
+= w * (setup.cross_conj[e] ? std::conj(v) : v);
}
// MapFromFPhi sums F exp(-2 pi i m.x): with F = conj(C) that is the conjugate of the series, which
// is real.
gemmi::FPhiGrid<float> hkl;
hkl.unit_cell = ft.Grid().unit_cell;
hkl.spacegroup = &gemmi::get_spacegroup_p1();
hkl.axis_order = gemmi::AxisOrder::XYZ;
hkl.half_l = true;
hkl.set_size_without_checking(nu, nv, nwh);
for (size_t i = 0; i < c.size(); i++)
hkl.data[i] = std::complex<float>(std::conj(c[i]));
const gemmi::Grid<float> map = MapFromFPhi(hkl);
double sum = 0, sum2 = 0;
size_t best = 0;
for (size_t i = 0; i < map.data.size(); i++) {
sum += map.data[i];
sum2 += static_cast<double>(map.data[i]) * map.data[i];
if (map.data[i] > map.data[best])
best = i;
}
const double count = static_cast<double>(map.data.size());
const double mean = sum / count, sd = std::sqrt(std::max(0.0, sum2 / count - mean * mean));
out.z = sd > 0 ? (map.data[best] - mean) / sd : 0.0;
const size_t u = best % map.nu, v = (best / map.nu) % map.nv, w = best / (static_cast<size_t>(map.nu) * map.nv);
out.t = gemmi::Fractional(static_cast<double>(u) / map.nu, static_cast<double>(v) / map.nv,
static_cast<double>(w) / map.nw);
return out;
}
// --- the likelihood ---
//
// The log-likelihood gain of a placement over a random structure, from the Rice (acentric) and Woolfson
// (centric) distributions of E_obs given E_calc with sigma_A, and sigma_A from the model's expected
// coordinate error, following Read (2001) Acta Cryst. D57, 1373-1382; the gain judged as a Z-score over
// the same orientation at other positions, as the likelihood-enhanced translation function's Z-score of
// McCoy et al. (2005) Acta Cryst. D61, 458-464.
namespace {
// ln I0(x), x >= 0. Following Abramowitz & Stegun (1964) 9.8.1-9.8.2
double log_bessel_i0(double x) {
if (x < 3.75) {
const double t = (x / 3.75) * (x / 3.75);
return std::log(1.0 + t * (3.5156229 + t * (3.0899424 + t * (1.2067492 + t * (0.2659732
+ t * (0.0360768 + t * 0.0045813))))));
}
const double u = 3.75 / x;
return x - 0.5 * std::log(x) + std::log(0.39894228 + u * (0.01328592 + u * (0.00225319 + u * (-0.00157565
+ u * (0.00916281 + u * (-0.02057706 + u * (0.02635537 + u * (-0.01647633 + u * 0.00392377))))))));
}
double log_cosh(double x) {
x = std::fabs(x);
return x + std::log1p(std::exp(-2.0 * x)) - std::log(2.0);
}
constexpr int LIKELIHOOD_NULL_POSITIONS = 100;
constexpr unsigned LIKELIHOOD_SEED = 20261012;
} // namespace
LikelihoodSetup MakeLikelihoodSetup(const SearchReflections &refl, const gemmi::UnitCell &cell,
const gemmi::SpaceGroup &sg) {
const gemmi::GroupOps gops = sg.operations();
const gemmi::Mat33 frac_t = cell.frac.mat.transpose();
LikelihoodSetup out;
out.ops = static_cast<int>(gops.sym_ops.size());
out.bins = refl.bins;
for (size_t i = 0; i < refl.hkl.size(); i++) {
const gemmi::Miller &h = refl.hkl[i];
if (gops.is_systematically_absent(h))
continue;
for (const gemmi::Op &op : gops.sym_ops) {
const gemmi::Miller k = rotated_index(op, h);
out.k.push_back(k);
out.s.push_back(frac_t.multiply(miller_vec(k)));
const double ht = (h[0] * op.tran[0] + h[1] * op.tran[1] + h[2] * op.tran[2])
/ static_cast<double>(gemmi::Op::DEN);
out.phase.push_back(std::polar(1.0, 2.0 * PI * ht));
for (int a = 0; a < 3; a++)
out.max_k[a] = std::max(out.max_k[a], std::abs(k[a]));
}
out.eo.push_back(std::sqrt(std::max(refl.e2[i], 0.0)));
out.epsilon.push_back(gops.epsilon_factor_without_centering(h));
out.s2.push_back(cell.calculate_1_d2(h));
out.centric.push_back(gops.is_reflection_centric(h));
out.bin.push_back(refl.bin[i]);
}
return out;
}
std::vector<std::complex<double>> CopyFactors(const SearchModel &model, const LikelihoodSetup &setup,
const gemmi::UnitCell &cell, const gemmi::Mat33 &rotation,
double d_min, size_t nthreads, const MolecularTransform *ft) {
std::vector<std::complex<double>> g(setup.s.size());
if (ft != nullptr && d_min >= ft->DMin()) {
const gemmi::Mat33 rt = rotation.transpose();
for (size_t e = 0; e < g.size(); e++)
g[e] = ft->At(rt.multiply(setup.s[e])) * setup.phase[e];
return g;
}
gemmi::Structure st = model.st;
std::vector<gemmi::Position> pos = ModelPositions(st.models[0]);
for (gemmi::Position &p : pos)
p = gemmi::Position(rotation.multiply(p));
SetModelPositions(st.models[0], pos);
st.cell = cell;
st.spacegroup_hm = "P 1";
st.setup_cell_images();
gemmi::DensityCalculator<Table, float> dc;
dc.d_min = d_min;
dc.rate = 1.5;
dc.set_grid_cell_and_spacegroup(st);
dc.set_refmac_compatible_blur(st.models[0]);
PutModelDensityOnGrid(dc, st.models[0], {}, nthreads);
const gemmi::FPhiGrid<float> f1 = MapToFPhi(dc.grid);
for (size_t e = 0; e < g.size(); e++) {
int i = setup.k[e][0], j = setup.k[e][1], l = setup.k[e][2];
bool conj = false;
if (l < 0) { // the half-l grid holds l >= 0; F(-h) = conj F(h)
i = -i; j = -j; l = -l;
conj = true;
}
i %= f1.nu; if (i < 0) i += f1.nu;
j %= f1.nv; if (j < 0) j += f1.nv;
const std::complex<double> v = f1.data[i + static_cast<size_t>(f1.nu) * (j + static_cast<size_t>(f1.nv) * l)];
g[e] = (conj ? std::conj(v) : v) * std::exp(dc.blur * setup.s[e].length_sq() / 4.0) * setup.phase[e];
}
return g;
}
void AddCopy(const LikelihoodSetup &setup, const std::vector<std::complex<double>> &g, const gemmi::Fractional &t,
std::vector<std::complex<double>> &f) {
// exp(2 pi i m t_a) for every index m an axis needs, so that each term is three products.
std::array<std::vector<std::complex<double>>, 3> table;
const double ta[3] = {t.x, t.y, t.z};
for (int a = 0; a < 3; a++) {
const int m = setup.max_k[a];
table[a].resize(2 * m + 1);
for (int i = -m; i <= m; i++)
table[a][i + m] = std::polar(1.0, 2.0 * PI * i * ta[a]);
}
const size_t nops = setup.ops;
for (size_t i = 0; i < setup.Terms(); i++)
for (size_t j = 0; j < nops; j++) {
const size_t e = i * nops + j;
const gemmi::Miller &k = setup.k[e];
f[i] += g[e] * table[0][k[0] + setup.max_k[0]] * table[1][k[1] + setup.max_k[1]]
* table[2][k[2] + setup.max_k[2]];
}
}
namespace {
// E_calc^2 of each term: |f|^2 / epsilon over its bin's mean (0 where the bin has none).
std::vector<double> calc_e2(const LikelihoodSetup &setup, const std::vector<std::complex<double>> &f) {
const size_t n = setup.Terms();
std::vector<double> bin_sum(setup.bins, 0.0);
std::vector<int> bin_n(setup.bins, 0);
std::vector<double> e2(n);
for (size_t i = 0; i < n; i++) {
e2[i] = std::norm(f[i]) / setup.epsilon[i];
bin_sum[setup.bin[i]] += e2[i];
bin_n[setup.bin[i]]++;
}
for (size_t i = 0; i < n; i++) {
const double mean = bin_sum[setup.bin[i]] / bin_n[setup.bin[i]];
e2[i] = mean > 0 ? e2[i] / mean : 0.0;
}
return e2;
}
// sigma_A's fall-off with resolution for the model's coordinate error, Luzzati's D, at each term.
std::vector<double> error_factor(const LikelihoodSetup &setup, double error_A) {
std::vector<double> out(setup.Terms());
for (size_t i = 0; i < out.size(); i++)
out[i] = std::exp(-2.0 * PI * PI * error_A * error_A * setup.s2[i] / 3.0);
return out;
}
double llg_sum(const LikelihoodSetup &setup, const std::vector<double> &ec2, const std::vector<double> &d,
double scale) {
double llg = 0;
for (size_t i = 0; i < setup.Terms(); i++) {
const double eo = setup.eo[i];
const double sa = scale * d[i];
const double v = 1.0 - sa * sa;
const double x = 2.0 * sa * eo * std::sqrt(ec2[i]) / v;
if (setup.centric[i])
llg += -0.5 * std::log(v) - sa * sa * (eo * eo + ec2[i]) / (2.0 * v) + log_cosh(x / 2.0);
else
llg += -std::log(v) - sa * sa * (eo * eo + ec2[i]) / v + log_bessel_i0(x);
}
return llg;
}
} // namespace
double LogLikelihoodGain(const LikelihoodSetup &setup, const std::vector<std::complex<double>> &f, double error_A,
double scale) {
return llg_sum(setup, calc_e2(setup, f), error_factor(setup, error_A), scale);
}
double RefinedLogLikelihoodGain(const LikelihoodSetup &setup, const std::vector<std::complex<double>> &f,
double error_A, double *scale) {
// Golden-section search for the scale in (0, 0.99), to a thousandth: the gain is 0 at 0 and has one
// maximum.
const std::vector<double> ec2 = calc_e2(setup, f), d_error = error_factor(setup, error_A);
const double golden = (std::sqrt(5.0) - 1.0) / 2.0;
double a = 0.0, b = 0.99;
double c = b - golden * (b - a), d = a + golden * (b - a);
double fc = llg_sum(setup, ec2, d_error, c), fd = llg_sum(setup, ec2, d_error, d);
for (int iter = 0; iter < 14; iter++) {
if (fc > fd) {
b = d; d = c; fd = fc;
c = b - golden * (b - a);
fc = llg_sum(setup, ec2, d_error, c);
} else {
a = c; c = d; fc = fd;
d = a + golden * (b - a);
fd = llg_sum(setup, ec2, d_error, d);
}
}
if (scale != nullptr)
*scale = fc > fd ? c : d;
return std::max(fc, fd);
}
double LikelihoodGain(const LikelihoodSetup &setup, const std::vector<std::complex<double>> &known, double known_llg,
const std::vector<std::complex<double>> &g, const gemmi::Fractional &t, double error_A,
double *llg, double *scale) {
std::vector<std::complex<double>> f = known;
AddCopy(setup, g, t, f);
const double all = RefinedLogLikelihoodGain(setup, f, error_A, scale);
if (llg != nullptr)
*llg = all;
return all - known_llg;
}
LikelihoodPeak LikelihoodTranslationZ(const LikelihoodSetup &setup, const std::vector<std::complex<double>> &known,
double known_llg, const std::vector<std::complex<double>> &g,
const gemmi::Fractional &t, double error_A, size_t nthreads) {
LikelihoodPeak out;
out.gain = LikelihoodGain(setup, known, known_llg, g, t, error_A, &out.llg, &out.scale);
// The positions drawn from the raw 32-bit stream, which every standard library makes the same.
std::mt19937 rng(LIKELIHOOD_SEED);
std::vector<gemmi::Fractional> where(LIKELIHOOD_NULL_POSITIONS);
for (gemmi::Fractional &w : where) {
const double x = rng() / 4294967296.0, y = rng() / 4294967296.0, z = rng() / 4294967296.0;
w = gemmi::Fractional(x, y, z);
}
std::vector<double> null(where.size());
const std::vector<double> d_error = error_factor(setup, error_A);
ParallelFor(static_cast<int>(where.size()), nthreads, [&](int i) {
std::vector<std::complex<double>> fr = known;
AddCopy(setup, g, where[i], fr);
null[i] = llg_sum(setup, calc_e2(setup, fr), d_error, out.scale);
});
double mean = 0, s2 = 0;
for (double x : null)
mean += x;
mean /= static_cast<double>(null.size());
for (double x : null)
s2 += (x - mean) * (x - mean);
const double sd = std::sqrt(s2 / static_cast<double>(null.size() - 1));
out.z = sd > 0 ? (out.llg - mean) / sd : 0.0;
return out;
}
// --- translational NCS ---
TranslationalNcs DetectTranslationalNcs(const std::vector<MergedReflection> &merged, const gemmi::UnitCell &cell,
const gemmi::SpaceGroup &sg) {
constexpr double D_MIN = 5.0, D_MAX = 10.0, MIN_DISTANCE_A = 15.0, MIN_HEIGHT = 0.2;
TranslationalNcs out;
const SearchReflections refl = MakeSearchReflections(merged, cell, sg, D_MIN, D_MAX);
const gemmi::GroupOps gops = sg.operations();
int max_m[3] = {0, 0, 0};
std::vector<std::pair<gemmi::Miller, double>> terms; // the whole sphere, l >= 0
for (size_t i = 0; i < refl.hkl.size(); i++)
for (const gemmi::Op &op : gops.sym_ops)
for (int sign : {1, -1}) {
gemmi::Miller k = rotated_index(op, refl.hkl[i]);
for (int &x : k)
x *= sign;
if (k[2] < 0)
continue;
terms.push_back({k, refl.e2[i]});
for (int a = 0; a < 3; a++)
max_m[a] = std::max(max_m[a], std::abs(k[a]));
}
if (terms.empty())
return out;
// Sampled at a third of D_MIN at least.
gemmi::FPhiGrid<float> hkl;
hkl.unit_cell = cell;
hkl.spacegroup = &gemmi::get_spacegroup_p1();
hkl.axis_order = gemmi::AxisOrder::XYZ;
hkl.half_l = true;
const int nu = good_even_size(std::max(2 * max_m[0] + 2, static_cast<int>(std::ceil(3 * cell.a / D_MIN))));
const int nv = good_even_size(std::max(2 * max_m[1] + 2, static_cast<int>(std::ceil(3 * cell.b / D_MIN))));
const int nw = good_even_size(std::max(2 * max_m[2] + 2, static_cast<int>(std::ceil(3 * cell.c / D_MIN))));
hkl.set_size_without_checking(nu, nv, nw / 2 + 1);
for (const auto &[k, e2] : terms) {
const size_t idx = ((k[0] % nu + nu) % nu) + static_cast<size_t>(nu) * (((k[1] % nv + nv) % nv)
+ static_cast<size_t>(nv) * k[2]);
hkl.data[idx] = std::complex<float>(static_cast<float>(e2), 0.0f);
}
const gemmi::Grid<float> map = MapFromFPhi(hkl);
const double origin = map.data[0];
for (int w = 0; w < map.nw; w++)
for (int v = 0; v < map.nv; v++)
for (int u = 0; u < map.nu; u++) {
const double value = map.data[u + static_cast<size_t>(map.nu) * (v + static_cast<size_t>(map.nv) * w)];
if (!(value > out.height * origin))
continue;
const gemmi::Fractional f(static_cast<double>(u) / map.nu, static_cast<double>(v) / map.nv,
static_cast<double>(w) / map.nw);
bool near_lattice = false;
for (const gemmi::Op::Tran &c : gops.cen_ops) {
gemmi::Fractional d(f.x - static_cast<double>(c[0]) / gemmi::Op::DEN,
f.y - static_cast<double>(c[1]) / gemmi::Op::DEN,
f.z - static_cast<double>(c[2]) / gemmi::Op::DEN);
d = gemmi::Fractional(d.x - std::round(d.x), d.y - std::round(d.y), d.z - std::round(d.z));
if (cell.orthogonalize_difference(d).length() < MIN_DISTANCE_A)
near_lattice = true;
}
if (near_lattice)
continue;
out.height = value / origin;
out.vector = f;
}
out.present = out.height >= MIN_HEIGHT;
return out;
}
bool TncsRelatesGroups(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b, const gemmi::Fractional &v) {
constexpr double TOLERANCE = 0.05;
const auto mod1 = [](double x) { return x - std::round(x); };
const gemmi::GroupOps bops = b.operations();
for (const gemmi::Op &op : a.operations().sym_ops) {
const gemmi::Op *other = bops.find_by_rotation(op.rot);
if (other == nullptr)
return false;
// The screw parts, (1/n) sum over k of R^k t, and the same of v.
double wa[3] = {0, 0, 0}, wb[3] = {0, 0, 0}, wv[3] = {0, 0, 0};
gemmi::Op power = gemmi::Op::identity();
int n = 0;
do {
for (int i = 0; i < 3; i++) {
double ta = 0, tb = 0, tv = 0;
for (int j = 0; j < 3; j++) {
const double r = static_cast<double>(power.rot[i][j]) / gemmi::Op::DEN;
ta += r * op.tran[j] / gemmi::Op::DEN;
tb += r * other->tran[j] / gemmi::Op::DEN;
tv += r * (j == 0 ? v.x : j == 1 ? v.y : v.z);
}
wa[i] += ta;
wb[i] += tb;
wv[i] += tv;
}
power = power.combine(gemmi::Op{op.rot, {0, 0, 0}, 'x'});
n++;
} while (!(power.rot == gemmi::Op::identity().rot) && n < 6);
bool bridged = false;
for (int m = 0; m < n && !bridged; m++) {
bridged = true;
for (int i = 0; i < 3; i++)
if (std::fabs(mod1((wa[i] - wb[i] - m * wv[i]) / n)) > TOLERANCE)
bridged = false;
}
if (!bridged)
return false;
}
return true;
}
// --- placements ---
gemmi::Structure PlaceSearchModel(const SearchModel &model, const gemmi::Mat33 &rotation, const gemmi::Fractional &t,
const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg) {
gemmi::Structure st = model.st;
const gemmi::Position shift = cell.orthogonalize(t);
std::vector<gemmi::Position> pos = ModelPositions(st.models[0]);
for (gemmi::Position &p : pos)
p = gemmi::Position(rotation.multiply(p) + shift);
SetModelPositions(st.models[0], pos);
st.cell = cell;
st.spacegroup_hm = sg.xhm();
st.setup_cell_images();
return st;
}
double ClashFraction(const gemmi::Structure &placed, size_t copies) {
const gemmi::UnitCell &cell = placed.cell;
const gemmi::SpaceGroup *sg = placed.find_spacegroup();
const size_t chains = placed.models[0].chains.size();
copies = std::max<size_t>(1, std::min(copies, chains));
// Each copy's C-alpha atoms, fractional, its centroid and its radius.
struct Copy {
std::vector<gemmi::Fractional> ca;
gemmi::Fractional c0;
double radius = 0;
};
std::vector<Copy> body(copies);
for (size_t i = 0; i < chains; i++)
for (const gemmi::Residue &r : placed.models[0].chains[i].residues)
for (const gemmi::Atom &a : r.atoms)
if (a.name == "CA")
body[i * copies / chains].ca.push_back(cell.fractionalize(a.pos));
size_t total = 0;
for (Copy &b : body) {
total += b.ca.size();
for (const gemmi::Fractional &f : b.ca)
b.c0 = b.c0 + f;
if (!b.ca.empty())
b.c0 = gemmi::Fractional(b.c0 * (1.0 / static_cast<double>(b.ca.size())));
for (const gemmi::Fractional &f : b.ca)
b.radius = std::max(b.radius, cell.orthogonalize_difference(f - b.c0).length());
}
if (total == 0 || sg == nullptr)
return 0.0;
constexpr double CLASH_A = 3.0;
// Every C-alpha of every image that can reach a copy - each copy under each operator and the lattice
// shifts about it - binned in CLASH_A cubes, sorted by cube; a C-alpha clashes where one of another
// image (not its own copy unmoved) sits within CLASH_A of it in its own cube or a neighbouring one.
struct Image {
long key;
gemmi::Position x;
int copy;
bool unmoved;
};
const auto cube = [](const gemmi::Position &x, int dx, int dy, int dz) {
return (static_cast<long>(std::floor(x.x / CLASH_A)) + dx + 4096) * 67108864L
+ (static_cast<long>(std::floor(x.y / CLASH_A)) + dy + 4096) * 8192L
+ (static_cast<long>(std::floor(x.z / CLASH_A)) + dz + 4096);
};
std::vector<Image> images;
for (size_t b = 0; b < copies; b++)
for (const gemmi::Op &op : sg->operations().all_ops_sorted()) {
std::vector<gemmi::Fractional> mate;
for (const gemmi::Fractional &f : body[b].ca) {
const std::array<double, 3> x = op.apply_to_xyz({{f.x, f.y, f.z}});
mate.push_back(gemmi::Fractional(x[0], x[1], x[2]));
}
gemmi::Fractional cm;
for (const gemmi::Fractional &f : mate)
cm = cm + f;
cm = gemmi::Fractional(cm * (1.0 / static_cast<double>(mate.size())));
// The lattice shifts that bring this image within reach of some copy.
std::vector<gemmi::Fractional> shifts;
for (size_t a = 0; a < copies; a++) {
const gemmi::Fractional nearest(std::round(body[a].c0.x - cm.x), std::round(body[a].c0.y - cm.y),
std::round(body[a].c0.z - cm.z));
for (int dx = -1; dx <= 1; dx++)
for (int dy = -1; dy <= 1; dy++)
for (int dz = -1; dz <= 1; dz++) {
const gemmi::Fractional shift(nearest.x + dx, nearest.y + dy, nearest.z + dz);
if (cell.orthogonalize_difference(cm + shift - body[a].c0).length()
<= body[a].radius + body[b].radius + CLASH_A
&& std::none_of(shifts.begin(), shifts.end(), [&](const gemmi::Fractional &o) {
return o.x == shift.x && o.y == shift.y && o.z == shift.z;
}))
shifts.push_back(shift);
}
}
for (const gemmi::Fractional &shift : shifts) {
const bool unmoved = op == gemmi::Op::identity() && shift.x == 0 && shift.y == 0 && shift.z == 0;
for (const gemmi::Fractional &m : mate) {
const gemmi::Position x = cell.orthogonalize(m + shift);
images.push_back({cube(x, 0, 0, 0), x, static_cast<int>(b), unmoved});
}
}
}
std::sort(images.begin(), images.end(), [](const Image &p, const Image &q) { return p.key < q.key; });
size_t clashing = 0;
for (size_t a = 0; a < copies; a++)
for (const gemmi::Fractional &f : body[a].ca) {
const gemmi::Position x = cell.orthogonalize(f);
bool clash = false;
for (int ex = -1; ex <= 1 && !clash; ex++)
for (int ey = -1; ey <= 1 && !clash; ey++)
for (int ez = -1; ez <= 1 && !clash; ez++) {
const long k = cube(x, ex, ey, ez);
auto it = std::lower_bound(images.begin(), images.end(), k,
[](const Image &p, long v) { return p.key < v; });
for (; it != images.end() && it->key == k; ++it)
if (!(it->unmoved && it->copy == static_cast<int>(a)) && x.dist_sq(it->x) < CLASH_A * CLASH_A) {
clash = true;
break;
}
}
clashing += clash;
}
return static_cast<double>(clashing) / static_cast<double>(total);
}
gemmi::Structure AssembleCopies(const std::vector<gemmi::Structure> &copies) {
static const std::string names = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789";
gemmi::Structure out = copies.at(0);
out.models[0].chains.clear();
size_t n = 0;
for (const gemmi::Structure &c : copies)
for (gemmi::Chain ch : c.models[0].chains) {
ch.name = names[n % names.size()] + (n < names.size() ? std::string() : std::to_string(n / names.size()));
out.models[0].chains.push_back(std::move(ch));
n++;
}
return out;
}
std::vector<gemmi::Structure> SplitCopies(const gemmi::Structure &assembly, int copies) {
std::vector<gemmi::Structure> body(copies, assembly);
const size_t chains = assembly.models[0].chains.size();
for (int k = 0; k < copies; k++) {
body[k].models[0].chains.clear();
for (size_t i = 0; i < chains; i++)
if (static_cast<int>(i * copies / chains) == k)
body[k].models[0].chains.push_back(assembly.models[0].chains[i]);
}
return body;
}
std::optional<SubgroupSetting> SubgroupBasis(const gemmi::UnitCell &g, const gemmi::SpaceGroup &g_sg,
const gemmi::UnitCell &h, const gemmi::SpaceGroup &h_sg) {
constexpr double TOLERANCE = 0.01; // relative, on lengths, and on the cosines of the angles
const gemmi::GroupOps gops = g_sg.operations();
const auto lattice = [&](const gemmi::Vec3 &v) {
for (const gemmi::Op::Tran &c : gops.cen_ops) {
const double d[3] = {v.x - static_cast<double>(c[0]) / gemmi::Op::DEN,
v.y - static_cast<double>(c[1]) / gemmi::Op::DEN,
v.z - static_cast<double>(c[2]) / gemmi::Op::DEN};
if (std::all_of(d, d + 3, [](double x) { return std::fabs(x - std::round(x)) < 1e-9; }))
return true;
}
return false;
};
const double length[3] = {h.a, h.b, h.c};
// Lattice vectors of g, in g's fractional coordinates, as long as each axis of h.
std::array<std::vector<gemmi::Vec3>, 3> axes;
for (int i = -4; i <= 4; i++)
for (int j = -4; j <= 4; j++)
for (int k = -4; k <= 4; k++) {
const gemmi::Vec3 v(i / 2.0, j / 2.0, k / 2.0);
if (!lattice(v))
continue;
const double len = g.orthogonalize_difference(gemmi::Fractional(v)).length();
for (int a = 0; a < 3; a++)
if (std::fabs(len - length[a]) < TOLERANCE * length[a])
axes[a].push_back(v);
}
const double cosine[3] = {std::cos(h.alpha * PI / 180), std::cos(h.beta * PI / 180), std::cos(h.gamma * PI / 180)};
const auto cos_between = [&](const gemmi::Vec3 &u, const gemmi::Vec3 &v) {
const gemmi::Position x = g.orthogonalize_difference(gemmi::Fractional(u));
const gemmi::Position y = g.orthogonalize_difference(gemmi::Fractional(v));
return x.dot(y) / (x.length() * y.length());
};
// The origin of h's setting in it: g's operators carried into h's axes, those of them h has shifted by o
// onto h's own - each translation t' + (1 - R) o the same as h's, modulo the lattice.
const gemmi::GroupOps hops = h_sg.operations();
const auto origin = [&](const gemmi::Mat33 &m) -> std::optional<gemmi::Vec3> {
const gemmi::Mat33 inv = m.inverse();
struct Pair { gemmi::Mat33 r; gemmi::Vec3 t, t_h; };
std::vector<Pair> pairs;
for (const gemmi::Op &op : gops.all_ops_sorted()) {
gemmi::Mat33 r;
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
r.a[i][j] = static_cast<double>(op.rot[i][j]) / gemmi::Op::DEN;
const gemmi::Mat33 rh = inv.multiply(r).multiply(m);
const gemmi::Vec3 th = inv.multiply(gemmi::Vec3(op.tran[0], op.tran[1], op.tran[2]) / gemmi::Op::DEN);
for (const gemmi::Op &hop : hops.all_ops_sorted()) {
bool same = true;
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
same = same && std::fabs(rh.a[i][j] - static_cast<double>(hop.rot[i][j]) / gemmi::Op::DEN) < 1e-6;
if (same)
pairs.push_back({rh, th, gemmi::Vec3(hop.tran[0], hop.tran[1], hop.tran[2]) / gemmi::Op::DEN});
}
}
for (int i = 0; i < 24; i++)
for (int j = 0; j < 24; j++)
for (int k = 0; k < 24; k++) {
const gemmi::Vec3 o(i / 24.0, j / 24.0, k / 24.0);
const bool all = std::all_of(pairs.begin(), pairs.end(), [&](const Pair &p) {
// Each of g's operators that lands on h's rotations has to land on one of h's operators.
const gemmi::Vec3 t = p.t + o - p.r.multiply(o);
for (const gemmi::Op &hop : hops.all_ops_sorted()) {
gemmi::Mat33 hr;
for (int a = 0; a < 3; a++)
for (int b = 0; b < 3; b++)
hr.a[a][b] = static_cast<double>(hop.rot[a][b]) / gemmi::Op::DEN;
if (!hr.approx(p.r, 1e-6))
continue;
const gemmi::Vec3 d = t - gemmi::Vec3(hop.tran[0], hop.tran[1], hop.tran[2]) / gemmi::Op::DEN;
if (std::fabs(d.x - std::round(d.x)) < 1e-6 && std::fabs(d.y - std::round(d.y)) < 1e-6
&& std::fabs(d.z - std::round(d.z)) < 1e-6)
return true;
}
return false;
});
if (all)
return o;
}
return std::nullopt;
};
std::optional<SubgroupSetting> best;
double best_error = INFINITY;
for (const gemmi::Vec3 &a : axes[0])
for (const gemmi::Vec3 &b : axes[1])
for (const gemmi::Vec3 &c : axes[2]) {
const double error = std::fabs(cos_between(b, c) - cosine[0]) + std::fabs(cos_between(a, c) - cosine[1])
+ std::fabs(cos_between(a, b) - cosine[2]);
const gemmi::Mat33 m(a.x, b.x, c.x, a.y, b.y, c.y, a.z, b.z, c.z);
if (error < 3 * TOLERANCE && m.determinant() > 0 && error < best_error)
if (const auto o = origin(m)) {
best_error = error;
best = SubgroupSetting{m, *o};
}
}
return best;
}
gemmi::Structure ExpandIntoSubgroup(const gemmi::Structure &placed, int copies, const gemmi::SpaceGroup &g_sg,
const SubgroupSetting &setting, const gemmi::UnitCell &h,
const gemmi::SpaceGroup &h_sg, int *copies_out) {
constexpr double SAME = 0.01; // fractional, between the centroids of two images
const gemmi::Mat33 to_h = setting.basis.inverse();
const gemmi::GroupOps hops = h_sg.operations();
std::vector<gemmi::Structure> out;
std::vector<gemmi::Fractional> centroids;
for (const gemmi::Structure &copy : SplitCopies(placed, copies))
for (const gemmi::Op &op : g_sg.operations().all_ops_sorted()) {
gemmi::Structure image = copy;
std::vector<gemmi::Position> pos = ModelPositions(image.models[0]);
std::vector<gemmi::Vec3> f;
gemmi::Vec3 c;
for (const gemmi::Position &p : pos) {
const gemmi::Fractional x = placed.cell.fractionalize(p);
const std::array<double, 3> y = op.apply_to_xyz({{x.x, x.y, x.z}});
f.push_back(to_h.multiply(gemmi::Vec3(y[0], y[1], y[2])) + setting.origin);
c += f.back();
}
c *= 1.0 / static_cast<double>(f.size());
const gemmi::Vec3 shift(-std::floor(c.x), -std::floor(c.y), -std::floor(c.z));
const gemmi::Fractional centre(c + shift);
// An image the subgroup's own symmetry makes of one already kept is the same copy.
const bool seen = std::any_of(centroids.begin(), centroids.end(), [&](const gemmi::Fractional &k) {
for (const gemmi::Op &hop : hops.all_ops_sorted()) {
const std::array<double, 3> y = hop.apply_to_xyz({{k.x, k.y, k.z}});
const double d[3] = {y[0] - centre.x, y[1] - centre.y, y[2] - centre.z};
if (std::all_of(d, d + 3, [](double x) { return std::fabs(x - std::round(x)) < SAME; }))
return true;
}
return false;
});
if (seen)
continue;
centroids.push_back(centre);
for (size_t i = 0; i < pos.size(); i++)
pos[i] = h.orthogonalize(gemmi::Fractional(f[i] + shift));
SetModelPositions(image.models[0], pos);
out.push_back(std::move(image));
}
gemmi::Structure st = AssembleCopies(out);
st.cell = h;
st.spacegroup_hm = h_sg.xhm();
st.setup_cell_images();
if (copies_out != nullptr)
*copies_out = static_cast<int>(out.size());
return st;
}
// The density of the Matthews coefficient over protein crystals: a Gumbel (extreme-value) distribution
// - the family Kantardjieff & Rupp fit to this highly skewed peak - with their mode (2.34 A^3/Da) and
// mean (2.69 A^3/Da) for 10471 protein crystal forms, which make its scale (mean - mode) / Euler's
// constant; without their resolution dependence.
// Following Kantardjieff & Rupp (2003) Protein Sci. 12, 1865-1871
std::vector<double> MatthewsProbabilities(double asu_volume, double mass_da) {
constexpr double MODE = 2.34, MEAN = 2.69, EULER_GAMMA = 0.5772156649;
constexpr double SCALE = (MEAN - MODE) / EULER_GAMMA;
std::vector<double> p;
for (int n = 1; asu_volume / (n * mass_da) >= VM_MIN; n++) {
const double z = (asu_volume / (n * mass_da) - MODE) / SCALE;
p.push_back(std::exp(-z - std::exp(-z)));
}
const double sum = std::accumulate(p.begin(), p.end(), 0.0);
for (double &x : p)
x /= sum;
return p;
}
std::vector<gemmi::SpaceGroup> SameRotationGroups(const gemmi::SpaceGroup &sg) {
const gemmi::GroupOps ops = sg.operations();
std::vector<gemmi::SpaceGroup> out{sg};
for (const gemmi::SpaceGroup &g : gemmi::spacegroup_tables::main) {
const gemmi::GroupOps gops = g.operations();
if (!same_rotations_and_centring(ops, gops) || gemmi::find_spacegroup_by_ops(gops) != &g)
continue;
// The table lists some groups a second time on another origin (P 21 21 2 and "P 21212(a)",
// I 2 2 2 and "I 2 2 2a"): the same crystal, which would tie with itself in the choice.
if (std::none_of(out.begin(), out.end(), [&](const gemmi::SpaceGroup &k) { return same_up_to_origin(k, g); }))
out.push_back(g);
}
return out;
}
// --- the search and the decision ---
namespace {
// The search and the decision for one search model.
MolecularReplacementResult search_model(const SearchModel &model, const std::vector<MRMerge> &merges, Logger &logger,
size_t nthreads) {
MolecularReplacementResult result;
result.residues_read = model.residues_read;
result.residues_kept = model.residues_kept;
result.mass_da = model.mass_da;
result.search_model = model.description;
auto t0 = std::chrono::steady_clock::now();
if (merges.empty() || merges[0].merged == nullptr) {
result.failure_reason = "no merge to search";
return result;
}
double data_d_min = 1e9;
for (const MergedReflection &r : *merges[0].merged)
if (r.d > 0)
data_d_min = std::min<double>(data_d_min, r.d);
result.d_min_search = std::max(SEARCH_D_MIN, data_d_min);
result.d_max_search = SEARCH_D_MAX;
result.d_min_score = std::max(SCORE_D_MIN, data_d_min);
logger.Info("Model prior: {}: {} of {} residues kept at pLDDT >= 70 ({:.1f} kDa, radius {:.1f} A)",
model.description, model.residues_kept, model.residues_read, model.mass_da / 1000.0, model.radius_A);
t0 = std::chrono::steady_clock::now();
std::optional<MolecularTransform> ft_made;
try {
ft_made.emplace(model, result.d_min_search, nthreads);
} catch (const std::exception &e) { // the GPU's FFT
result.failure_reason = std::string("the molecular transform: ") + e.what();
return result;
}
const MolecularTransform &ft = *ft_made;
#ifdef JFJOCH_USE_CUDA
// The rotation and translation functions on the GPU where there is one; the rest of the search is
// the same either way.
std::unique_ptr<MolecularReplacementGPU> gpu;
if (get_gpu_count() > 0) {
try {
gpu = std::make_unique<MolecularReplacementGPU>(ft);
} catch (const JFJochException &e) {
result.failure_reason = std::string("GPU: ") + e.what();
return result;
}
}
result.gpu_search = gpu != nullptr;
#endif
result.t_transform = seconds_since(t0);
result.tncs = DetectTranslationalNcs(*merges[0].merged, gemmi::UnitCell(merges[0].cell), merges[0].space_group);
logger.Info("Model prior: native Patterson's highest peak off the lattice {:.2f} of the origin at {:.3f} {:.3f} {:.3f}{}",
result.tncs.height, result.tncs.vector.x, result.tncs.vector.y, result.tncs.vector.z,
result.tncs.present ? " - tNCS" : "");
const std::vector<gemmi::Mat33> sphere = UniformRotations(ROTATION_STEP_DEG);
ScoreEngines scorers;
std::vector<gemmi::Mat33> orientations; // of each candidate's placement, before the rigid body
// One merge: the rotation function over its point group, then for each of its groups the
// translation function and the best placements through packing, the rigid body and R-free.
const auto search_merge = [&](size_t m) {
const MRMerge &merge = merges[m];
const gemmi::UnitCell cell(merge.cell);
const double d_min = result.d_min_score;
// --- rotation ---
t0 = std::chrono::steady_clock::now();
const SearchReflections refl = MakeSearchReflections(*merge.merged, cell, merge.space_group,
result.d_min_search, result.d_max_search);
const std::vector<gemmi::Mat33> pg = CartesianPointGroup(merge.space_group, cell);
const std::vector<gemmi::Mat33> rotations = RotationAsu(sphere, pg);
// The likelihood's reflections: down to LLG_D_MIN, or as far as the lowest-resolution LLG_TERMS of
// them reach where there are more - at least the search's own.
std::vector<double> ds;
for (const MergedReflection &r : *merge.merged)
if (std::isfinite(r.I) && r.d <= result.d_max_search)
ds.push_back(r.d);
std::sort(ds.begin(), ds.end(), std::greater<>());
double d_llg = std::max(LLG_D_MIN, data_d_min);
if (ds.size() > LLG_TERMS)
d_llg = std::max(d_llg, std::min(result.d_min_search, ds[LLG_TERMS - 1]));
const SearchReflections lrefl = MakeSearchReflections(*merge.merged, cell, merge.space_group, d_llg,
result.d_max_search);
#ifdef JFJOCH_USE_CUDA
const std::vector<double> rf = gpu ? gpu->RotationFunction(refl, cell, merge.space_group, rotations)
: RotationFunction(ft, refl, cell, merge.space_group, rotations, nthreads);
#else
const std::vector<double> rf = RotationFunction(ft, refl, cell, merge.space_group, rotations, nthreads);
#endif
double rf_mean = 0, rf_sd = 0;
for (double x : rf)
rf_mean += x;
rf_mean /= static_cast<double>(rf.size());
for (double x : rf)
rf_sd += (x - rf_mean) * (x - rf_mean);
rf_sd = std::sqrt(rf_sd / static_cast<double>(rf.size()));
// The highest orientations, one per peak: a rotation within two steps of one already taken
// (in any of its symmetry classes) belongs to that one's peak.
std::vector<size_t> order(rf.size());
std::iota(order.begin(), order.end(), 0);
std::stable_sort(order.begin(), order.end(), [&](size_t a, size_t b) { return rf[a] > rf[b]; });
std::vector<size_t> peaks;
for (size_t i : order) {
if (std::none_of(peaks.begin(), peaks.end(), [&](size_t p) {
return RotationDistanceDeg(rotations[i], rotations[p], pg) < 2 * ROTATION_STEP_DEG;
}))
peaks.push_back(i);
if (static_cast<int>(peaks.size()) == TRANSLATION_ROTATIONS)
break;
}
// Each peak refined to the function's local maximum, off the sampling: a large model turned by
// half a sampling step is a few angstroms out at its surface, which costs the translation
// function much of its peak. Hill climbing about the three axes, the step halved from half the
// sampling to an eighth.
std::vector<gemmi::Mat33> peak_rotations;
std::vector<double> peak_rf;
for (size_t p : peaks) {
peak_rotations.push_back(rotations[p]);
peak_rf.push_back(rf[p]);
}
for (double step = ROTATION_STEP_DEG / 2; step >= ROTATION_STEP_DEG / 8 - 1e-9; step /= 2) {
std::vector<char> active(peaks.size(), 1);
for (int iter = 0; iter < PEAK_REFINE_ITERATIONS
&& std::find(active.begin(), active.end(), 1) != active.end(); iter++) {
std::vector<gemmi::Mat33> trial;
for (size_t k = 0; k < peaks.size(); k++)
if (active[k])
for (int axis = 0; axis < 3; axis++)
for (double sign : {1.0, -1.0})
trial.push_back(axis_rotation(axis, sign * step).multiply(peak_rotations[k]));
#ifdef JFJOCH_USE_CUDA
const std::vector<double> v = gpu ? gpu->RotationFunction(refl, cell, merge.space_group, trial)
: RotationFunction(ft, refl, cell, merge.space_group, trial, nthreads);
#else
const std::vector<double> v = RotationFunction(ft, refl, cell, merge.space_group, trial, nthreads);
#endif
size_t next = 0;
for (size_t k = 0; k < peaks.size(); k++) {
if (!active[k])
continue;
size_t best = next;
for (size_t j = next; j < next + 6; j++)
if (v[j] > v[best])
best = j;
if (v[best] > peak_rf[k]) {
peak_rf[k] = v[best];
peak_rotations[k] = trial[best];
} else {
active[k] = 0;
}
next += 6;
}
}
}
result.t_rotation += seconds_since(t0);
logger.Info("Model prior: {} in {}: {} reflections {:.1f}-{:.1f} A, {} orientations, best RF Z {:.1f}",
merge.name, merge.space_group.xhm(), refl.hkl.size(), result.d_max_search, result.d_min_search,
rotations.size(), rf_sd > 0 ? (peak_rf.front() - rf_mean) / rf_sd : 0.0);
// --- translation, for every group and every peak orientation ---
t0 = std::chrono::steady_clock::now();
const std::vector<gemmi::SpaceGroup> groups = SameRotationGroups(merge.space_group);
std::vector<TranslationPeak> tf_all;
for (const gemmi::SpaceGroup &sg : groups) {
const TranslationSetup setup = MakeTranslationSetup(refl, cell, sg);
std::vector<TranslationPeak> tf(peaks.size());
#ifdef JFJOCH_USE_CUDA
if (gpu)
tf = gpu->TranslationFunction(setup, peak_rotations);
else
#endif
ParallelFor(static_cast<int>(peaks.size()), nthreads, [&](int i) {
tf[i] = TranslationFunction(ft, setup, peak_rotations[i]);
});
tf_all.insert(tf_all.end(), tf.begin(), tf.end());
}
result.t_translation += seconds_since(t0);
// Each group searched on a thread of its own - they share nothing but the read-only search and the
// score engines - its log kept and replayed in order.
struct GroupResult {
MRCandidate cand;
gemmi::Mat33 orientation;
std::vector<gemmi::Mat33> copy_rotation; // of every copy placed
Logger log = Logger::Buffered();
};
const auto search_group = [&](size_t g, const std::vector<gemmi::Mat33> &seeds) {
GroupResult out;
Logger &glog = out.log;
const gemmi::SpaceGroup &sg = groups[g];
MRCandidate cand;
cand.merge = static_cast<int>(m);
cand.space_group = sg;
const gemmi::GroupOps gops = sg.operations();
const double asu_volume = cell.volume / static_cast<double>(gops.order());
cand.vm = asu_volume / model.mass_da;
const std::vector<double> pm = MatthewsProbabilities(asu_volume, model.mass_da);
for (size_t n = 0; n < pm.size(); n++) {
if (pm[n] > pm[cand.copies_likely - 1])
cand.copies_likely = static_cast<int>(n) + 1;
if (pm[n] >= MATTHEWS_P_MIN)
cand.copies_max = static_cast<int>(n) + 1;
}
// The orientations searched: the rotation function's peaks and the seeds, which are tried first.
std::vector<gemmi::Mat33> rots = peak_rotations;
std::vector<double> rfs = peak_rf;
std::vector<TranslationPeak> tfs(tf_all.begin() + g * peaks.size(), tf_all.begin() + (g + 1) * peaks.size());
std::vector<size_t> by_tf(peaks.size());
std::iota(by_tf.begin(), by_tf.end(), 0);
std::stable_sort(by_tf.begin(), by_tf.end(), [&](size_t a, size_t b) { return tfs[a].z > tfs[b].z; });
if (!seeds.empty()) {
const TranslationSetup setup = MakeTranslationSetup(refl, cell, sg);
std::vector<TranslationPeak> tf_seed(seeds.size());
#ifdef JFJOCH_USE_CUDA
const std::vector<double> rf_seed = gpu ? gpu->RotationFunction(refl, cell, merge.space_group, seeds)
: RotationFunction(ft, refl, cell, merge.space_group, seeds, nthreads);
if (gpu)
tf_seed = gpu->TranslationFunction(setup, seeds);
else
#else
const std::vector<double> rf_seed = RotationFunction(ft, refl, cell, merge.space_group, seeds, nthreads);
#endif
for (size_t k = 0; k < seeds.size(); k++)
tf_seed[k] = TranslationFunction(ft, setup, seeds[k]);
std::vector<size_t> first;
for (size_t k = 0; k < seeds.size(); k++) {
first.push_back(rots.size());
rots.push_back(seeds[k]);
rfs.push_back(rf_seed[k]);
tfs.push_back(tf_seed[k]);
}
by_tf.insert(by_tf.begin(), first.begin(), first.end());
}
const TranslationPeak *tf = tfs.data();
const LikelihoodSetup lsetup = MakeLikelihoodSetup(lrefl, cell, sg);
const std::vector<std::complex<double>> no_known(lsetup.Terms(), 0.0);
// The further copies' peaks are ranked by the likelihood on the search's own reflections, read off
// the molecular transform.
const LikelihoodSetup rsetup = MakeLikelihoodSetup(refl, cell, sg);
// --- the first copy: the best few through packing, the rigid body and R-free ---
RigidBodyEngines engines;
bool have = false;
gemmi::Mat33 orientation;
// Where not even one copy fits the cell, only the best translation-function placement is
// checked for packing.
const bool fits = cand.vm >= VM_MIN;
// Each copy placed: its orientation and its centroid, where the rigid body left them.
std::vector<gemmi::Structure> copies;
std::vector<gemmi::Mat33> copy_rotation;
std::vector<gemmi::Fractional> copy_position;
// CopyFactors of each, for the likelihood and for ranking (rsetup), as they are needed.
std::vector<std::vector<std::complex<double>>> copy_g, copy_rg;
// CopyFactors of the orientations searched, as they are needed.
std::vector<std::vector<std::complex<double>>> rot_g(rots.size());
// Down the list until RESCORED_PER_GROUP placements that pack have been refined: a placement
// clashing with its own mates over twice the limit is not refined - the rigid body moves a
// search placement by a few degrees and a fraction of an angstrom, and the clash fraction moves
// by a percent or two with it - and does not count, so that a group is CLASH only where no
// orientation the rotation function offers packs in it.
int refined = 0;
for (size_t k = 0; k < by_tf.size() && refined < (fits ? RESCORED_PER_GROUP : 1); k++) {
const size_t i = by_tf[k];
auto placed = std::make_shared<gemmi::Structure>(
PlaceSearchModel(model, rots[i], tf[i].t, cell, sg));
double clash = ClashFraction(*placed);
PlacementScore score{NAN, NAN, NAN, NAN};
gemmi::Mat33 moved;
if (fits && clash <= 2 * CLASH_MAX) {
if (engines.get() == nullptr)
engines = make_engines(*placed, sg, *merge.merged, d_min, 1, logger);
score = refine_placement(*placed, sg, *merge.merged, d_min, engines.get(), nthreads, &scorers,
&moved);
clash = ClashFraction(*placed);
refined++;
} else if (fits && have) {
continue;
} else if (!fits) {
refined++;
}
// The refined placement's likelihood, judged against the same orientation elsewhere.
// The refined placement with the lowest R-free; the best translation-function one where
// none was refined.
if (!have || (std::isfinite(score.r_free) && !(cand.r_free <= score.r_free))) {
have = true;
cand.rf_z = rf_sd > 0 ? (rfs[i] - rf_mean) / rf_sd : 0.0;
cand.tf_z = tf[i].z;
cand.copy_tf_z = {tf[i].z};
cand.clash_fraction = clash;
cand.r_work = score.r_work;
cand.r_free = score.r_free;
cand.k_sol = score.k_sol;
cand.b_sol = score.b_sol;
cand.placed = placed;
cand.copies = 1;
orientation = rots[i];
copies = {*placed};
copy_rotation = {moved.multiply(rots[i])};
copy_position = {centroid(*placed)};
copy_g = {{}};
copy_rg = {{}};
}
glog.Info("Model prior: {:<12} RF Z {:5.1f} TF Z {:5.1f} clash {:4.1f}% R-work {:.3f} R-free {:.3f}",
sg.xhm(), rf_sd > 0 ? (rfs[i] - rf_mean) / rf_sd : 0.0, tf[i].z, 100 * clash,
score.r_work, score.r_free);
}
// The likelihood of the placement kept, judged against the same orientation elsewhere - only where
// its translation function's Z-score does not already stand out.
if (have && std::isfinite(cand.r_free) && cand.tf_z < TF_Z_PLACED) {
copy_g[0] = CopyFactors(model, lsetup, cell, copy_rotation[0], d_llg, nthreads, &ft);
const LikelihoodPeak lk = LikelihoodTranslationZ(lsetup, no_known, 0.0, copy_g[0], copy_position[0],
model.error_A, nthreads);
cand.llg = lk.llg;
cand.llg_z = lk.z;
cand.copy_llg_z = {lk.z};
glog.Info("Model prior: {:<12} the placement kept: LLG {:.1f}, Z {:.1f} over its orientation elsewhere",
sg.xhm(), lk.llg, lk.z);
}
// --- further copies, one at a time, with the copies placed as a known partial structure ---
// Each copy is refined on its own and all of them scored together.
std::optional<TranslationSetup> setup2;
// The first `count` copies refined in turn, each with all the others fixed, and the assembly scored.
const auto refine_beside_others = [&](int count) {
for (int c = 0; c < count; c++) {
std::vector<gemmi::Structure> rest;
for (size_t o = 0; o < copies.size(); o++)
if (static_cast<int>(o) != c)
rest.push_back(copies[o]);
const gemmi::Structure fixed = AssembleCopies(rest);
gemmi::Mat33 moved;
refine_rigid_body(copies[c], sg, *merge.merged, d_min, engines.get(), nthreads, &moved,
&fixed.models[0]);
copy_rotation[c] = moved.multiply(copy_rotation[c]);
copy_position[c] = centroid(copies[c]);
copy_g[c].clear();
copy_rg[c].clear();
}
const gemmi::Structure assembly = AssembleCopies(copies);
const PlacementScore held{NAN, NAN, cand.k_sol, cand.b_sol};
const PlacementScore score = score_placement(assembly, sg, *merge.merged, d_min, nthreads, &scorers,
&held);
cand.clash_fraction = ClashFraction(assembly, cand.copies);
cand.r_work = score.r_work;
cand.r_free = score.r_free;
cand.k_sol = score.k_sol;
cand.b_sol = score.b_sol;
cand.placed = std::make_shared<gemmi::Structure>(assembly);
glog.Info("Model prior: {:<12} {} copies, refined beside each other: R-work {:.3f} R-free {:.3f}",
sg.xhm(), cand.copies, score.r_work, score.r_free);
};
while (fits && cand.copies < cand.copies_max && cand.clash_fraction <= CLASH_MAX
&& std::isfinite(cand.r_free)) {
// The copies placed are searched for again in their own orientations: copies related by a
// translation alone (tNCS) share one.
std::vector<gemmi::Mat33> rots2 = rots;
rots2.insert(rots2.end(), copy_rotation.begin(), copy_rotation.end());
if (!setup2)
setup2 = MakeTranslationSetup(refl, cell, sg, true);
const std::vector<std::complex<double>> known = KnownStructureFactors(ft, *setup2, cell, copy_rotation,
copy_position);
std::vector<TranslationPeak> tf2(rots2.size());
#ifdef JFJOCH_USE_CUDA
if (gpu)
tf2 = gpu->TranslationFunction(*setup2, rots2, &known);
else
#endif
ParallelFor(static_cast<int>(rots2.size()), nthreads, [&](int i) {
tf2[i] = TranslationFunction(ft, *setup2, rots2[i], &known);
});
std::vector<size_t> by_tf2(rots2.size());
std::iota(by_tf2.begin(), by_tf2.end(), 0);
std::stable_sort(by_tf2.begin(), by_tf2.end(), [&](size_t a, size_t b) { return tf2[a].z > tf2[b].z; });
struct Trial {
size_t i = 0;
double clash = 1;
double gain = 0;
LikelihoodPeak refined;
std::vector<std::complex<double>> g;
PlacementScore score{NAN, NAN, NAN, NAN};
gemmi::Structure copy, assembly;
gemmi::Mat33 rotation;
};
std::optional<Trial> best;
const int n = cand.copies + 1;
// The highest peaks that do not sit on a copy already placed (or a mate of one), ranked by the
// likelihood gain of the copy beside the known ones, and the best few of those through the
// rigid body. A copy that only repeats a known one - moved by one of the group's origin shifts,
// it doubles F on half the reflections and cancels it on the other half - gains nothing.
std::vector<std::complex<double>> rknown(rsetup.Terms(), 0.0);
for (size_t c = 0; c < copies.size(); c++) {
if (copy_rg[c].empty())
copy_rg[c] = CopyFactors(model, rsetup, cell, copy_rotation[c], result.d_min_search, nthreads, &ft);
AddCopy(rsetup, copy_rg[c], copy_position[c], rknown);
}
double known_scale = 0;
const double rknown_llg = RefinedLogLikelihoodGain(rsetup, rknown, model.error_A, &known_scale);
// Down the highest peaks until COPY_TRIALS that gain likelihood - a copy that lowers it is not
// one, and a peak on a copy already placed (or a mate of one) does - and pack. The gain that
// ranks them is taken at the known copies' sigma_A.
std::vector<Trial> trials;
for (size_t k = 0; k < by_tf2.size() && k < 2 * COPY_TRIALS && trials.size() < COPY_TRIALS; k++) {
Trial trial;
trial.i = by_tf2[k];
std::vector<std::complex<double>> &g = trial.i < rots.size() ? rot_g[trial.i]
: copy_rg[trial.i - rots.size()];
if (g.empty())
g = CopyFactors(model, rsetup, cell, rots2[trial.i], result.d_min_search, nthreads, &ft);
std::vector<std::complex<double>> f = rknown;
AddCopy(rsetup, g, tf2[trial.i].t, f);
trial.gain = LogLikelihoodGain(rsetup, f, model.error_A, known_scale) - rknown_llg;
if (!(trial.gain > 0))
continue;
trial.copy = PlaceSearchModel(model, rots2[trial.i], tf2[trial.i].t, cell, sg);
std::vector<gemmi::Structure> all = copies;
all.push_back(trial.copy);
trial.clash = ClashFraction(AssembleCopies(all), n);
if (trial.clash > 2 * CLASH_MAX)
continue;
trials.push_back(std::move(trial));
}
std::stable_sort(trials.begin(), trials.end(),
[](const Trial &a, const Trial &b) { return a.gain > b.gain; });
if (trials.size() > RESCORED_PER_GROUP)
trials.resize(RESCORED_PER_GROUP);
const gemmi::Structure others = AssembleCopies(copies);
// The known copies on the likelihood's own reflections, made where a trial needs its Z-score.
std::optional<std::vector<std::complex<double>>> lknown;
double known_llg = 0;
for (Trial &trial : trials) {
std::vector<gemmi::Structure> all = copies;
gemmi::Mat33 moved;
refine_rigid_body(trial.copy, sg, *merge.merged, d_min, engines.get(), nthreads, &moved,
&others.models[0]);
all.push_back(trial.copy);
trial.assembly = AssembleCopies(all);
trial.clash = ClashFraction(trial.assembly, n);
const PlacementScore held{NAN, NAN, cand.k_sol, cand.b_sol};
trial.score = score_placement(trial.assembly, sg, *merge.merged, d_min, nthreads, &scorers, &held);
trial.rotation = moved.multiply(rots2[trial.i]);
trial.refined.gain = trial.refined.llg = trial.refined.z = NAN;
if (tf2[trial.i].z < TF_Z_PLACED) {
if (!lknown) {
lknown.emplace(lsetup.Terms(), 0.0);
for (size_t c = 0; c < copies.size(); c++) {
if (copy_g[c].empty())
copy_g[c] = CopyFactors(model, lsetup, cell, copy_rotation[c], d_llg, nthreads, &ft);
AddCopy(lsetup, copy_g[c], copy_position[c], *lknown);
}
known_llg = RefinedLogLikelihoodGain(lsetup, *lknown, model.error_A);
}
trial.g = CopyFactors(model, lsetup, cell, trial.rotation, d_llg, nthreads, &ft);
trial.refined = LikelihoodTranslationZ(lsetup, *lknown, known_llg, trial.g, centroid(trial.copy),
model.error_A, nthreads);
}
glog.Info("Model prior: {:<12} copy {} TF Z {:5.1f} LLG gain {:7.1f}, refined {:7.1f} Z {:5.1f} clash {:4.1f}% R-work {:.3f} R-free {:.3f}",
sg.xhm(), n, tf2[trial.i].z, trial.gain, trial.refined.gain,
trial.refined.z, 100 * trial.clash, trial.score.r_work, trial.score.r_free);
// The first, in that order, that stands out, packs and lowers R-free.
if ((tf2[trial.i].z >= TF_Z_PLACED || trial.refined.z >= TF_Z_PLACED) && trial.clash <= CLASH_MAX
&& trial.score.r_free <= cand.r_free - COPY_R_FREE_GAIN) {
best = trial;
break;
}
}
if (!best)
break;
copies.push_back(best->copy);
copy_rotation.push_back(best->rotation);
copy_position.push_back(centroid(best->copy));
copy_g.push_back(best->g);
copy_rg.push_back({});
cand.copies = n;
cand.copy_tf_z.push_back(tf2[best->i].z);
cand.tf_z = std::max(cand.tf_z, tf2[best->i].z);
cand.copy_llg_z.push_back(best->refined.z);
cand.llg_z = std::fmax(cand.llg_z, best->refined.z);
cand.llg = best->refined.llg; // NaN where the translation function stood out
cand.clash_fraction = best->clash;
cand.r_work = best->score.r_work;
cand.r_free = best->score.r_free;
cand.k_sol = best->score.k_sol;
cand.b_sol = best->score.b_sol;
cand.placed = std::make_shared<gemmi::Structure>(best->assembly);
// The first copy, refined alone, refined again beside the second: next to a copy related to it
// by a translation alone, one copy by itself fits the data no better than at random.
if (n == 2)
refine_beside_others(1);
}
// Every copy but the last refined again beside all the others, which each was refined beside
// only some of.
if (cand.copies >= 3)
refine_beside_others(cand.copies - 1);
// The copies' R-factors were taken with the bulk solvent of the first copy alone held; all of
// them scored with the solvent fitted to them.
if (cand.copies >= 2) {
const PlacementScore score = score_placement(*cand.placed, sg, *merge.merged, d_min, nthreads, &scorers);
cand.r_work = score.r_work;
cand.r_free = score.r_free;
cand.k_sol = score.k_sol;
cand.b_sol = score.b_sol;
}
if (cand.vm < VM_MIN)
cand.verdict = "MATTHEWS";
else if (cand.clash_fraction > CLASH_MAX)
cand.verdict = "CLASH";
// The likelihood's Z-score stands in for the translation function's only where the copies placed
// are as many as the cell most probably holds: beside content the search has not found, a copy's
// likelihood over random positions of its own is no measure of the placement.
else if (!(cand.tf_z >= TF_Z_PLACED || (cand.llg_z >= TF_Z_PLACED && cand.copies >= cand.copies_likely)))
cand.verdict = "NOT_PLACED";
else
cand.verdict = "PLACED";
out.copy_rotation = copy_rotation;
out.cand = std::move(cand);
out.orientation = orientation;
return out;
};
// Their time in wall-clock seconds, counted as rescoring (the further copies' translation functions
// included).
const auto t_groups = std::chrono::steady_clock::now();
std::vector<std::future<GroupResult>> running;
for (size_t g = 0; g < groups.size(); g++)
running.push_back(std::async(std::launch::async, search_group, g, std::vector<gemmi::Mat33>()));
std::vector<GroupResult> found;
for (auto &r : running)
found.push_back(r.get());
// With tNCS, every group that placed fewer copies than the best one searched again with the best one's
// copies' orientations tried first: the groups share the rotations, so a copy the search placed in one is in
// the same orientation in the others, and one group's search failing - next to copies related by a
// translation, the first copy alone fits no better than at random, and the copy the search keeps
// can be the wrong one - is then not taken for evidence against it. The better of the two searches
// is kept.
int lead = -1;
for (size_t g = 0; g < found.size(); g++)
if (found[g].cand.verdict == "PLACED" && (lead < 0 || found[g].cand.r_free < found[lead].cand.r_free))
lead = static_cast<int>(g);
if (lead >= 0 && result.tncs.present) {
std::vector<std::future<GroupResult>> again(found.size());
for (size_t g = 0; g < found.size(); g++)
if (static_cast<int>(g) != lead && found[g].cand.vm >= VM_MIN
&& found[g].cand.copies < found[lead].cand.copies)
again[g] = std::async(std::launch::async, search_group, g, found[lead].copy_rotation);
for (size_t g = 0; g < found.size(); g++)
if (again[g].valid()) {
GroupResult r = again[g].get();
found[g].log.Info("Model prior: {:<12} searched again from {}'s copies: {} copies, R-free {:.3f}, was {:.3f}",
groups[g].xhm(), groups[lead].xhm(), r.cand.copies, r.cand.r_free,
found[g].cand.r_free);
if (std::isfinite(r.cand.r_free) && !(found[g].cand.r_free <= r.cand.r_free)) {
r.log.ReplayInto(found[g].log);
r.log = found[g].log;
found[g] = std::move(r);
}
}
}
result.t_rescore += seconds_since(t_groups);
for (GroupResult &gr : found) {
gr.log.ReplayInto(logger);
orientations.push_back(gr.orientation);
result.candidates.push_back(std::move(gr.cand));
}
};
// --- the decision ---
// Within one merge the candidates are scored on the same reflections, so their R-free values can
// be compared directly: the placed candidate with the lowest R-free is chosen where it leads every
// other candidate of that merge that packs by R_FREE_MARGIN.
// The null is built for the candidate the choice would fall on, and only for it: a placed candidate
// that does not beat its own null is set aside, and the choice falls on the next.
const auto choose_in = [&](int merge) -> int {
int best = -1;
while (true) {
best = -1;
for (size_t i = 0; i < result.candidates.size(); i++) {
const MRCandidate &c = result.candidates[i];
if (c.merge == merge && c.verdict == "PLACED"
&& (best < 0 || c.r_free < result.candidates[best].r_free))
best = static_cast<int>(i);
}
if (best < 0)
return -1;
MRCandidate &c = result.candidates[best];
const auto t_null = std::chrono::steady_clock::now();
RigidBodyEngines engines = make_engines(*c.placed, c.space_group, *merges[merge].merged,
result.d_min_score, NULL_REPLICATES, logger);
// One copy: its replicates placed by the translation function, as it was.
std::function<std::vector<gemmi::Fractional>(const std::vector<gemmi::Mat33> &)> translate;
if (c.copies == 1)
translate = [&](const std::vector<gemmi::Mat33> &turned) {
const SearchReflections refl = MakeSearchReflections(*merges[merge].merged,
gemmi::UnitCell(merges[merge].cell),
merges[merge].space_group,
result.d_min_search, result.d_max_search);
const TranslationSetup setup = MakeTranslationSetup(refl, gemmi::UnitCell(merges[merge].cell),
c.space_group);
std::vector<TranslationPeak> tf(turned.size());
#ifdef JFJOCH_USE_CUDA
if (gpu)
tf = gpu->TranslationFunction(setup, turned);
else
#endif
ParallelFor(static_cast<int>(turned.size()), nthreads, [&](int i) {
tf[i] = TranslationFunction(ft, setup, turned[i]);
});
std::vector<gemmi::Fractional> out;
for (const TranslationPeak &p : tf)
out.push_back(p.t);
return out;
};
c.null_sigma = null_sigma(*c.placed, c.copies, orientations[best], c.space_group, *merges[merge].merged,
result.d_min_score, engines.get(), nthreads, &scorers, translate);
result.t_null += seconds_since(t_null);
if (c.null_sigma >= NULL_SIGMA_MIN)
break;
c.verdict = "NULL";
}
for (size_t i = 0; i < result.candidates.size(); i++) {
const MRCandidate &c = result.candidates[i];
if (static_cast<int>(i) == best || c.merge != merge || c.clash_fraction > CLASH_MAX)
continue;
// With tNCS one copy alone fits no better than at random, so a group whose search stopped at
// fewer copies than the chosen one's has not been shown unable to hold them - where the tNCS
// vector makes up the groups' screw translations, or where they are each other's enantiomorph.
if (result.tncs.present && c.copies < result.candidates[best].copies
&& (TncsRelatesGroups(c.space_group, result.candidates[best].space_group, result.tncs.vector)
|| enantiomorphs(c.space_group, result.candidates[best].space_group))) {
// Unless the chosen copies, carried into this group as they are and refined, build a crystal
// that fits clearly worse: where the tNCS vector does make up the difference between the groups,
// they build the same one.
const MRCandidate &b = result.candidates[best];
const PlacementScore carried = carry_copies(*b.placed, b.copies, c.space_group, *merges[merge].merged,
result.d_min_score, nthreads, &scorers, logger);
logger.Info("Model prior: {}'s {} copies carried into {}: R-free {:.3f} against {:.3f}",
b.space_group.xhm(), b.copies, c.space_group.xhm(), carried.r_free, b.r_free);
if (!(carried.r_free - b.r_free >= R_FREE_MARGIN)) {
result.decision = fmt::format("UNDECIDED: {} placed fewer copies than {} beside tNCS",
c.space_group.xhm(), b.space_group.xhm());
return -2;
}
}
if (!(c.r_free - result.candidates[best].r_free >= R_FREE_MARGIN)) {
result.decision = fmt::format("UNDECIDED: {} and {} within {:.2f} in R-free",
result.candidates[best].space_group.xhm(), c.space_group.xhm(),
R_FREE_MARGIN);
return -2;
}
}
return best;
};
// A merge whose every group the model cannot be in at all - one copy does not fit the cell, or
// every placement clashes with its own symmetry mates.
const auto cannot_hold = [&](int merge) {
bool any = false;
for (const MRCandidate &c : result.candidates)
if (c.merge == merge) {
any = true;
if (c.verdict != "MATTHEWS" && c.verdict != "CLASH")
return false;
}
return any;
};
// The twin hypothesis' merge is searched only where it can decide something: where the answer's
// merge cannot hold the model in any of its groups.
search_merge(0);
const int answer = choose_in(0);
if (answer >= 0) {
result.chosen = answer;
} else if (answer == -1 && merges.size() > 1 && cannot_hold(0)) {
search_merge(1);
const int twin = choose_in(1);
if (twin >= 0)
result.chosen = twin;
else if (twin == -1)
result.decision = "NOT_PLACED: no group of either merge holds a placed model";
} else if (answer == -1) {
result.decision = "NOT_PLACED: no group of the answer's merge holds a placed model";
}
// The twin comparison, where the answer's merge chose and the twin hypothesis' merge is there: on that
// merge, against its data as twinned by its law, the chosen placement with all its symmetry images - the
// group's crystal - and the model searched for afresh in the subgroup, as a crystal of the subgroup. A
// crystal of the group twinned by the law reads like one of the subgroup, and its placement read in the
// subgroup fits as well as it fits in the group; a crystal of the subgroup has copies the group's
// symmetry does not relate, which only a search in the subgroup places. The subgroup is chosen where its
// own placement fits the twinned data better by R_FREE_MARGIN in R-free - on reflections none of the
// copies, nor the fraction, was fitted to - and passes the same choice as any merge's (its null, and a
// lead over the subgroup's other groups).
if (result.chosen >= 0 && result.candidates[result.chosen].merge == 0 && merges.size() > 1
&& merges[1].twin_law) {
const MRCandidate g = result.candidates[result.chosen];
const MRMerge &tm = merges[1];
const gemmi::UnitCell h_cell(tm.cell);
if (const auto setting = SubgroupBasis(gemmi::UnitCell(merges[0].cell), g.space_group, h_cell, tm.space_group)) {
const gemmi::Op &law = *tm.twin_law;
int n = 0;
const gemmi::Structure expanded = ExpandIntoSubgroup(*g.placed, g.copies, g.space_group, *setting, h_cell,
tm.space_group, &n);
const PlacementScore as_group = score_placement(expanded, tm.space_group, *tm.merged, result.d_min_score,
nthreads, &scorers, nullptr, &law);
result.candidates[result.chosen].r_free_twin = as_group.r_free_twin;
const std::string decision = result.decision;
search_merge(1);
const int h = choose_in(1);
if (h >= 0) {
MRCandidate &c = result.candidates[h];
const PlacementScore t = score_placement(*c.placed, c.space_group, *tm.merged, result.d_min_score,
nthreads, &scorers, nullptr, &law);
c.r_work_twin = t.r_work_twin;
c.r_free_twin = t.r_free_twin;
c.twin_fraction = t.twin_fraction;
logger.Info("Model prior: twin comparison on the {} merge, twinned by {}: {} read as its {} copies "
"there R-free {:.3f}; searched for in {}, {} copies, R-free {:.3f} at a fraction {:.2f}",
tm.space_group.xhm(), law.as_hkl().triplet('h'), g.space_group.xhm(), n,
as_group.r_free_twin, c.space_group.xhm(), c.copies, t.r_free_twin, t.twin_fraction);
if (t.r_free_twin <= as_group.r_free_twin - R_FREE_MARGIN)
result.chosen = h;
}
result.decision = decision;
}
}
if (result.chosen >= 0) {
const MRCandidate &c = result.candidates[result.chosen];
const gemmi::SpaceGroup &data = merges[0].space_group;
if (c.merge == 0 && c.space_group.xhm() == data.xhm())
result.decision = "AGREES";
else if (c.merge == 0 && enantiomorphs(data, c.space_group))
result.decision = "ENANTIOMORPH";
else if (c.merge == 0)
result.decision = "SPACE_GROUP";
else
result.decision = "TWIN_SUBGROUP";
}
result.ok = true;
logger.Info("Model prior: {} (model {:.1f} s, transform {:.2f} s, rotation {:.2f} s, translation {:.2f} s, "
"rescoring {:.2f} s, null {:.2f} s)",
result.chosen >= 0 ? result.decision + " - " + result.candidates[result.chosen].space_group.xhm()
: result.decision,
result.t_model, result.t_transform, result.t_rotation, result.t_translation, result.t_rescore,
result.t_null);
return result;
}
} // namespace
MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector<MRMerge> &merges,
Logger &logger, size_t nthreads) {
const auto t0 = std::chrono::steady_clock::now();
SearchModel model = PrepareSearchModel(model_path);
if (!model.failure_reason.empty()) {
MolecularReplacementResult result;
result.model_path = model_path;
result.residues_read = model.residues_read;
result.failure_reason = model.failure_reason;
return result;
}
model.description = "the whole model";
const std::string pae_path = PaePathFor(model_path);
const std::vector<std::vector<int>> domains = pae_path.empty() ? std::vector<std::vector<int>>()
: PaeDomains(pae_path, model);
const double t_model = seconds_since(t0);
if (!pae_path.empty())
logger.Info("Model prior: {} rigid domain(s) in the predicted aligned error {}", domains.size(), pae_path);
MolecularReplacementResult result = search_model(model, merges, logger, nthreads);
result.model_path = model_path;
result.pae_domains = static_cast<int>(domains.size());
result.t_model += t_model;
// Where the whole model is placed nowhere and its prediction falls apart into domains whose relative
// placement it is not confident of, each domain is searched by itself; the domains that are placed have
// to agree on the space group, and the one with the lowest R-free is taken where that is the answer's
// own.
if (result.ok && result.chosen < 0 && result.decision.starts_with("NOT_PLACED") && domains.size() >= 2) {
const auto t_domains = std::chrono::steady_clock::now();
// Each domain on a thread of its own, its log kept and replayed in order.
std::vector<Logger> logs;
for (size_t i = 0; i < domains.size(); i++)
logs.push_back(Logger::Buffered());
std::vector<std::future<MolecularReplacementResult>> running;
for (size_t i = 0; i < domains.size(); i++)
running.push_back(std::async(std::launch::async, [&, i] {
return search_model(DomainModel(model, domains[i]), merges, logs[i], nthreads);
}));
std::vector<MolecularReplacementResult> found;
for (size_t i = 0; i < domains.size(); i++) {
MolecularReplacementResult r = running[i].get();
logs[i].ReplayInto(logger);
if (r.ok && (r.chosen >= 0 || r.decision.starts_with("UNDECIDED")))
found.push_back(std::move(r));
}
const double domain_seconds = seconds_since(t_domains);
int best = -1;
bool agree = true;
for (size_t i = 0; i < found.size(); i++) {
if (found[i].chosen < 0) {
agree = false;
continue;
}
if (best >= 0 && found[i].candidates[found[i].chosen].space_group.xhm()
!= found[best].candidates[found[best].chosen].space_group.xhm())
agree = false;
if (best < 0 || found[i].candidates[found[i].chosen].r_free < found[best].candidates[found[best].chosen].r_free)
best = static_cast<int>(i);
}
// A domain is part of the model: it may confirm the answer's group, not change it.
if (best >= 0 && agree && found[best].decision != "AGREES") {
result.decision = fmt::format("UNDECIDED: a PAE domain alone chooses {}",
found[best].candidates[found[best].chosen].space_group.xhm());
} else if (best >= 0 && agree) {
MolecularReplacementResult chosen = std::move(found[best]);
chosen.model_path = model_path;
chosen.pae_domains = result.pae_domains;
chosen.t_model = result.t_model;
chosen.t_transform = result.t_transform;
chosen.t_rotation = result.t_rotation;
chosen.t_translation = result.t_translation;
chosen.t_rescore = result.t_rescore;
chosen.t_null = result.t_null;
result = std::move(chosen);
result.t_domains = domain_seconds;
} else {
result.t_domains = domain_seconds;
if (!found.empty())
result.decision = "UNDECIDED: the PAE domains placed disagree";
}
logger.Info("Model prior: with the PAE domains searched by themselves: {}{}", result.decision,
result.chosen >= 0 ? " - " + result.candidates[result.chosen].space_group.xhm() + ", "
+ result.search_model : "");
}
return result;
}
void CheckEnantiomorphHand(MolecularReplacementResult &result, const std::vector<MergedReflection> &merged,
const gemmi::SpaceGroup &data_space_group, Logger &logger, size_t nthreads) {
const MRCandidate &c = result.candidates[result.chosen];
const auto t0 = std::chrono::steady_clock::now();
Logger quiet = Logger::Buffered();
ModelValidationResult v;
try {
v = ReadModelAnomalousMap(merged, *c.placed, data_space_group, c.k_sol, c.b_sol, quiet, nthreads);
} catch (const std::runtime_error &e) { // a GPU failure: no reading, so no relabelling
logger.Warning("Model prior: the anomalous map could not be read ({})", e.what());
}
auto &h = result.hand_check;
h.performed = true;
h.model_fits = c.null_sigma >= NULL_SIGMA_MIN;
h.fit_sigma = c.null_sigma;
h.adopted = v.ok && h.model_fits && !v.anomalous_hands_disagree;
h.anomalous_inverted = v.anomalous_hands_disagree;
h.anomalous_pairs = v.anomalous_pairs;
h.anomalous_scatterer_mean_sigma = v.anomalous_scatterer_mean_sigma;
h.deepest_site = v.anomalous_deepest_site;
h.deepest_sigma = v.anomalous_deepest_sigma;
h.seconds = seconds_since(t0);
if (!h.adopted) {
result.decision = "ENANTIOMORPH_VETOED";
logger.Info("Model prior: the enantiomorph {} is not taken into the answer - {} ({:.3f} s)",
c.space_group.xhm(),
h.anomalous_inverted ? fmt::format("the anomalous density at the placed model's atoms is "
"inverted ({} at {:.1f} sigma)", h.deepest_site, h.deepest_sigma)
: fmt::format("the placed model does not fit against its null ({:+.1f} sd)",
h.fit_sigma),
h.seconds);
}
}
Reference in New Issue View Git Blame Copy Permalink