leonarski_fandClaude Opus 5.5
6c8ead466a
The search was CPU-bound: on a 54-chain partial complex ~80% of all CPU samples were the score's FitModelScale grid (and libm's exp under it), with the GPU busy 13 of 37 s. - score_placement fits its scale on the GPU where its structure factors are made there (ModelScaleGPU's FitSolvent / Fit, as the rigid body already does); equal to rounding, the CPU build unchanged. - A model whose expected log-likelihood gain (sum sigma_A^4 / 2, McCoy et al. 2017; sigma_A^2 the part of the scattering one copy is by Matthews, times Luzzati's D^2) is below 10 is not searched for. Over the 193-set survey the lowest eLLG of a placed model is 26; every PAE domain below 24 placed nowhere, so no decision changes. - The likelihood setups made where first needed, one serving both where the reflections are the same. - After a candidate fails its null, the other placed candidates' nulls side by side. - The first copy's placements refined side by side on an engine each; the groups' translation functions side by side; CopyFactors and KnownStructureFactors over threads; the further copies' CopyFactors and gains for the peaks walked all at once; the known-structure TF setup made while the first copy is refined. All exact. Survey replay (193 sets): no decision changes; two R values differ in the 3rd decimal (GPU scale). Alternating runs, MODEL_PRIOR_SECONDS base -> new: 54-chain complex 38.3 -> 14.7, 4mtk 18.3 -> 9.7, 9v16 13.7 -> 8.3 (shared machine, indicative). Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3
2913 lines
145 KiB
C++
2913 lines
145 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;
|
|
|
|
double seconds_since(std::chrono::steady_clock::time_point t0) {
|
|
return std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count();
|
|
}
|
|
|
|
// A rotation drawn uniformly from SO(3), through a uniform random unit quaternion.
|
|
// Following Shoemake (1992) Graphics Gems III, 124-132
|
|
gemmi::Mat33 random_rotation(std::mt19937 &rng) {
|
|
std::uniform_real_distribution<double> u(0.0, 1.0);
|
|
const double u1 = u(rng), t2 = 2 * PI * u(rng), t3 = 2 * PI * u(rng);
|
|
const double r1 = std::sqrt(1 - u1), r2 = std::sqrt(u1);
|
|
const double x = r1 * std::sin(t2), y = r1 * std::cos(t2), z = r2 * std::sin(t3), w = r2 * std::cos(t3);
|
|
return {1 - 2 * (y * y + z * z), 2 * (x * y - z * w), 2 * (x * z + y * w),
|
|
2 * (x * y + z * w), 1 - 2 * (x * x + z * z), 2 * (y * z - x * w),
|
|
2 * (x * z - y * w), 2 * (y * z + x * w), 1 - 2 * (x * x + y * y)};
|
|
}
|
|
|
|
// A rotation by `deg` about Cartesian axis 0, 1 or 2.
|
|
gemmi::Mat33 axis_rotation(int axis, double deg) {
|
|
const double c = std::cos(deg * PI / 180.0), s = std::sin(deg * PI / 180.0);
|
|
if (axis == 0)
|
|
return {1, 0, 0, 0, c, -s, 0, s, c};
|
|
if (axis == 1)
|
|
return {c, 0, s, 0, 1, 0, -s, 0, c};
|
|
return {c, -s, 0, s, c, 0, 0, 0, 1};
|
|
}
|
|
|
|
double trace(const gemmi::Mat33 &m) {
|
|
return m.a[0][0] + m.a[1][1] + m.a[2][2];
|
|
}
|
|
|
|
gemmi::Vec3 miller_vec(const std::array<int, 3> &h) {
|
|
return {static_cast<double>(h[0]), static_cast<double>(h[1]), static_cast<double>(h[2])};
|
|
}
|
|
|
|
// h R for a symmetry operator's rotation: the index the operator takes h to.
|
|
gemmi::Miller rotated_index(const gemmi::Op &op, const gemmi::Miller &h) {
|
|
const auto k = op.apply_to_hkl_without_division(h);
|
|
return {{k[0] / gemmi::Op::DEN, k[1] / gemmi::Op::DEN, k[2] / gemmi::Op::DEN}};
|
|
}
|
|
|
|
long hkl_key(const gemmi::Miller &h) {
|
|
return (static_cast<long>(h[0]) + 512) * 1048576L + (static_cast<long>(h[1]) + 512) * 1024L + (h[2] + 512);
|
|
}
|
|
|
|
int good_even_size(int n) {
|
|
while (n % 2 != 0 || !gemmi::has_small_factorization(n))
|
|
++n;
|
|
return n;
|
|
}
|
|
|
|
// --- placements scored against the data: R-work / R-free after the bulk solvent and the scale ---
|
|
|
|
struct PlacementScore {
|
|
double r_work = 1, r_free = 1;
|
|
double k_sol = NAN, b_sol = NAN;
|
|
// Against the data as twinned by a law (score_placement's `law`): at the fraction that gives the lowest
|
|
// R-work.
|
|
double r_work_twin = NAN, r_free_twin = NAN, twin_fraction = NAN;
|
|
};
|
|
|
|
gemmi::AsuData<gemmi::ValueSigma<float>> working_amplitudes(const std::vector<MergedReflection> &merged,
|
|
const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
|
|
double d_min) {
|
|
const gemmi::GroupOps gops = sg.operations();
|
|
const gemmi::ReciprocalAsu asu(&sg);
|
|
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
|
|
fobs.unit_cell_ = cell;
|
|
fobs.spacegroup_ = &sg;
|
|
for (const MergedReflection &r : merged) {
|
|
if (std::isnan(r.F) || r.rfree_flag || r.d < d_min)
|
|
continue;
|
|
gemmi::Miller h{{r.h, r.k, r.l}};
|
|
if (!asu.is_in(h))
|
|
h = asu.to_asu(h, gops).first;
|
|
fobs.v.push_back({h, {r.F, 1.0f}});
|
|
}
|
|
fobs.ensure_asu();
|
|
fobs.ensure_sorted();
|
|
return fobs;
|
|
}
|
|
|
|
// The model's Fcalc and bulk-solvent mask to d_min on the CPU, as model validation makes them.
|
|
void model_factors(const gemmi::Structure &st, double d_min, gemmi::AsuData<std::complex<float>> &fmodel,
|
|
gemmi::AsuData<std::complex<float>> &fmask, size_t nthreads) {
|
|
gemmi::DensityCalculator<Table, float> dc;
|
|
dc.d_min = d_min;
|
|
dc.rate = 1.5;
|
|
dc.set_grid_cell_and_spacegroup(st);
|
|
dc.set_refmac_compatible_blur(st.models[0]);
|
|
gemmi::Grid<float> mask;
|
|
mask.unit_cell = dc.grid.unit_cell;
|
|
mask.spacegroup = dc.grid.spacegroup;
|
|
mask.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
|
|
const std::vector<size_t> leaders = OrbitLeaders(mask, nthreads);
|
|
PutMaskOnGrid(mask, st.models[0], leaders, nthreads);
|
|
fmask = MapToFPhi(mask).prepare_asu_data(d_min, 0);
|
|
PutModelDensityOnGrid(dc, st.models[0], leaders, nthreads);
|
|
fmodel = MapToFPhi(dc.grid).prepare_asu_data(d_min, dc.blur, false, false, false);
|
|
}
|
|
|
|
// R-work and R-free of the scaled Fcalc over the reflections to d_min.
|
|
PlacementScore r_factors(const gemmi::AsuData<std::complex<float>> &fmodel, const gemmi::Scaling<float> &scaling,
|
|
const gemmi::SpaceGroup &sg, const std::vector<MergedReflection> &merged, double d_min) {
|
|
const gemmi::GroupOps gops = sg.operations();
|
|
const gemmi::ReciprocalAsu asu(&sg);
|
|
std::unordered_map<long, std::pair<double, bool>> obs;
|
|
for (const MergedReflection &r : merged) {
|
|
if (std::isnan(r.F) || r.d < d_min)
|
|
continue;
|
|
gemmi::Miller h{{r.h, r.k, r.l}};
|
|
if (!asu.is_in(h))
|
|
h = asu.to_asu(h, gops).first;
|
|
obs[hkl_key(h)] = {r.F, r.rfree_flag};
|
|
}
|
|
double num_w = 0, den_w = 0, num_f = 0, den_f = 0;
|
|
for (const auto &hv : fmodel.v) {
|
|
const auto it = obs.find(hkl_key(hv.hkl));
|
|
if (it == obs.end())
|
|
continue;
|
|
const double fo = it->second.first, fc = std::abs(hv.value);
|
|
if (it->second.second) {
|
|
num_f += std::fabs(fo - fc);
|
|
den_f += fo;
|
|
} else {
|
|
num_w += std::fabs(fo - fc);
|
|
den_w += fo;
|
|
}
|
|
}
|
|
PlacementScore out;
|
|
out.k_sol = scaling.k_sol;
|
|
out.b_sol = scaling.b_sol;
|
|
out.r_work = den_w > 0 ? num_w / den_w : 1;
|
|
out.r_free = den_f > 0 ? num_f / den_f : 1;
|
|
return out;
|
|
}
|
|
|
|
// Model validation's structure-factor engines on the GPU, one per group, copy count and resolution, made
|
|
// where the first placement of that kind is scored, 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) {
|
|
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);
|
|
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.
|
|
// `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, 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::Position> pos = ModelPositions(rep.models[0]);
|
|
gemmi::Vec3 c;
|
|
for (const gemmi::Position &p : pos)
|
|
c += p;
|
|
c *= 1.0 / static_cast<double>(pos.size());
|
|
const gemmi::Vec3 to = translate ? gemmi::Vec3(placed.cell.orthogonalize(where[i])) : c;
|
|
for (gemmi::Position &p : pos)
|
|
p = gemmi::Position(draws[i].multiply(gemmi::Vec3(p) - c) + to);
|
|
SetModelPositions(rep.models[0], pos);
|
|
r_work[i] = refine_placement(rep, sg, coarse, null_d_min, pool, nthreads, scorers).r_work;
|
|
}
|
|
}));
|
|
for (auto &w : workers)
|
|
w.get();
|
|
const double mean = std::accumulate(r_work.begin(), r_work.end(), 0.0) / NULL_REPLICATES;
|
|
double s2 = 0;
|
|
for (double x : r_work)
|
|
s2 += (x - mean) * (x - mean);
|
|
const double sd = std::sqrt(s2 / (NULL_REPLICATES - 1));
|
|
return sd > 0 ? (mean - real) / sd : 0.0;
|
|
}
|
|
|
|
bool same_rotations_and_centring(const gemmi::GroupOps &a, const gemmi::GroupOps &b) {
|
|
if (a.sym_ops.size() != b.sym_ops.size() || a.cen_ops.size() != b.cen_ops.size())
|
|
return false;
|
|
for (const gemmi::Op &op : a.sym_ops)
|
|
if (b.find_by_rotation(op.rot) == nullptr)
|
|
return false;
|
|
for (const gemmi::Op::Tran &c : a.cen_ops)
|
|
if (std::find(b.cen_ops.begin(), b.cen_ops.end(), c) == b.cen_ops.end())
|
|
return false;
|
|
return true;
|
|
}
|
|
|
|
// Whether b is a on a shifted origin, the shift on the twelfths every crystallographic one lies on.
|
|
bool same_up_to_origin(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b) {
|
|
if (a.number != b.number)
|
|
return false;
|
|
const gemmi::GroupOps bops = b.operations();
|
|
for (int i = 0; i < 12; i++)
|
|
for (int j = 0; j < 12; j++)
|
|
for (int k = 0; k < 12; k++) {
|
|
gemmi::Op shift = gemmi::Op::identity();
|
|
shift.tran = {i * gemmi::Op::DEN / 12, j * gemmi::Op::DEN / 12, k * gemmi::Op::DEN / 12};
|
|
gemmi::GroupOps aops = a.operations();
|
|
aops.change_basis_forward(shift);
|
|
if (aops.is_same_as(bops))
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
bool enantiomorphs(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b) {
|
|
if (!a.is_enantiomorphic())
|
|
return false;
|
|
gemmi::GroupOps e = a.operations();
|
|
e.change_basis_forward(a.change_of_hand_op());
|
|
const gemmi::SpaceGroup *other = gemmi::find_spacegroup_by_ops(e);
|
|
return other != nullptr && other->number == b.number;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
// --- the model ---
|
|
|
|
namespace {
|
|
|
|
// The atoms' centroid moved to the origin, and the radius, mass and coordinate error of what is left.
|
|
void finish_search_model(SearchModel &m) {
|
|
std::vector<gemmi::Position> pos = ModelPositions(m.st.models[0]);
|
|
double error2 = 0;
|
|
for (const gemmi::Chain &ch : m.st.models[0].chains)
|
|
for (const gemmi::Residue &r : ch.residues)
|
|
for (const gemmi::Atom &a : r.atoms)
|
|
error2 += 3.0 * a.b_iso / (8.0 * PI * PI); // the B from pLDDT, back to its rmsd^2
|
|
m.error_A = std::sqrt(error2 / static_cast<double>(pos.size()));
|
|
gemmi::Vec3 c;
|
|
for (const gemmi::Position &p : pos)
|
|
c += p;
|
|
c *= 1.0 / static_cast<double>(pos.size());
|
|
m.radius_A = 0;
|
|
for (gemmi::Position &p : pos) {
|
|
p = gemmi::Position(gemmi::Vec3(p) - c);
|
|
m.radius_A = std::max(m.radius_A, p.length());
|
|
}
|
|
SetModelPositions(m.st.models[0], pos);
|
|
m.mass_da = DALTON_PER_RESIDUE * m.residues_kept;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
SearchModel PrepareSearchModel(const std::string &path) {
|
|
SearchModel out;
|
|
try {
|
|
out.st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
|
|
} catch (const std::exception &e) {
|
|
out.failure_reason = std::string("cannot read the model: ") + e.what();
|
|
return out;
|
|
}
|
|
if (out.st.models.empty()) {
|
|
out.failure_reason = "the model has no atoms";
|
|
return out;
|
|
}
|
|
out.st.models.resize(1);
|
|
gemmi::remove_hydrogens(out.st);
|
|
gemmi::Model &model = out.st.models[0];
|
|
for (gemmi::Chain &ch : model.chains) {
|
|
std::vector<gemmi::Residue> kept;
|
|
for (gemmi::Residue &r : ch.residues) {
|
|
for (const gemmi::Atom &a : r.atoms)
|
|
if (a.b_iso > 100.0f) {
|
|
out.failure_reason = "the B-factor column holds values above 100, so it is not a pLDDT: "
|
|
"not a predicted model";
|
|
return out;
|
|
}
|
|
// The residues a prediction is confident of: pLDDT 70 and above, as Oeffner et al. (2022)
|
|
// Acta Cryst. D78, 1303-1314 trim an AlphaFold model for molecular replacement. Each kept with
|
|
// its position in the whole prediction, which its predicted aligned error is indexed by.
|
|
if (!r.atoms.empty() && r.atoms[0].b_iso >= 70.0f) {
|
|
kept.push_back(std::move(r));
|
|
out.residue_index.push_back(out.residues_read);
|
|
}
|
|
out.residues_read++;
|
|
}
|
|
ch.residues = std::move(kept);
|
|
out.residues_kept += static_cast<int>(ch.residues.size());
|
|
}
|
|
std::erase_if(model.chains, [](const gemmi::Chain &ch) { return ch.residues.empty(); });
|
|
if (out.residues_kept == 0) {
|
|
out.failure_reason = "no residue of the model has a pLDDT of 70 or more";
|
|
return out;
|
|
}
|
|
// B from pLDDT through the coordinate error it predicts: rmsd = 1.5 exp(4 (0.7 - pLDDT/100)) A
|
|
// and B = 8 pi^2 rmsd^2 / 3. Following Oeffner et al. (2022) Acta Cryst. D78, 1303-1314
|
|
for (gemmi::Chain &ch : model.chains)
|
|
for (gemmi::Residue &r : ch.residues)
|
|
for (gemmi::Atom &a : r.atoms) {
|
|
const double rmsd = 1.5 * std::exp(4.0 * (0.7 - a.b_iso / 100.0));
|
|
a.b_iso = static_cast<float>(8.0 * PI * PI * rmsd * rmsd / 3.0);
|
|
}
|
|
finish_search_model(out);
|
|
return out;
|
|
}
|
|
|
|
// --- domains from the predicted aligned error ---
|
|
|
|
std::string PaePathFor(const std::string &model_path) {
|
|
const size_t slash = model_path.find_last_of("/\\");
|
|
const std::string dir = slash == std::string::npos ? std::string() : model_path.substr(0, slash + 1);
|
|
const std::string name = model_path.substr(dir.size());
|
|
const size_t at = name.rfind("-model_v");
|
|
if (at == std::string::npos)
|
|
return {};
|
|
const size_t dot = name.find('.', at);
|
|
const std::string version = name.substr(at + 8, dot == std::string::npos ? std::string::npos : dot - at - 8);
|
|
return dir + name.substr(0, at) + "-predicted_aligned_error_v" + version + ".json";
|
|
}
|
|
|
|
// A prediction cut into domains by its predicted aligned error before molecular replacement, as Oeffner et
|
|
// al. (2022) Acta Cryst. D78, 1303-1314 do; here by average linkage.
|
|
std::vector<std::vector<int>> PaeDomains(const std::string &pae_path, const SearchModel &model) {
|
|
// Two domains a prediction does not place relative to each other with confidence are this far apart in
|
|
// mean PAE (within one domain it is a few angstroms); a domain is at least this many residues.
|
|
constexpr double PAE_DOMAIN_A = 10.0;
|
|
constexpr int DOMAIN_MIN_RESIDUES = 40;
|
|
// The clustering starts from runs of this many consecutive residues kept.
|
|
constexpr int PIECE = 5;
|
|
std::vector<std::vector<double>> pae;
|
|
try {
|
|
std::ifstream in(pae_path);
|
|
if (!in)
|
|
return {};
|
|
const nlohmann::json j = nlohmann::json::parse(in);
|
|
const nlohmann::json &m = j.is_array() ? j.at(0).at("predicted_aligned_error") : j.at("predicted_aligned_error");
|
|
for (const auto &row : m)
|
|
pae.push_back(row.get<std::vector<double>>());
|
|
} catch (const std::exception &) {
|
|
return {};
|
|
}
|
|
const size_t n = pae.size();
|
|
if (n != static_cast<size_t>(model.residues_read)
|
|
|| std::any_of(pae.begin(), pae.end(), [n](const std::vector<double> &r) { return r.size() != n; }))
|
|
return {};
|
|
// The pieces: runs of consecutive kept residues, at most PIECE long.
|
|
std::vector<std::vector<int>> members;
|
|
for (size_t i = 0; i < model.residue_index.size(); i++) {
|
|
const int r = model.residue_index[i];
|
|
if (members.empty() || static_cast<int>(members.back().size()) == PIECE || members.back().back() != r - 1)
|
|
members.push_back({});
|
|
members.back().push_back(r);
|
|
}
|
|
// Average linkage on the symmetrised PAE: the two clusters with the lowest mean PAE between them merged
|
|
// until that is PAE_DOMAIN_A or more.
|
|
const size_t m = members.size();
|
|
std::vector<double> sum(m * m, 0.0);
|
|
for (size_t a = 0; a < m; a++)
|
|
for (size_t b = 0; b < m; b++)
|
|
for (int i : members[a])
|
|
for (int k : members[b])
|
|
sum[a * m + b] += (pae[i][k] + pae[k][i]) / 2.0;
|
|
std::vector<char> active(m, 1);
|
|
while (true) {
|
|
double best = PAE_DOMAIN_A;
|
|
size_t ba = 0, bb = 0;
|
|
for (size_t a = 0; a < m; a++)
|
|
for (size_t b = a + 1; b < m; b++)
|
|
if (active[a] && active[b]) {
|
|
const double mean = sum[a * m + b] / static_cast<double>(members[a].size() * members[b].size());
|
|
if (mean < best) {
|
|
best = mean;
|
|
ba = a;
|
|
bb = b;
|
|
}
|
|
}
|
|
if (!(best < PAE_DOMAIN_A))
|
|
break;
|
|
for (size_t c = 0; c < m; c++) {
|
|
sum[ba * m + c] += sum[bb * m + c];
|
|
sum[c * m + ba] += sum[c * m + bb];
|
|
}
|
|
members[ba].insert(members[ba].end(), members[bb].begin(), members[bb].end());
|
|
active[bb] = 0;
|
|
}
|
|
std::vector<std::vector<int>> out;
|
|
for (size_t a = 0; a < m; a++)
|
|
if (active[a] && static_cast<int>(members[a].size()) >= DOMAIN_MIN_RESIDUES) {
|
|
std::sort(members[a].begin(), members[a].end());
|
|
out.push_back(members[a]);
|
|
}
|
|
std::stable_sort(out.begin(), out.end(), [](const std::vector<int> &a, const std::vector<int> &b) {
|
|
return a.size() > b.size();
|
|
});
|
|
return out;
|
|
}
|
|
|
|
SearchModel DomainModel(const SearchModel &model, const std::vector<int> &residues) {
|
|
SearchModel out = model;
|
|
out.residue_index.clear();
|
|
out.residues_kept = 0;
|
|
size_t i = 0;
|
|
for (gemmi::Chain &ch : out.st.models[0].chains) {
|
|
std::vector<gemmi::Residue> kept;
|
|
for (gemmi::Residue &r : ch.residues) {
|
|
const int index = model.residue_index[i++];
|
|
if (std::binary_search(residues.begin(), residues.end(), index)) {
|
|
kept.push_back(std::move(r));
|
|
out.residue_index.push_back(index);
|
|
}
|
|
}
|
|
ch.residues = std::move(kept);
|
|
out.residues_kept += static_cast<int>(ch.residues.size());
|
|
}
|
|
std::erase_if(out.st.models[0].chains, [](const gemmi::Chain &ch) { return ch.residues.empty(); });
|
|
finish_search_model(out);
|
|
// The domain's residues as ranges of the prediction's numbering, from 1.
|
|
std::string ranges;
|
|
for (size_t k = 0; k < residues.size(); k++) {
|
|
if (k > 0 && residues[k] == residues[k - 1] + 1 && k + 1 < residues.size() && residues[k + 1] == residues[k] + 1)
|
|
continue;
|
|
if (k > 0 && residues[k] == residues[k - 1] + 1)
|
|
ranges += fmt::format("-{}", residues[k] + 1);
|
|
else
|
|
ranges += fmt::format("{}{}", ranges.empty() ? "" : ",", residues[k] + 1);
|
|
}
|
|
out.description = fmt::format("PAE domain of {} residues ({})", out.residues_kept, ranges);
|
|
return out;
|
|
}
|
|
|
|
// --- rotations ---
|
|
|
|
std::vector<gemmi::Mat33> UniformRotations(double step_deg) {
|
|
const double step = step_deg * PI / 180.0;
|
|
const double golden = PI * (3.0 - std::sqrt(5.0));
|
|
const int shells = std::max(1, static_cast<int>(std::lround(PI / step)));
|
|
std::vector<gemmi::Mat33> out{gemmi::Mat33()};
|
|
for (int i = 1; i <= shells; i++) {
|
|
const double kappa = PI * i / shells;
|
|
const double reach = step / (2.0 * std::sin(kappa / 2.0)); // the axis spacing this shell needs
|
|
int n = std::max(1, static_cast<int>(std::lround(4.0 * PI / (reach * reach))));
|
|
// At 180 deg an axis and its opposite are the same rotation: one hemisphere.
|
|
const bool half = i == shells;
|
|
if (half)
|
|
n = (n + 1) / 2;
|
|
for (int k = 0; k < n; k++) {
|
|
const double z = half ? 1.0 - (k + 0.5) / n : 1.0 - (2.0 * k + 1.0) / n;
|
|
const double r = std::sqrt(std::max(0.0, 1.0 - z * z));
|
|
const double x = r * std::cos(golden * k), y = r * std::sin(golden * k);
|
|
const double s = std::sin(kappa), c1 = 1.0 - std::cos(kappa);
|
|
out.push_back({1 + c1 * (x * x - 1), -s * z + c1 * x * y, s * y + c1 * x * z,
|
|
s * z + c1 * x * y, 1 + c1 * (y * y - 1), -s * x + c1 * y * z,
|
|
-s * y + c1 * x * z, s * x + c1 * y * z, 1 + c1 * (z * z - 1)});
|
|
}
|
|
}
|
|
return out;
|
|
}
|
|
|
|
std::vector<gemmi::Mat33> CartesianPointGroup(const gemmi::SpaceGroup &sg, const gemmi::UnitCell &cell) {
|
|
std::vector<gemmi::Mat33> out;
|
|
for (const gemmi::Op &op : sg.operations().sym_ops) {
|
|
gemmi::Mat33 f;
|
|
for (int i = 0; i < 3; i++)
|
|
for (int j = 0; j < 3; j++)
|
|
f.a[i][j] = static_cast<double>(op.rot[i][j]) / gemmi::Op::DEN;
|
|
out.push_back(cell.orth.mat.multiply(f).multiply(cell.frac.mat));
|
|
}
|
|
return out;
|
|
}
|
|
|
|
std::vector<gemmi::Mat33> RotationAsu(const std::vector<gemmi::Mat33> &all,
|
|
const std::vector<gemmi::Mat33> &point_group) {
|
|
std::vector<gemmi::Mat33> out;
|
|
for (const gemmi::Mat33 &r : all) {
|
|
const double t = trace(r);
|
|
if (std::all_of(point_group.begin(), point_group.end(),
|
|
[&](const gemmi::Mat33 &s) { return t >= trace(s.multiply(r)) - 1e-9; }))
|
|
out.push_back(r);
|
|
}
|
|
return out;
|
|
}
|
|
|
|
double RotationDistanceDeg(const gemmi::Mat33 &a, const gemmi::Mat33 &b,
|
|
const std::vector<gemmi::Mat33> &point_group) {
|
|
double best = -1.0;
|
|
for (const gemmi::Mat33 &s : point_group)
|
|
best = std::max(best, trace(a.multiply(s.multiply(b).transpose())));
|
|
return std::acos(std::clamp((best - 1.0) / 2.0, -1.0, 1.0)) * 180.0 / PI;
|
|
}
|
|
|
|
// --- the molecular transform ---
|
|
|
|
MolecularTransform::MolecularTransform(const SearchModel &model, double d_min, size_t nthreads) {
|
|
// On a P1 box six times the model's radius: the transform of a body of radius r turns its phase
|
|
// by pi r / box per grid step, and trilinear interpolation reads it to a few percent at pi/6.
|
|
// Interpolating the molecular transform of one model in place of a Fourier summation per
|
|
// orientation, following Navaza (1994) Acta Cryst. A50, 157-163
|
|
box_ = std::max(6.0 * model.radius_A, 120.0);
|
|
gemmi::Structure st = model.st;
|
|
st.cell = gemmi::UnitCell(box_, box_, box_, 90, 90, 90);
|
|
st.spacegroup_hm = "P 1";
|
|
st.setup_cell_images();
|
|
gemmi::DensityCalculator<Table, float> dc;
|
|
dc.d_min = d_min;
|
|
dc.rate = 1.5;
|
|
dc.set_grid_cell_and_spacegroup(st);
|
|
dc.set_refmac_compatible_blur(st.models[0]);
|
|
PutModelDensityOnGrid(dc, st.models[0], {}, nthreads);
|
|
#ifdef JFJOCH_USE_CUDA
|
|
if (get_gpu_count() > 0 && MapToFPhiGPUFits(dc.grid))
|
|
grid_ = MapToFPhiGPU(dc.grid);
|
|
else
|
|
#endif
|
|
grid_ = MapToFPhi(dc.grid);
|
|
d_min_ = d_min;
|
|
blur_ = dc.blur;
|
|
}
|
|
|
|
std::complex<double> MolecularTransform::At(const gemmi::Vec3 &s) const {
|
|
const double u[3] = {s.x * box_, s.y * box_, s.z * box_};
|
|
const int i0[3] = {static_cast<int>(std::floor(u[0])), static_cast<int>(std::floor(u[1])),
|
|
static_cast<int>(std::floor(u[2]))};
|
|
const double f[3] = {u[0] - i0[0], u[1] - i0[1], u[2] - i0[2]};
|
|
const auto get = [&](int i, int j, int l) {
|
|
bool conj = false;
|
|
if (l < 0) { // the half-l grid holds l >= 0; F(-h) = conj F(h)
|
|
i = -i; j = -j; l = -l;
|
|
conj = true;
|
|
}
|
|
i %= grid_.nu; if (i < 0) i += grid_.nu;
|
|
j %= grid_.nv; if (j < 0) j += grid_.nv;
|
|
const std::complex<float> v = grid_.data[i + static_cast<size_t>(grid_.nu) * (j + static_cast<size_t>(grid_.nv) * l)];
|
|
return conj ? std::conj(std::complex<double>(v)) : std::complex<double>(v);
|
|
};
|
|
std::complex<double> out = 0;
|
|
for (int dx = 0; dx < 2; dx++)
|
|
for (int dy = 0; dy < 2; dy++)
|
|
for (int dz = 0; dz < 2; dz++) {
|
|
const double w = (dx ? f[0] : 1 - f[0]) * (dy ? f[1] : 1 - f[1]) * (dz ? f[2] : 1 - f[2]);
|
|
out += w * get(i0[0] + dx, i0[1] + dy, i0[2] + dz);
|
|
}
|
|
return out * std::exp(blur_ * s.length_sq() / 4.0);
|
|
}
|
|
|
|
// --- the observations ---
|
|
|
|
SearchReflections MakeSearchReflections(const std::vector<MergedReflection> &merged, const gemmi::UnitCell &cell,
|
|
const gemmi::SpaceGroup &sg, double d_min, double d_max) {
|
|
const gemmi::GroupOps gops = sg.operations();
|
|
const gemmi::ReciprocalAsu asu(&sg);
|
|
struct Row { gemmi::Miller h; double i_eps; double s2; double mult; };
|
|
std::vector<Row> rows;
|
|
for (const MergedReflection &r : merged) {
|
|
if (!std::isfinite(r.I) || r.d < d_min || r.d > d_max)
|
|
continue;
|
|
gemmi::Miller h{{r.h, r.k, r.l}};
|
|
if (!asu.is_in(h))
|
|
h = asu.to_asu(h, gops).first;
|
|
if (gops.is_systematically_absent(h))
|
|
continue;
|
|
// The reflections the whole sphere holds for this one: its images under the point group and
|
|
// their Friedel mates, each counted once.
|
|
std::vector<gemmi::Miller> images;
|
|
for (const gemmi::Op &op : gops.sym_ops)
|
|
for (int sign : {1, -1}) {
|
|
gemmi::Miller k = rotated_index(op, h);
|
|
for (int &x : k)
|
|
x *= sign;
|
|
if (std::find(images.begin(), images.end(), k) == images.end())
|
|
images.push_back(k);
|
|
}
|
|
rows.push_back({h, r.I / gops.epsilon_factor_without_centering(h), cell.calculate_1_d2(h),
|
|
static_cast<double>(images.size())});
|
|
}
|
|
std::stable_sort(rows.begin(), rows.end(), [](const Row &a, const Row &b) { return a.s2 < b.s2; });
|
|
SearchReflections out;
|
|
const size_t n = rows.size();
|
|
out.bins = static_cast<int>(std::min<size_t>(SEARCH_BINS, std::max<size_t>(n / 100, 1)));
|
|
std::vector<double> sum(out.bins, 0.0);
|
|
std::vector<int> count(out.bins, 0);
|
|
for (size_t i = 0; i < n; i++) {
|
|
const int b = static_cast<int>(i * out.bins / n);
|
|
out.hkl.push_back(rows[i].h);
|
|
out.multiplicity.push_back(rows[i].mult);
|
|
out.bin.push_back(b);
|
|
sum[b] += rows[i].i_eps;
|
|
count[b]++;
|
|
}
|
|
for (size_t i = 0; i < n; i++) {
|
|
const double mean = sum[out.bin[i]] / count[out.bin[i]];
|
|
out.e2.push_back(mean > 0 ? rows[i].i_eps / mean : 0.0);
|
|
}
|
|
return out;
|
|
}
|
|
|
|
// --- the rotation function ---
|
|
|
|
std::vector<gemmi::Vec3> CopyVectors(const SearchReflections &refl, const gemmi::UnitCell &cell,
|
|
const gemmi::SpaceGroup &sg) {
|
|
const gemmi::GroupOps gops = sg.operations();
|
|
const size_t nops = gops.sym_ops.size();
|
|
const gemmi::Mat33 frac_t = cell.frac.mat.transpose();
|
|
std::vector<gemmi::Vec3> s(refl.hkl.size() * nops);
|
|
for (size_t i = 0; i < refl.hkl.size(); i++)
|
|
for (size_t j = 0; j < nops; j++)
|
|
s[i * nops + j] = frac_t.multiply(miller_vec(rotated_index(gops.sym_ops[j], refl.hkl[i])));
|
|
return s;
|
|
}
|
|
|
|
std::vector<double> RotationFunction(const MolecularTransform &ft, const SearchReflections &refl,
|
|
const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
|
|
const std::vector<gemmi::Mat33> &rotations, size_t nthreads) {
|
|
const size_t n = refl.hkl.size();
|
|
const std::vector<gemmi::Vec3> s = CopyVectors(refl, cell, sg);
|
|
const size_t nops = s.size() / n;
|
|
double w_mean = 0;
|
|
for (double e2 : refl.e2)
|
|
w_mean += e2;
|
|
w_mean /= static_cast<double>(n);
|
|
double w_ss = 0;
|
|
for (double e2 : refl.e2)
|
|
w_ss += (e2 - w_mean) * (e2 - w_mean);
|
|
|
|
std::vector<double> out(rotations.size(), 0.0);
|
|
ParallelFor(static_cast<int>(rotations.size()), nthreads, [&](int r) {
|
|
// The model turned by R has F_R(s) = F_model(R^T s).
|
|
const gemmi::Mat33 rt = rotations[r].transpose();
|
|
std::vector<double> ic(n, 0.0), bin_sum(refl.bins, 0.0);
|
|
std::vector<int> bin_n(refl.bins, 0);
|
|
for (size_t i = 0; i < n; i++) {
|
|
for (size_t j = 0; j < nops; j++)
|
|
ic[i] += std::norm(ft.At(rt.multiply(s[i * nops + j])));
|
|
bin_sum[refl.bin[i]] += ic[i];
|
|
bin_n[refl.bin[i]]++;
|
|
}
|
|
double x_sum = 0, x_ss = 0, xy = 0;
|
|
for (size_t i = 0; i < n; i++) {
|
|
const double x = ic[i] / (bin_sum[refl.bin[i]] / bin_n[refl.bin[i]]);
|
|
x_sum += x;
|
|
x_ss += x * x;
|
|
xy += x * (refl.e2[i] - w_mean);
|
|
}
|
|
const double sxx = x_ss - x_sum * x_sum / static_cast<double>(n);
|
|
out[r] = sxx > 0 && w_ss > 0 ? xy / std::sqrt(sxx * w_ss) : 0.0;
|
|
});
|
|
return out;
|
|
}
|
|
|
|
// --- the translation function ---
|
|
//
|
|
// The full-symmetry translation function as one Fourier synthesis, following Crowther & Blow (1967)
|
|
// Acta Cryst. 23, 544-548: with the model at t, copy j adds G_j(h) exp(2 pi i (h R_j).t) exp(2 pi i h.t_j)
|
|
// to F(h), so |F(h)|^2 summed against E^2 - 1 over the sphere is, for each pair of copies j != k, a
|
|
// Fourier series in t with coefficient (E^2 - 1) g_j conj(g_k), g_j = G_j(h) exp(2 pi i h.t_j), at the
|
|
// index h R_j - h R_k. Every pair's terms go on one grid, and one FFT gives the function at every t.
|
|
// The pairs j = k are the rotation function's, constant in t.
|
|
|
|
TranslationSetup MakeTranslationSetup(const SearchReflections &refl, const gemmi::UnitCell &cell,
|
|
const gemmi::SpaceGroup &sg, bool known_part) {
|
|
const gemmi::GroupOps gops = sg.operations();
|
|
const int nops = static_cast<int>(gops.sym_ops.size());
|
|
const gemmi::Mat33 frac_t = cell.frac.mat.transpose();
|
|
TranslationSetup out;
|
|
out.ops = nops;
|
|
out.bins = refl.bins;
|
|
for (int j = 0; j < nops; j++)
|
|
for (int k = 0; k < nops; k++)
|
|
if (j != k) {
|
|
out.pair_j.push_back(j);
|
|
out.pair_k.push_back(k);
|
|
}
|
|
std::vector<gemmi::Miller> rotated; // terms x ops
|
|
for (size_t i = 0; i < refl.hkl.size(); i++) {
|
|
const gemmi::Miller &h = refl.hkl[i];
|
|
if (gops.is_systematically_absent(h))
|
|
continue;
|
|
for (const gemmi::Op &op : gops.sym_ops) {
|
|
rotated.push_back(rotated_index(op, h));
|
|
out.s.push_back(frac_t.multiply(miller_vec(rotated.back())));
|
|
const double ht = (h[0] * op.tran[0] + h[1] * op.tran[1] + h[2] * op.tran[2])
|
|
/ static_cast<double>(gemmi::Op::DEN);
|
|
out.phase.push_back(std::polar(1.0, 2.0 * PI * ht));
|
|
}
|
|
out.weight.push_back(refl.multiplicity[i] * (refl.e2[i] - 1.0));
|
|
out.e2.push_back(refl.e2[i]);
|
|
out.multiplicity.push_back(refl.multiplicity[i]);
|
|
out.bin.push_back(refl.bin[i]);
|
|
}
|
|
int max_m[3] = {0, 0, 0};
|
|
for (size_t t = 0; t < out.Terms(); t++)
|
|
for (size_t p = 0; p < out.pair_j.size(); p++)
|
|
for (int a = 0; a < 3; a++)
|
|
max_m[a] = std::max(max_m[a], std::abs(rotated[t * nops + out.pair_j[p]][a]
|
|
- rotated[t * nops + out.pair_k[p]][a]));
|
|
if (known_part)
|
|
for (const gemmi::Miller &k : rotated)
|
|
for (int a = 0; a < 3; a++)
|
|
max_m[a] = std::max(max_m[a], std::abs(k[a]));
|
|
out.nu = good_even_size(2 * max_m[0] + 2);
|
|
out.nv = good_even_size(2 * max_m[1] + 2);
|
|
out.nw = good_even_size(2 * max_m[2] + 2);
|
|
for (size_t t = 0; t < out.Terms(); t++)
|
|
for (size_t p = 0; p < out.pair_j.size(); p++) {
|
|
const gemmi::Miller &kj = rotated[t * nops + out.pair_j[p]], &kk = rotated[t * nops + out.pair_k[p]];
|
|
const int l = kj[2] - kk[2];
|
|
out.target_u.push_back(l < 0 ? -1 : ((kj[0] - kk[0]) % out.nu + out.nu) % out.nu);
|
|
out.target_v.push_back(((kj[1] - kk[1]) % out.nv + out.nv) % out.nv);
|
|
out.target_l.push_back(l);
|
|
}
|
|
if (known_part)
|
|
for (size_t t = 0; t < out.Terms(); t++)
|
|
for (int j = 0; j < nops; j++) {
|
|
const gemmi::Miller &k = rotated[t * nops + j];
|
|
const auto add = [&](int sign) {
|
|
out.cross_term.push_back(static_cast<int>(t));
|
|
out.cross_op.push_back(j);
|
|
out.cross_u.push_back(((sign * k[0]) % out.nu + out.nu) % out.nu);
|
|
out.cross_v.push_back(((sign * k[1]) % out.nv + out.nv) % out.nv);
|
|
out.cross_l.push_back(sign * k[2]);
|
|
out.cross_conj.push_back(sign < 0);
|
|
};
|
|
if (k[2] >= 0)
|
|
add(1);
|
|
if (k[2] <= 0)
|
|
add(-1);
|
|
}
|
|
return out;
|
|
}
|
|
|
|
std::vector<std::complex<double>> KnownStructureFactors(const MolecularTransform &ft, const TranslationSetup &setup,
|
|
const gemmi::UnitCell &cell,
|
|
const std::vector<gemmi::Mat33> &rotations,
|
|
const std::vector<gemmi::Fractional> &positions,
|
|
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.
|
|
else if (!(cand.tf_z >= TF_Z_PLACED || (cand.llg_z >= TF_Z_PLACED && cand.copies >= cand.copies_likely)))
|
|
cand.verdict = "NOT_PLACED";
|
|
else
|
|
cand.verdict = "PLACED";
|
|
out.copy_rotation = copy_rotation;
|
|
out.cand = std::move(cand);
|
|
out.orientation = orientation;
|
|
return out;
|
|
};
|
|
// Their time in wall-clock seconds, counted as rescoring (the further copies' translation functions
|
|
// included).
|
|
const auto t_groups = std::chrono::steady_clock::now();
|
|
std::vector<std::future<GroupResult>> running;
|
|
for (size_t g = 0; g < groups.size(); g++)
|
|
running.push_back(std::async(std::launch::async, search_group, g, std::vector<gemmi::Mat33>()));
|
|
std::vector<GroupResult> found;
|
|
for (auto &r : running)
|
|
found.push_back(r.get());
|
|
// With tNCS, every group that placed fewer copies than the best one searched again with the best one's
|
|
// copies' orientations tried first: the groups share the rotations, so a copy the search placed in one is in
|
|
// the same orientation in the others, and one group's search failing - next to copies related by a
|
|
// translation, the first copy alone fits no better than at random, and the copy the search keeps
|
|
// can be the wrong one - is then not taken for evidence against it. The better of the two searches
|
|
// is kept.
|
|
int lead = -1;
|
|
for (size_t g = 0; g < found.size(); g++)
|
|
if (found[g].cand.verdict == "PLACED" && (lead < 0 || found[g].cand.r_free < found[lead].cand.r_free))
|
|
lead = static_cast<int>(g);
|
|
if (lead >= 0 && result.tncs.present) {
|
|
std::vector<std::future<GroupResult>> again(found.size());
|
|
for (size_t g = 0; g < found.size(); g++)
|
|
if (static_cast<int>(g) != lead && found[g].cand.vm >= VM_MIN
|
|
&& found[g].cand.copies < found[lead].cand.copies)
|
|
again[g] = std::async(std::launch::async, search_group, g, found[lead].copy_rotation);
|
|
for (size_t g = 0; g < found.size(); g++)
|
|
if (again[g].valid()) {
|
|
GroupResult r = again[g].get();
|
|
found[g].log.Info("Model prior: {:<12} searched again from {}'s copies: {} copies, R-free {:.3f}, was {:.3f}",
|
|
groups[g].xhm(), groups[lead].xhm(), r.cand.copies, r.cand.r_free,
|
|
found[g].cand.r_free);
|
|
if (std::isfinite(r.cand.r_free) && !(found[g].cand.r_free <= r.cand.r_free)) {
|
|
r.log.ReplayInto(found[g].log);
|
|
r.log = found[g].log;
|
|
found[g] = std::move(r);
|
|
}
|
|
}
|
|
}
|
|
result.t_rescore += seconds_since(t_groups);
|
|
for (GroupResult &gr : found) {
|
|
gr.log.ReplayInto(logger);
|
|
orientations.push_back(gr.orientation);
|
|
result.candidates.push_back(std::move(gr.cand));
|
|
}
|
|
};
|
|
|
|
// --- the decision ---
|
|
// Within one merge the candidates are scored on the same reflections, so their R-free values can
|
|
// be compared directly: the placed candidate with the lowest R-free is chosen where it leads every
|
|
// other candidate of that merge that packs by R_FREE_MARGIN.
|
|
// The null is built for the candidate the choice would fall on, and only for it: a placed candidate
|
|
// that does not beat its own null is set aside, and the choice falls on the next.
|
|
const auto choose_in = [&](int merge) -> int {
|
|
// 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, 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;
|
|
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 (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];
|
|
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);
|
|
}
|
|
}
|