leonarski_fandClaude Opus 5.5
50e3b59fb7
The ENANTIOMORPH hand check ran a whole ValidateAgainstModel on the placed model: model read and frame probe, a rigid body, a 9-replicate null (~0.5 s of the ~0.6-0.9 s), sigma_A, three maps written. The decision needs only --model's rule: the model fits its null (the MR choice has already required its own null of the same construction) and the anomalous map at the model's atoms is not inverted. ReadModelAnomalousMap (ModelValidation) computes F_calc and F_mask to the data's resolution (GPU engine where there is one), carries them through the bulk solvent MR scored the placement with (the overall scale and anisotropic B do not move phases), and reads the anomalous map with the same code as the validation: its coefficient construction and its reading at the atoms are factored out of Validate (anomalous_coefficients, read_anomalous_map) and used by both, so plain --model is unchanged. Hand check 0.6/0.6/0.9 s -> 0.027/0.028/0.040 s on 11if/7ris/9zmu; same ADOPTED/VETOED on all seven cases tried (the three demos, two anomalous lysozyme sets and their mirrored-detector runs), with the scatterer means and deepest troughs within 0.1 sigma of the full validation's. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
1075 lines
50 KiB
C++
1075 lines
50 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 <numeric>
|
|
#include <random>
|
|
#include <unordered_map>
|
|
|
|
#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 "RigidBodyRefine.h"
|
|
#ifdef JFJOCH_USE_CUDA
|
|
#include "MolecularReplacementGPU.h"
|
|
#include "RigidBodyGPU.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;
|
|
// 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 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; at or above VM_TWO_COPIES two copies would.
|
|
constexpr double VM_MIN = 1.5;
|
|
constexpr double VM_TWO_COPIES = 4.0;
|
|
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;
|
|
constexpr double NULL_SIGMA_MIN = 3.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)};
|
|
}
|
|
|
|
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;
|
|
};
|
|
|
|
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, scaled to the working amplitudes as model
|
|
// validation scales them (ModelValidation.cpp, fit_model), and the R-factors over the reflections
|
|
// to d_min.
|
|
PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGroup &sg,
|
|
const std::vector<MergedReflection> &merged, double d_min, 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);
|
|
gemmi::AsuData<std::complex<float>> fmask = MapToFPhi(mask).prepare_asu_data(d_min, 0);
|
|
PutModelDensityOnGrid(dc, st.models[0], leaders, nthreads);
|
|
gemmi::AsuData<std::complex<float>> fmodel = MapToFPhi(dc.grid).prepare_asu_data(d_min, dc.blur, false, false, false);
|
|
|
|
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);
|
|
FitModelScale(scaling, {.b_hi = SolventBMax(st.models[0])}, nthreads);
|
|
scaling.scale_data(fmodel, &fmask);
|
|
|
|
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;
|
|
}
|
|
|
|
// 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.
|
|
PlacementScore refine_placement(gemmi::Structure &st, const gemmi::SpaceGroup &sg,
|
|
const std::vector<MergedReflection> &merged, double d_min, RigidBodyGPUPool *pool,
|
|
size_t nthreads) {
|
|
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);
|
|
if (!rb.converged)
|
|
SetModelPositions(st.models[0], start);
|
|
return score_placement(st, sg, merged, d_min, nthreads);
|
|
}
|
|
|
|
// 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.
|
|
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) {
|
|
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).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<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());
|
|
for (gemmi::Position &p : pos)
|
|
p = gemmi::Position(draws[i].multiply(gemmi::Vec3(p) - c) + c);
|
|
SetModelPositions(rep.models[0], pos);
|
|
r_work[i] = refine_placement(rep, sg, coarse, null_d_min, pool, nthreads).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;
|
|
}
|
|
|
|
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 ---
|
|
|
|
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) {
|
|
for (const 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;
|
|
}
|
|
out.residues_read += static_cast<int>(ch.residues.size());
|
|
// 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.
|
|
std::erase_if(ch.residues, [](const gemmi::Residue &r) { return r.atoms.empty() || r.atoms[0].b_iso < 70.0f; });
|
|
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
|
|
std::vector<gemmi::Position> pos;
|
|
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);
|
|
pos.push_back(a.pos);
|
|
}
|
|
gemmi::Vec3 c;
|
|
for (const gemmi::Position &p : pos)
|
|
c += p;
|
|
c *= 1.0 / static_cast<double>(pos.size());
|
|
for (gemmi::Position &p : pos) {
|
|
p = gemmi::Position(gemmi::Vec3(p) - c);
|
|
out.radius_A = std::max(out.radius_A, p.length());
|
|
}
|
|
SetModelPositions(model, pos);
|
|
out.mass_da = DALTON_PER_RESIDUE * out.residues_kept;
|
|
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);
|
|
grid_ = MapToFPhi(dc.grid);
|
|
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) {
|
|
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.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]));
|
|
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);
|
|
}
|
|
return out;
|
|
}
|
|
|
|
TranslationPeak TranslationFunction(const MolecularTransform &ft, const TranslationSetup &setup,
|
|
const gemmi::Mat33 &rotation) {
|
|
TranslationPeak out;
|
|
const size_t terms = setup.Terms(), nops = setup.ops, npairs = setup.pair_j.size();
|
|
if (nops < 2 || terms == 0)
|
|
return out; // 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]]);
|
|
}
|
|
}
|
|
// 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;
|
|
}
|
|
|
|
// --- 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) {
|
|
const gemmi::UnitCell &cell = placed.cell;
|
|
const gemmi::SpaceGroup *sg = placed.find_spacegroup();
|
|
std::vector<gemmi::Fractional> ca;
|
|
for (const gemmi::Chain &ch : placed.models[0].chains)
|
|
for (const gemmi::Residue &r : ch.residues)
|
|
for (const gemmi::Atom &a : r.atoms)
|
|
if (a.name == "CA")
|
|
ca.push_back(cell.fractionalize(a.pos));
|
|
if (ca.empty() || sg == nullptr)
|
|
return 0.0;
|
|
gemmi::Fractional c0;
|
|
for (const gemmi::Fractional &f : ca)
|
|
c0 = c0 + f;
|
|
c0 = gemmi::Fractional(c0 * (1.0 / static_cast<double>(ca.size())));
|
|
double radius = 0;
|
|
for (const gemmi::Fractional &f : ca)
|
|
radius = std::max(radius, cell.orthogonalize_difference(f - c0).length());
|
|
constexpr double CLASH_A = 3.0;
|
|
std::vector<char> clash(ca.size(), 0);
|
|
for (const gemmi::Op &op : sg->operations().all_ops_sorted()) {
|
|
std::vector<gemmi::Fractional> mate;
|
|
for (const gemmi::Fractional &f : 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())));
|
|
const gemmi::Fractional nearest(std::round(c0.x - cm.x), std::round(c0.y - cm.y), std::round(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 (op == gemmi::Op::identity() && shift.x == 0 && shift.y == 0 && shift.z == 0)
|
|
continue;
|
|
if (cell.orthogonalize_difference(cm + shift - c0).length() > 2 * radius + CLASH_A)
|
|
continue;
|
|
for (size_t a = 0; a < ca.size(); a++) {
|
|
if (clash[a])
|
|
continue;
|
|
for (const gemmi::Fractional &m : mate)
|
|
if (cell.orthogonalize_difference(m + shift - ca[a]).length_sq() < CLASH_A * CLASH_A) {
|
|
clash[a] = 1;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return static_cast<double>(std::count(clash.begin(), clash.end(), 1)) / static_cast<double>(ca.size());
|
|
}
|
|
|
|
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) {
|
|
if (g.number == sg.number && g.xhm() == sg.xhm())
|
|
continue;
|
|
const gemmi::GroupOps gops = g.operations();
|
|
if (!same_rotations_and_centring(ops, gops) || gemmi::find_spacegroup_by_ops(gops) != &g)
|
|
continue;
|
|
out.push_back(g);
|
|
}
|
|
return out;
|
|
}
|
|
|
|
// --- the search and the decision ---
|
|
|
|
MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector<MRMerge> &merges,
|
|
Logger &logger, size_t nthreads) {
|
|
MolecularReplacementResult result;
|
|
result.model_path = model_path;
|
|
auto t0 = std::chrono::steady_clock::now();
|
|
const SearchModel model = PrepareSearchModel(model_path);
|
|
result.residues_read = model.residues_read;
|
|
result.residues_kept = model.residues_kept;
|
|
result.mass_da = model.mass_da;
|
|
if (!model.failure_reason.empty()) {
|
|
result.failure_reason = model.failure_reason;
|
|
return result;
|
|
}
|
|
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);
|
|
result.t_model = seconds_since(t0);
|
|
logger.Info("Model prior: {} of {} residues kept at pLDDT >= 70 ({:.1f} kDa, radius {:.1f} A)",
|
|
model.residues_kept, model.residues_read, model.mass_da / 1000.0, model.radius_A);
|
|
|
|
t0 = std::chrono::steady_clock::now();
|
|
const MolecularTransform ft(model, result.d_min_search, nthreads);
|
|
#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);
|
|
const std::vector<gemmi::Mat33> sphere = UniformRotations(ROTATION_STEP_DEG);
|
|
|
|
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);
|
|
#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;
|
|
}
|
|
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 ? (rf[peaks.front()] - rf_mean) / rf_sd : 0.0);
|
|
|
|
// --- translation, for every group and every peak orientation ---
|
|
t0 = std::chrono::steady_clock::now();
|
|
const std::vector<gemmi::SpaceGroup> groups = SameRotationGroups(merge.space_group);
|
|
std::vector<gemmi::Mat33> peak_rotations;
|
|
for (size_t p : peaks)
|
|
peak_rotations.push_back(rotations[p]);
|
|
std::vector<TranslationPeak> tf_all;
|
|
for (const gemmi::SpaceGroup &sg : groups) {
|
|
const TranslationSetup setup = MakeTranslationSetup(refl, cell, sg);
|
|
std::vector<TranslationPeak> tf(peaks.size());
|
|
#ifdef JFJOCH_USE_CUDA
|
|
if (gpu)
|
|
tf = gpu->TranslationFunction(setup, peak_rotations);
|
|
else
|
|
#endif
|
|
ParallelFor(static_cast<int>(peaks.size()), nthreads, [&](int i) {
|
|
tf[i] = TranslationFunction(ft, setup, peak_rotations[i]);
|
|
});
|
|
tf_all.insert(tf_all.end(), tf.begin(), tf.end());
|
|
}
|
|
result.t_translation += seconds_since(t0);
|
|
|
|
for (size_t g = 0; g < groups.size(); g++) {
|
|
const gemmi::SpaceGroup &sg = groups[g];
|
|
MRCandidate cand;
|
|
cand.merge = static_cast<int>(m);
|
|
cand.space_group = sg;
|
|
const gemmi::GroupOps gops = sg.operations();
|
|
cand.vm = cell.volume / static_cast<double>(gops.order()) / model.mass_da;
|
|
cand.copies_likely = cand.vm >= VM_TWO_COPIES ? 2 : 1;
|
|
const TranslationPeak *tf = &tf_all[g * 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 tf[a].z > tf[b].z; });
|
|
|
|
// --- the best few through packing, the rigid body and R-free ---
|
|
t0 = std::chrono::steady_clock::now();
|
|
RigidBodyEngines engines;
|
|
bool have = false;
|
|
gemmi::Mat33 orientation;
|
|
// Where one copy cannot be the content of the cell - too big for it, or small enough for two -
|
|
// the group is not decided on here, and only its best translation-function placement is
|
|
// checked for packing.
|
|
const bool one_copy = cand.vm >= VM_MIN && cand.copies_likely == 1;
|
|
for (int k = 0; k < std::min<int>(one_copy ? RESCORED_PER_GROUP : 1, static_cast<int>(by_tf.size())); k++) {
|
|
const size_t i = by_tf[k];
|
|
auto placed = std::make_shared<gemmi::Structure>(
|
|
PlaceSearchModel(model, rotations[peaks[i]], tf[i].t, cell, sg));
|
|
// 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.
|
|
double clash = ClashFraction(*placed);
|
|
PlacementScore score{NAN, NAN, NAN, NAN};
|
|
if (one_copy && clash <= 2 * CLASH_MAX) {
|
|
if (engines.get() == nullptr)
|
|
engines = make_engines(*placed, sg, *merge.merged, d_min, 1, logger);
|
|
score = refine_placement(*placed, sg, *merge.merged, d_min, engines.get(), nthreads);
|
|
clash = ClashFraction(*placed);
|
|
}
|
|
// 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 ? (rf[peaks[i]] - rf_mean) / rf_sd : 0.0;
|
|
cand.tf_z = tf[i].z;
|
|
cand.clash_fraction = clash;
|
|
cand.r_work = score.r_work;
|
|
cand.r_free = score.r_free;
|
|
cand.k_sol = score.k_sol;
|
|
cand.b_sol = score.b_sol;
|
|
cand.placed = placed;
|
|
orientation = rotations[peaks[i]];
|
|
}
|
|
logger.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 ? (rf[peaks[i]] - rf_mean) / rf_sd : 0.0, tf[i].z, 100 * clash,
|
|
score.r_work, score.r_free);
|
|
}
|
|
result.t_rescore += seconds_since(t0);
|
|
|
|
if (cand.vm < VM_MIN)
|
|
cand.verdict = "MATTHEWS";
|
|
else if (cand.copies_likely > 1)
|
|
cand.verdict = "MULTI_COPY";
|
|
else if (cand.clash_fraction > CLASH_MAX)
|
|
cand.verdict = "CLASH";
|
|
else if (!(cand.tf_z >= TF_Z_PLACED))
|
|
cand.verdict = "NOT_PLACED";
|
|
else
|
|
cand.verdict = "PLACED";
|
|
orientations.push_back(orientation);
|
|
result.candidates.push_back(std::move(cand));
|
|
}
|
|
};
|
|
|
|
// --- the decision ---
|
|
// Within one merge the candidates are scored on the same reflections, so their R-free values can
|
|
// be compared directly: the placed candidate with the lowest R-free is chosen where it leads every
|
|
// other candidate of that merge that packs by R_FREE_MARGIN.
|
|
// The null is built for the candidate the choice would fall on, and only for it: a placed candidate
|
|
// that does not beat its own null is set aside, and the choice falls on the next.
|
|
const auto choose_in = [&](int merge) -> int {
|
|
int best = -1;
|
|
while (true) {
|
|
best = -1;
|
|
for (size_t i = 0; i < result.candidates.size(); i++) {
|
|
const MRCandidate &c = result.candidates[i];
|
|
if (c.merge == merge && c.verdict == "PLACED"
|
|
&& (best < 0 || c.r_free < result.candidates[best].r_free))
|
|
best = static_cast<int>(i);
|
|
}
|
|
if (best < 0)
|
|
return -1;
|
|
MRCandidate &c = result.candidates[best];
|
|
const auto t_null = std::chrono::steady_clock::now();
|
|
RigidBodyEngines engines = make_engines(*c.placed, c.space_group, *merges[merge].merged,
|
|
result.d_min_score, NULL_REPLICATES, logger);
|
|
c.null_sigma = null_sigma(*c.placed, orientations[best], c.space_group, *merges[merge].merged,
|
|
result.d_min_score, engines.get(), nthreads);
|
|
result.t_null += seconds_since(t_null);
|
|
if (c.null_sigma >= NULL_SIGMA_MIN)
|
|
break;
|
|
c.verdict = "NULL";
|
|
}
|
|
for (size_t i = 0; i < result.candidates.size(); i++) {
|
|
const MRCandidate &c = result.candidates[i];
|
|
if (static_cast<int>(i) == best || c.merge != merge || c.clash_fraction > CLASH_MAX)
|
|
continue;
|
|
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";
|
|
}
|
|
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;
|
|
}
|