leonarski_fandClaude Opus 5.5
af796d0a6a
Completes c6912e2a0 (whose message was a placeholder). What the two commits do:
- tNCS by half a lattice vector v: copies v apart weaken the reflections with h.v half-integral, and a
placement moved by (v + o)/2 (o an origin shift of the group) fits the strong ones exactly as well.
The candidate the choice falls on is compared with every such move of itself (rigid body refined,
packing) by R over the weakened reflections only; a move better by 0.01 is taken instead, and where
the two best are within 0.01 the merge is UNDECIDED.
- With tNCS, the group the merge was made in, searched to fewer copies than the chosen group, keeps
the merge UNDECIDED (not for its enantiomorph, which the hand check decides): its first copy alone
can sit at half the vector, and no placement of the chosen group, which may sit there itself, rules
it out.
- P1: the first copy sits at the origin, which no function of its position can single out (the
translation function is flat), so one copy stands out where its rotation-function Z reaches 8; the
null then judges it as before.
Checked against the deposited coordinates (placement comparison, symmetry/origin/setting matched):
two half-vector-shifted placements that had been decided (a SPACE_GROUP in the wrong setting, an
AGREES off by a/4) are UNDECIDED; two UNDECIDED tNCS sets are now AGREES with confirmed placements;
the nine confirmed tNCS placements and the ENANTIOMORPH stay; a P1 single copy with a perfect-match
prediction is placed (confirmed), a wrong PAE-domain placement in P1 (RF Z 6.4) is not.
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3
3076 lines
155 KiB
C++
3076 lines
155 KiB
C++
// 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 "ModelScaleGPU.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;
|
|
// A model whose expected log-likelihood gain (ExpectedLogLikelihoodGain) is below this is not searched for. The
|
|
// lowest of a model the search did place, over 193 survey sets: 26; every PAE domain below 24 placed nowhere.
|
|
constexpr double ELLG_MIN = 10.0;
|
|
// tNCS by half a lattice vector v: copies v apart make the reflections with h.v half-integral weak, and a
|
|
// placement moved by (v + o) / 2, o an origin shift of the group, fits the strong ones as well as the right one
|
|
// does; only the weak ones tell them apart. A placement is compared with every such move of itself on those
|
|
// reflections (R_tncs_weak), and one of them that fits them better by TNCS_SHIFT_MARGIN is taken instead; where the
|
|
// two best are closer than that, the placement is not decided on. A vector is taken as half-integral within
|
|
// TNCS_HALF_TOLERANCE.
|
|
constexpr double TNCS_SHIFT_MARGIN = 0.01;
|
|
constexpr double TNCS_HALF_TOLERANCE = 0.05;
|
|
|
|
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;
|
|
// Over the reflections a tNCS vector v weakens - h.v half-integral - work and free together (score_placement's
|
|
// `tncs`).
|
|
double r_tncs_weak = 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::Fractional *tncs = nullptr) {
|
|
const gemmi::GroupOps gops = sg.operations();
|
|
const gemmi::ReciprocalAsu asu(&sg);
|
|
struct Obs {
|
|
double f;
|
|
bool free;
|
|
bool weak; // h.v half-integral
|
|
};
|
|
std::unordered_map<long, Obs> 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;
|
|
bool weak = false;
|
|
if (tncs != nullptr) {
|
|
const double hv = h[0] * tncs->x + h[1] * tncs->y + h[2] * tncs->z;
|
|
weak = std::fabs(hv - std::round(hv)) > 0.25;
|
|
}
|
|
obs[hkl_key(h)] = {r.F, r.rfree_flag, weak};
|
|
}
|
|
double num_w = 0, den_w = 0, num_f = 0, den_f = 0, num_t = 0, den_t = 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.f, fc = std::abs(hv.value);
|
|
if (it->second.free) {
|
|
num_f += std::fabs(fo - fc);
|
|
den_f += fo;
|
|
} else {
|
|
num_w += std::fabs(fo - fc);
|
|
den_w += fo;
|
|
}
|
|
if (it->second.weak) {
|
|
num_t += std::fabs(fo - fc);
|
|
den_t += 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;
|
|
if (tncs != nullptr)
|
|
out.r_tncs_weak = den_t > 0 ? num_t / den_t : 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, each with the scale fit on the GPU beside it - or none,
|
|
// where there is no GPU, the cell is too small for its gridding or the two 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:
|
|
struct Engine {
|
|
std::unique_ptr<ModelStructureFactorsGPU> sf;
|
|
// The scale fit of the reflections the engine's F_calc are matched to, on its card: FitModelScale and
|
|
// gemmi's fit at a held solvent as ModelScaleGPU makes them - the same fits, equal to rounding.
|
|
std::mutex scale_m;
|
|
std::unique_ptr<CudaStream> stream;
|
|
std::unique_ptr<ModelScaleGPU> scale;
|
|
CudaDevicePtr<std::complex<float>> fcalc, fmask;
|
|
};
|
|
|
|
Engine *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);
|
|
const size_t points = sf->Reflections();
|
|
const size_t bytes = sf->DeviceBytes() + ModelScaleGPU::DeviceBytes(points)
|
|
+ 2 * points * sizeof(std::complex<float>);
|
|
if (sf->Supported() && bytes <= sf->DeviceTotalMemory() / 4) {
|
|
sf->Reserve();
|
|
const int device = RigidBodyGPUEngine::CurrentDevice();
|
|
set_gpu(sf->Device());
|
|
e.engine = std::make_unique<Engine>();
|
|
e.engine->stream = std::make_unique<CudaStream>();
|
|
e.engine->scale = std::make_unique<ModelScaleGPU>(*e.engine->stream, points);
|
|
e.engine->fcalc = CudaDevicePtr<std::complex<float>>(points);
|
|
e.engine->fmask = CudaDevicePtr<std::complex<float>>(points);
|
|
e.engine->sf = std::move(sf);
|
|
set_gpu(device);
|
|
}
|
|
} catch (const JFJochException &) {
|
|
e.engine.reset();
|
|
cuda_clear_error();
|
|
}
|
|
}
|
|
return e.engine.get();
|
|
}
|
|
|
|
private:
|
|
struct Entry {
|
|
bool tried = false;
|
|
std::unique_ptr<Engine> engine;
|
|
};
|
|
std::map<std::tuple<std::string, size_t, double>, Entry> engines_;
|
|
std::mutex m_;
|
|
};
|
|
|
|
// The scale of `scaling`'s points fitted on the engine's card: FitModelScale's grid over the solvent box, or
|
|
// with `solvent`, the overall scale and anisotropic B at its k_sol and b_sol held.
|
|
void fit_scale_gpu(ScoreEngines::Engine &e, gemmi::Scaling<float> &scaling, const PlacementScore *solvent,
|
|
double b_hi, size_t nthreads) {
|
|
const size_t n = scaling.points.size();
|
|
std::vector<std::array<int, 3>> hkl(n);
|
|
std::vector<double> stol2(n);
|
|
std::vector<float> fobs(n), sigma(n);
|
|
std::vector<std::complex<float>> fc(n), fm(n);
|
|
for (size_t i = 0; i < n; i++) {
|
|
const auto &p = scaling.points[i];
|
|
hkl[i] = {p.hkl[0], p.hkl[1], p.hkl[2]};
|
|
stol2[i] = p.stol2;
|
|
fobs[i] = p.fobs;
|
|
sigma[i] = p.sigma;
|
|
fc[i] = p.fcmol;
|
|
fm[i] = p.fmask;
|
|
}
|
|
std::vector<std::array<double, 6>> constraints;
|
|
for (const gemmi::Vec6 &c : scaling.constraint_matrix)
|
|
constraints.push_back(c);
|
|
double frac[9];
|
|
for (int i = 0; i < 3; i++)
|
|
for (int j = 0; j < 3; j++)
|
|
frac[3 * i + j] = scaling.cell.frac.mat[i][j];
|
|
|
|
std::unique_lock lock(e.scale_m);
|
|
const int device = RigidBodyGPUEngine::CurrentDevice();
|
|
set_gpu(e.sf->Device());
|
|
e.scale->SetPoints(hkl, stol2, fobs, sigma, constraints, frac);
|
|
if (cudaMemcpyAsync(e.fcalc.get(), fc.data(), n * sizeof(std::complex<float>), cudaMemcpyHostToDevice,
|
|
*e.stream) != cudaSuccess
|
|
|| cudaMemcpyAsync(e.fmask.get(), fm.data(), n * sizeof(std::complex<float>), cudaMemcpyHostToDevice,
|
|
*e.stream) != cudaSuccess)
|
|
throw JFJochException(JFJochExceptionCategory::GPUCUDAError, "Model prior: the scale fit's upload failed");
|
|
const auto *d_fc = reinterpret_cast<const float2 *>(e.fcalc.get());
|
|
const auto *d_fm = reinterpret_cast<const float2 *>(e.fmask.get());
|
|
ModelScaleParams p;
|
|
if (solvent != nullptr) {
|
|
scaling.k_sol = solvent->k_sol;
|
|
scaling.b_sol = solvent->b_sol;
|
|
p = e.scale->Fit(d_fc, d_fm, solvent->k_sol, solvent->b_sol);
|
|
} else {
|
|
try {
|
|
const ModelSolventFit fit = e.scale->FitSolvent(d_fc, d_fm, b_hi);
|
|
scaling.k_sol = fit.k_sol;
|
|
scaling.b_sol = fit.b_sol;
|
|
p = fit.scale;
|
|
} catch (const ModelScaleGPUTooFewReflections &) {
|
|
// FitModelScale walks such a grid in order, each point from where the last ended.
|
|
set_gpu(device);
|
|
lock.unlock();
|
|
FitModelScale(scaling, {.b_hi = b_hi}, nthreads);
|
|
return;
|
|
}
|
|
}
|
|
set_gpu(device);
|
|
scaling.k_overall = p.k_overall;
|
|
scaling.b_star = {p.b_star[0], p.b_star[1], p.b_star[2], p.b_star[3], p.b_star[4], p.b_star[5]};
|
|
}
|
|
#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, const gemmi::Fractional *tncs = nullptr) {
|
|
gemmi::AsuData<std::complex<float>> fmask, fmodel;
|
|
#ifdef JFJOCH_USE_CUDA
|
|
ScoreEngines::Engine *gpu = engines != nullptr ? engines->get(st, sg, d_min) : nullptr;
|
|
if (gpu)
|
|
gpu->sf->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);
|
|
#ifdef JFJOCH_USE_CUDA
|
|
if (gpu)
|
|
fit_scale_gpu(*gpu, scaling, solvent, SolventBMax(st.models[0]), nthreads);
|
|
else
|
|
#endif
|
|
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, tncs);
|
|
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;
|
|
}
|
|
|
|
// The tNCS vector where it is half a lattice vector - every component within TNCS_HALF_TOLERANCE of 0 or 1/2 -
|
|
// on those halves; empty otherwise.
|
|
std::optional<gemmi::Fractional> half_tncs_vector(const TranslationalNcs &tncs) {
|
|
if (!tncs.present)
|
|
return std::nullopt;
|
|
double v[3] = {tncs.vector.x, tncs.vector.y, tncs.vector.z};
|
|
for (double &x : v) {
|
|
const double h = std::round(2 * x) / 2;
|
|
if (std::fabs(x - h) > TNCS_HALF_TOLERANCE)
|
|
return std::nullopt;
|
|
x = h - std::floor(h);
|
|
}
|
|
return gemmi::Fractional(v[0], v[1], v[2]);
|
|
}
|
|
|
|
// Whether moving the origin by `s` (a multiple of 1/24 in each component) leaves `sg`'s operators as they are:
|
|
// a placement moved by s is then the same crystal.
|
|
bool is_origin_shift(const gemmi::SpaceGroup &sg, const gemmi::Fractional &s) {
|
|
gemmi::Op shift = gemmi::Op::identity();
|
|
const double x[3] = {s.x, s.y, s.z};
|
|
for (int a = 0; a < 3; a++)
|
|
shift.tran[a] = static_cast<int>(std::lround(x[a] * gemmi::Op::DEN)) % gemmi::Op::DEN;
|
|
gemmi::GroupOps ops = sg.operations();
|
|
ops.change_basis_forward(shift);
|
|
return ops.is_same_as(sg.operations());
|
|
}
|
|
|
|
} // 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,
|
|
size_t nthreads) {
|
|
const size_t terms = setup.Terms(), nops = setup.ops;
|
|
std::vector<gemmi::Mat33> rt;
|
|
std::vector<gemmi::Position> x;
|
|
for (size_t c = 0; c < rotations.size(); c++) {
|
|
rt.push_back(rotations[c].transpose());
|
|
x.push_back(cell.orthogonalize(positions[c]));
|
|
}
|
|
std::vector<std::complex<double>> out(terms, 0.0);
|
|
ParallelChunks(static_cast<int>(terms), ThreadsForWork(terms * nops * rotations.size(), nthreads),
|
|
[&](int lo, int hi) {
|
|
for (int t = lo; t < hi; t++)
|
|
for (size_t c = 0; c < rotations.size(); c++)
|
|
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[c].multiply(sj)) * setup.phase[t * nops + j]
|
|
* std::polar(1.0, 2.0 * PI * sj.dot(x[c]));
|
|
}
|
|
});
|
|
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();
|
|
ParallelChunks(static_cast<int>(g.size()), ThreadsForWork(g.size(), nthreads), [&](int lo, int hi) {
|
|
for (int e = lo; e < hi; 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 © : 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 expected log-likelihood gain of one copy of the model on the reflections of `merge` the search is made
|
|
// on: sum sigma_A^4 / 2, with sigma_A^2 the part of the scattering the copy is - one over the copy number the
|
|
// cell most probably holds, by Matthews - times Luzzati's D^2 for the coordinate error the pLDDTs predict.
|
|
// Following McCoy et al. (2017) Proc. Natl. Acad. Sci. USA 114, 3637-3641
|
|
double ExpectedLogLikelihoodGain(const SearchModel &model, const MRMerge &merge, double d_min, double d_max) {
|
|
const gemmi::UnitCell cell(merge.cell);
|
|
const SearchReflections refl = MakeSearchReflections(*merge.merged, cell, merge.space_group, d_min, d_max);
|
|
const double asu_volume = cell.volume / static_cast<double>(merge.space_group.operations().order());
|
|
const std::vector<double> pm = MatthewsProbabilities(asu_volume, model.mass_da);
|
|
size_t likely = 0;
|
|
for (size_t n = 0; n < pm.size(); n++)
|
|
if (pm[n] > pm[likely])
|
|
likely = n;
|
|
const double fraction = 1.0 / static_cast<double>(likely + 1);
|
|
double ellg = 0;
|
|
for (const gemmi::Miller &h : refl.hkl) {
|
|
const double d = std::exp(-2.0 * PI * PI * model.error_A * model.error_A * cell.calculate_1_d2(h) / 3.0);
|
|
ellg += fraction * fraction * d * d * d * d / 2.0;
|
|
}
|
|
return ellg;
|
|
}
|
|
|
|
// 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);
|
|
// What a placement of the model can gain at best: a model too small a part of the cell, or too far from
|
|
// it, for its placement to stand out is not searched for.
|
|
const double ellg = ExpectedLogLikelihoodGain(model, merges[0], result.d_min_search, result.d_max_search);
|
|
logger.Info("Model prior: {}: expected LLG {:.1f}", model.description, ellg);
|
|
if (ellg < ELLG_MIN) {
|
|
result.ok = true;
|
|
result.decision = fmt::format("NOT_PLACED: an expected LLG of {:.1f}, below {:.0f}", ellg, ELLG_MIN);
|
|
return result;
|
|
}
|
|
|
|
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);
|
|
// Each group's on a thread of its own.
|
|
std::vector<std::future<std::vector<TranslationPeak>>> tf_running;
|
|
for (const gemmi::SpaceGroup &sg : groups)
|
|
tf_running.push_back(std::async(std::launch::async, [&, sg] {
|
|
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]);
|
|
});
|
|
return tf;
|
|
}));
|
|
std::vector<TranslationPeak> tf_all;
|
|
for (auto &r : tf_running) {
|
|
const std::vector<TranslationPeak> tf = r.get();
|
|
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();
|
|
// The further copies' peaks are ranked by the likelihood on the search's own reflections, read off
|
|
// the molecular transform. Each made where it is first needed; where the likelihood's reflections are
|
|
// the search's own, one serves both.
|
|
std::optional<LikelihoodSetup> lsetup_made, rsetup_made;
|
|
const auto lsetup = [&]() -> const LikelihoodSetup & {
|
|
if (!lsetup_made)
|
|
lsetup_made.emplace(MakeLikelihoodSetup(lrefl, cell, sg));
|
|
return *lsetup_made;
|
|
};
|
|
const auto rsetup = [&]() -> const LikelihoodSetup & {
|
|
if (d_llg == result.d_min_search)
|
|
return lsetup();
|
|
if (!rsetup_made)
|
|
rsetup_made.emplace(MakeLikelihoodSetup(refl, cell, sg));
|
|
return *rsetup_made;
|
|
};
|
|
|
|
// --- 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.
|
|
// The placements that go through it are known before any does, so they run through it side by side,
|
|
// on an engine each, and are then taken in order.
|
|
struct First {
|
|
size_t i = 0;
|
|
std::shared_ptr<gemmi::Structure> placed;
|
|
double clash = 1;
|
|
bool refine = false;
|
|
PlacementScore score{NAN, NAN, NAN, NAN};
|
|
gemmi::Mat33 moved;
|
|
};
|
|
std::vector<First> first;
|
|
int refined = 0;
|
|
for (size_t k = 0; k < by_tf.size() && refined < (fits ? RESCORED_PER_GROUP : 1); k++) {
|
|
First f;
|
|
f.i = by_tf[k];
|
|
f.placed = std::make_shared<gemmi::Structure>(PlaceSearchModel(model, rots[f.i], tf[f.i].t, cell, sg));
|
|
f.clash = ClashFraction(*f.placed);
|
|
f.refine = fits && f.clash <= 2 * CLASH_MAX;
|
|
if (f.refine || !fits)
|
|
refined++;
|
|
first.push_back(std::move(f));
|
|
}
|
|
// The further copies' translation function, made while the first copy is refined.
|
|
std::future<TranslationSetup> setup2_made;
|
|
if (fits && cand.copies_max > 1)
|
|
setup2_made = std::async(std::launch::async, [&] { return MakeTranslationSetup(refl, cell, sg, true); });
|
|
std::vector<First *> to_refine;
|
|
for (First &f : first)
|
|
if (f.refine)
|
|
to_refine.push_back(&f);
|
|
if (!to_refine.empty()) {
|
|
engines = make_engines(*to_refine[0]->placed, sg, *merge.merged, d_min, to_refine.size(), logger);
|
|
std::vector<std::future<void>> running;
|
|
for (First *f : to_refine)
|
|
running.push_back(std::async(std::launch::async, [&, f] {
|
|
f->score = refine_placement(*f->placed, sg, *merge.merged, d_min, engines.get(), nthreads,
|
|
&scorers, &f->moved);
|
|
f->clash = ClashFraction(*f->placed);
|
|
}));
|
|
for (auto &r : running)
|
|
r.get();
|
|
}
|
|
for (const First &f : first) {
|
|
if (!f.refine && fits && have)
|
|
continue;
|
|
const size_t i = f.i;
|
|
const PlacementScore &score = f.score;
|
|
// 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 = f.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 = f.placed;
|
|
cand.copies = 1;
|
|
orientation = rots[i];
|
|
copies = {*f.placed};
|
|
copy_rotation = {f.moved.multiply(rots[i])};
|
|
copy_position = {centroid(*f.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 * f.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 std::vector<std::complex<double>> no_known(lsetup().Terms(), 0.0);
|
|
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 = setup2_made.get();
|
|
const std::vector<std::complex<double>> known = KnownStructureFactors(ft, *setup2, cell, copy_rotation,
|
|
copy_position, nthreads);
|
|
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;
|
|
// The CopyFactors the walk down the peaks may need, all made at once.
|
|
const auto factors_of = [&](size_t i) -> std::vector<std::complex<double>> & {
|
|
return i < rots.size() ? rot_g[i] : copy_rg[i - rots.size()];
|
|
};
|
|
std::vector<size_t> missing;
|
|
for (size_t k = 0; k < by_tf2.size() && k < 2 * COPY_TRIALS; k++)
|
|
if (factors_of(by_tf2[k]).empty())
|
|
missing.push_back(by_tf2[k]);
|
|
rsetup();
|
|
ParallelFor(static_cast<int>(missing.size()), nthreads, [&](int m) {
|
|
factors_of(missing[m]) = CopyFactors(model, rsetup(), cell, rots2[missing[m]],
|
|
result.d_min_search, nthreads, &ft);
|
|
});
|
|
// Their gains, all taken at once.
|
|
std::vector<double> gains(std::min(by_tf2.size(), 2 * COPY_TRIALS));
|
|
ParallelFor(static_cast<int>(gains.size()), nthreads, [&](int k) {
|
|
std::vector<std::complex<double>> f = rknown;
|
|
AddCopy(rsetup(), factors_of(by_tf2[k]), tf2[by_tf2[k]].t, f);
|
|
gains[k] = LogLikelihoodGain(rsetup(), f, model.error_A, known_scale) - rknown_llg;
|
|
});
|
|
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];
|
|
trial.gain = gains[k];
|
|
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. In P1 the first copy sits
|
|
// at the origin, which is arbitrary there and which no function of its position can single out: one copy
|
|
// stands out by its orientation alone, where the rotation function's Z-score reaches the same 8.
|
|
else if (!(cand.tf_z >= TF_Z_PLACED || (cand.llg_z >= TF_Z_PLACED && cand.copies >= cand.copies_likely)
|
|
|| (gops.order() == 1 && cand.copies == 1 && cand.rf_z >= TF_Z_PLACED)))
|
|
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 {
|
|
// The null of candidate i.
|
|
const auto run_null = [&](int i) {
|
|
const MRCandidate &c = result.candidates[i];
|
|
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;
|
|
};
|
|
return null_sigma(*c.placed, c.copies, orientations[i], c.space_group, *merges[merge].merged,
|
|
result.d_min_score, engines.get(), nthreads, &scorers, translate);
|
|
};
|
|
// The first candidate's null alone; where it fails, the nulls of all the other placed candidates side
|
|
// by side, each taken where the choice comes to it.
|
|
std::vector<double> nulls(result.candidates.size(), NAN);
|
|
bool first = true;
|
|
// With tNCS by half a lattice vector, candidate i against its moves by half that vector (TNCS_SHIFT_MARGIN):
|
|
// the better taken into it, and whether the two best are too close to tell apart.
|
|
// The tNCS vector is the answer's merge's, in its cell.
|
|
const std::optional<gemmi::Fractional> half = merge == 0 ? half_tncs_vector(result.tncs) : std::nullopt;
|
|
std::vector<char> shift_checked(result.candidates.size(), 0);
|
|
const auto check_shifts = [&](int i) -> bool {
|
|
MRCandidate &c = result.candidates[i];
|
|
const gemmi::Fractional v = *half;
|
|
const std::vector<MergedReflection> &merged = *merges[merge].merged;
|
|
const double d_min = result.d_min_score;
|
|
struct Option {
|
|
gemmi::Fractional shift;
|
|
std::shared_ptr<gemmi::Structure> placed;
|
|
PlacementScore score;
|
|
double clash = 0;
|
|
};
|
|
std::vector<Option> options{{gemmi::Fractional(0, 0, 0), c.placed,
|
|
score_placement(*c.placed, c.space_group, merged, d_min, nthreads, &scorers,
|
|
nullptr, nullptr, &v), c.clash_fraction}};
|
|
RigidBodyEngines engines;
|
|
for (int o = 0; o < 8; o++) {
|
|
const gemmi::Fractional shift((v.x + 0.5 * (o & 1)) / 2, (v.y + 0.5 * ((o >> 1) & 1)) / 2,
|
|
(v.z + 0.5 * ((o >> 2) & 1)) / 2);
|
|
// The same crystal as one already tried, or as no move at all.
|
|
if (std::any_of(options.begin(), options.end(), [&](const Option &p) {
|
|
return is_origin_shift(c.space_group, gemmi::Fractional(shift.x - p.shift.x, shift.y - p.shift.y,
|
|
shift.z - p.shift.z));
|
|
}))
|
|
continue;
|
|
auto moved = std::make_shared<gemmi::Structure>(*c.placed);
|
|
const gemmi::Position by = moved->cell.orthogonalize(shift);
|
|
std::vector<gemmi::Position> pos = ModelPositions(moved->models[0]);
|
|
for (gemmi::Position &p : pos)
|
|
p = gemmi::Position(p + by);
|
|
SetModelPositions(moved->models[0], pos);
|
|
Option opt{shift, moved, {}, ClashFraction(*moved, c.copies)};
|
|
if (opt.clash <= 2 * CLASH_MAX) {
|
|
if (engines.get() == nullptr)
|
|
engines = make_engines(*moved, c.space_group, merged, d_min, 1, logger);
|
|
refine_rigid_body(*moved, c.space_group, merged, d_min, engines.get(), nthreads);
|
|
opt.clash = ClashFraction(*moved, c.copies);
|
|
opt.score = score_placement(*moved, c.space_group, merged, d_min, nthreads, &scorers, nullptr,
|
|
nullptr, &v);
|
|
}
|
|
logger.Info("Model prior: {} moved by {:.2f} {:.2f} {:.2f}: clash {:.1f}%, R-free {:.3f}, R on the "
|
|
"reflections tNCS weakens {:.3f} against {:.3f}", c.space_group.xhm(), shift.x, shift.y,
|
|
shift.z, 100 * opt.clash, opt.score.r_free, opt.score.r_tncs_weak,
|
|
options[0].score.r_tncs_weak);
|
|
options.push_back(std::move(opt));
|
|
}
|
|
std::vector<const Option *> packing;
|
|
for (const Option &p : options)
|
|
if (p.clash <= CLASH_MAX && std::isfinite(p.score.r_tncs_weak))
|
|
packing.push_back(&p);
|
|
std::stable_sort(packing.begin(), packing.end(), [](const Option *a, const Option *b) {
|
|
return a->score.r_tncs_weak < b->score.r_tncs_weak;
|
|
});
|
|
if (packing.empty() || packing[0] == &options[0])
|
|
return packing.size() > 1 && packing[1]->score.r_tncs_weak - packing[0]->score.r_tncs_weak < TNCS_SHIFT_MARGIN;
|
|
const Option &b = *packing[0];
|
|
if (options[0].score.r_tncs_weak - b.score.r_tncs_weak < TNCS_SHIFT_MARGIN
|
|
|| (packing.size() > 1 && packing[1]->score.r_tncs_weak - b.score.r_tncs_weak < TNCS_SHIFT_MARGIN))
|
|
return true;
|
|
logger.Info("Model prior: {} taken moved by {:.2f} {:.2f} {:.2f}", c.space_group.xhm(), b.shift.x,
|
|
b.shift.y, b.shift.z);
|
|
c.placed = b.placed;
|
|
c.clash_fraction = b.clash;
|
|
c.r_work = b.score.r_work;
|
|
c.r_free = b.score.r_free;
|
|
c.k_sol = b.score.k_sol;
|
|
c.b_sol = b.score.b_sol;
|
|
return false;
|
|
};
|
|
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];
|
|
if (half && !shift_checked[best]) {
|
|
shift_checked[best] = 1;
|
|
if (check_shifts(best)) {
|
|
result.decision = fmt::format("UNDECIDED: {}'s placement and its move by half the tNCS vector fit "
|
|
"alike", c.space_group.xhm());
|
|
return -2;
|
|
}
|
|
continue; // its R-free may have changed
|
|
}
|
|
if (std::isnan(nulls[best])) {
|
|
const auto t_null = std::chrono::steady_clock::now();
|
|
std::vector<int> todo{best};
|
|
if (!first)
|
|
for (size_t i = 0; i < result.candidates.size(); i++)
|
|
if (static_cast<int>(i) != best && result.candidates[i].merge == merge
|
|
&& result.candidates[i].verdict == "PLACED")
|
|
todo.push_back(static_cast<int>(i));
|
|
std::vector<std::future<double>> running;
|
|
for (int i : todo)
|
|
running.push_back(std::async(std::launch::async, run_null, i));
|
|
for (size_t k = 0; k < todo.size(); k++)
|
|
nulls[todo[k]] = running[k].get();
|
|
result.t_null += seconds_since(t_null);
|
|
first = false;
|
|
}
|
|
c.null_sigma = nulls[best];
|
|
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];
|
|
// With tNCS, the group the merge was made in, searched to fewer copies than the chosen one - whether it
|
|
// packs or not: its first copy alone can sit where half the tNCS vector puts it - is not ruled out by
|
|
// any placement of the chosen one, which may sit there itself. (Its enantiomorph is the hand check's.)
|
|
if (result.tncs.present && static_cast<int>(i) != best && c.merge == merge
|
|
&& c.space_group.xhm() == merges[merge].space_group.xhm()
|
|
&& result.candidates[best].space_group.xhm() != c.space_group.xhm()
|
|
&& !enantiomorphs(c.space_group, result.candidates[best].space_group)
|
|
&& c.copies < result.candidates[best].copies) {
|
|
result.decision = fmt::format("UNDECIDED: {} placed fewer copies than {} beside tNCS",
|
|
c.space_group.xhm(), result.candidates[best].space_group.xhm());
|
|
return -2;
|
|
}
|
|
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);
|
|
}
|
|
}
|