From f800669839d416c47d3968069f64bdc97b05761c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 11 Oct 2026 20:50:47 +0200 Subject: [PATCH 1/3] jfjoch_mr_replay: --af-model's molecular replacement re-run on a run's merged MTZ files A developer utility, not installed: the search, the decision and the hand check run on the answer's p.mtz (and the twin hypothesis' p_twin_.mtz) that a rugnux run wrote, so a survey of the search over every dataset with a predicted model takes minutes instead of a reprocessing per set. The hand check of an ENANTIOMORPH decision moves from the rugnux driver into CheckEnantiomorphHand, which both call; unchanged. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- .../MolecularReplacement.cpp | 35 ++++++ .../MolecularReplacement.h | 8 ++ rugnux/RugnuxScaleMerge.cpp | 44 +------ tools/CMakeLists.txt | 5 + tools/jfjoch_mr_replay.cpp | 109 ++++++++++++++++++ 5 files changed, 160 insertions(+), 41 deletions(-) create mode 100644 tools/jfjoch_mr_replay.cpp diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp index 72ed74c6c..dc2556006 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.cpp +++ b/image_analysis/structure_refinement/MolecularReplacement.cpp @@ -25,6 +25,7 @@ #include "ModelFFT.h" #include "ModelGrid.h" #include "ModelScaling.h" +#include "ModelValidation.h" #include "RigidBodyRefine.h" #ifdef JFJOCH_USE_CUDA #include "MolecularReplacementGPU.h" @@ -1091,3 +1092,37 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path result.t_null); return result; } + +void CheckEnantiomorphHand(MolecularReplacementResult &result, const std::vector &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); + } +} diff --git a/image_analysis/structure_refinement/MolecularReplacement.h b/image_analysis/structure_refinement/MolecularReplacement.h index fa5031d17..c20a6e44d 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.h +++ b/image_analysis/structure_refinement/MolecularReplacement.h @@ -207,3 +207,11 @@ struct MolecularReplacementResult { // answer's groups can hold the model at all (Matthews, packing) and a group of the subgroup can. MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector &merges, Logger &logger, size_t nthreads); + +// The hand check of an ENANTIOMORPH decision, filled into result.hand_check: taken into the answer by +// --model's rule for a model's enantiomorph - the model has to fit against its own null, which the choice +// already required of it, and the anomalous map at its atoms (`merged`, the answer's merge in +// `data_space_group`) must not come out inverted, read by model validation's own code on the placement +// and the bulk solvent MR scored. Where it is not taken, the decision becomes ENANTIOMORPH_VETOED. +void CheckEnantiomorphHand(MolecularReplacementResult &result, const std::vector &merged, + const gemmi::SpaceGroup &data_space_group, Logger &logger, size_t nthreads); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index a3391b14d..07c9095b4 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -4000,47 +4000,9 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { "The predicted model was not placed: " + mr.failure_reason}); if (mr.chosen >= 0) { const MRCandidate &c = mr.candidates[mr.chosen]; - if (mr.decision == "ENANTIOMORPH") { - // Taken into the answer by --model's rule for a model's enantiomorph: the model has to fit - // against its own null - which the MR choice already required of it - and the anomalous map - // at its atoms must not come out inverted, read by model validation's own code on the - // placement and the bulk solvent MR scored. - const auto t_hand = std::chrono::steady_clock::now(); - Logger quiet = Logger::Buffered(); - const auto data_sg = experiment_.GetSpaceGroupOrP1(); - ModelValidationResult v; - try { - v = ReadModelAnomalousMap(sm.merged, *c.placed, data_sg, c.k_sol, c.b_sol, quiet, - static_cast(config_.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()); - } - v.r_work_sigma = c.null_sigma; - v.model_fits = c.null_sigma >= 3.0; - v.adopted_model_enantiomorph = v.ok && v.model_fits && !v.anomalous_hands_disagree; - auto &h = mr.hand_check; - h.performed = true; - h.adopted = v.adopted_model_enantiomorph; - h.model_fits = v.model_fits; - h.fit_sigma = v.r_work_sigma; - 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 = std::chrono::duration(std::chrono::steady_clock::now() - t_hand).count(); - if (!h.adopted) { - mr.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); - } - } + if (mr.decision == "ENANTIOMORPH") + CheckEnantiomorphHand(mr, sm.merged, experiment_.GetSpaceGroupOrP1(), logger, + static_cast(config_.nthreads)); if (mr.decision == "ENANTIOMORPH") { // A change of label only (see ValidateAgainstModel): the two groups index identically. experiment_.SetSpaceGroup(c.space_group); diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index b6754ec31..d8744bb51 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -12,6 +12,11 @@ TARGET_LINK_LIBRARIES(jfjoch_extract_hkl JFJochReader) # rugnux, the single offline analysis CLI, lives in rugnux/ next to its library. +# rugnux --af-model's molecular replacement re-run on a run's merged MTZ files (developer utility, not +# installed): a survey of the search over many datasets in minutes instead of a reprocessing each. +ADD_EXECUTABLE(jfjoch_mr_replay jfjoch_mr_replay.cpp) +TARGET_LINK_LIBRARIES(jfjoch_mr_replay JFJochStructureRefinement) + # In-place re-compress /entry/data/data of a _data file from bitshuffle/LZ4 to bitshuffle/zstd. ADD_EXECUTABLE(jfjoch_recompress jfjoch_recompress.cpp) TARGET_LINK_LIBRARIES(jfjoch_recompress JFJochHDF5Wrappers) diff --git a/tools/jfjoch_mr_replay.cpp b/tools/jfjoch_mr_replay.cpp new file mode 100644 index 000000000..2d751dcea --- /dev/null +++ b/tools/jfjoch_mr_replay.cpp @@ -0,0 +1,109 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +// Developer utility: rugnux --af-model's molecular replacement re-run on the merged reflections a rugnux +// run wrote, without reprocessing the images - the answer's p.mtz and, where the run made one, the twin +// hypothesis' p_twin_.mtz. The search, the decision and the hand check are rugnux's own code; the +// output is the MODEL PRIOR keys and table of the rugnux report. +// +// jfjoch_mr_replay [ ] +// +// The answer group is the one the run searched in: where it took the model's enantiomorph into its +// label, the file carries the relabelled group and the original has to be given here. + +#include +#include +#include +#include + +#include + +#include "../common/Logger.h" +#include "../image_analysis/structure_refinement/MolecularReplacement.h" + +namespace { + +std::vector read_merged(const std::string &path, UnitCell &cell) { + gemmi::Mtz mtz; + mtz.read_file(path); + const gemmi::UnitCell &c = mtz.cell; + cell = UnitCell(c.a, c.b, c.c, c.alpha, c.beta, c.gamma); + const auto col = [&](const char *label) -> int { + const gemmi::Mtz::Column *x = mtz.column_with_label(label); + return x ? static_cast(x->idx) : -1; + }; + const int i_mean = col("IMEAN"), s_mean = col("SIGIMEAN"), ip = col("I(+)"), sip = col("SIGI(+)"), + im = col("I(-)"), sim = col("SIGI(-)"), f = col("F"), sf = col("SIGF"), fp = col("F(+)"), + sfp = col("SIGF(+)"), fm = col("F(-)"), sfm = col("SIGF(-)"), free = col("FreeR_flag"); + const size_t nc = mtz.columns.size(); + const auto get = [&](size_t row, int i) { return i < 0 ? NAN : mtz.data[row * nc + i]; }; + std::vector out; + for (size_t row = 0; row < static_cast(mtz.nreflections); row++) { + MergedReflection r; + r.h = static_cast(get(row, 0)); + r.k = static_cast(get(row, 1)); + r.l = static_cast(get(row, 2)); + r.I = get(row, i_mean); + r.sigma = get(row, s_mean); + r.I_plus = get(row, ip); + r.sigma_plus = get(row, sip); + r.I_minus = get(row, im); + r.sigma_minus = get(row, sim); + r.F = get(row, f); + r.sigmaF = get(row, sf); + r.F_plus = get(row, fp); + r.sigmaF_plus = get(row, sfp); + r.F_minus = get(row, fm); + r.sigmaF_minus = get(row, sfm); + r.rfree_flag = get(row, free) == 0.0f; + r.d = static_cast(1.0 / std::sqrt(c.calculate_1_d2(gemmi::Miller{{r.h, r.k, r.l}}))); + out.push_back(r); + } + return out; +} + +} // namespace + +int main(int argc, char **argv) { + if (argc != 4 && argc != 6) { + std::cerr << "Usage: jfjoch_mr_replay [ ]\n"; + return 1; + } + Logger logger("jfjoch_mr_replay"); + const gemmi::SpaceGroup *answer_sg = gemmi::find_spacegroup_by_name(argv[2]); + UnitCell cell{}, twin_cell{}; + const std::vector merged = read_merged(argv[3], cell); + std::vector twin; + std::vector merges{{"MERGE", &merged, cell, *answer_sg}}; + if (argc == 6) { + twin = read_merged(argv[5], twin_cell); + merges.push_back({"TWIN_SUBGROUP_MERGE", &twin, twin_cell, *gemmi::find_spacegroup_by_name(argv[4])}); + } + const size_t nthreads = std::max(1u, std::thread::hardware_concurrency()); + MolecularReplacementResult mr = RunMolecularReplacement(argv[1], merges, logger, nthreads); + if (mr.chosen >= 0 && mr.decision == "ENANTIOMORPH") + CheckEnantiomorphHand(mr, merged, *answer_sg, logger, nthreads); + + std::cout << "MODEL_PRIOR= " << (mr.ok ? "PERFORMED" : "NOT_PERFORMED") << "\n"; + if (!mr.ok) { + std::cout << "MODEL_PRIOR_REASON= " << mr.failure_reason << "\n"; + return 0; + } + std::cout << "MODEL_PRIOR_DECISION= " << mr.decision << "\n"; + if (mr.chosen >= 0) + std::cout << "MODEL_PRIOR_SPACE_GROUP= " << mr.candidates[mr.chosen].space_group.xhm() << "\n"; + if (mr.hand_check.performed) + std::cout << "MODEL_PRIOR_HAND_CHECK= " << (mr.hand_check.adopted ? "ADOPTED" : "VETOED") << "\n"; + std::cout << fmt::format("MODEL_PRIOR_SECONDS= {:.2f}\n", mr.t_model + mr.t_transform + mr.t_rotation + + mr.t_translation + mr.t_rescore + mr.t_null + + mr.hand_check.seconds); + for (const MRCandidate &c : mr.candidates) { + const auto r = [](double v) { return std::isfinite(v) ? fmt::format("{:.3f}", v) : std::string("-"); }; + std::cout << fmt::format( + "# {:<19s} {:<11s} {:5.2f} {:5.1f} {:5.1f} {:4.1f}% {:>6s} {:>6s} {:>10s} {}\n", + c.merge == 0 ? "MERGE" : "TWIN_SUBGROUP_MERGE", c.space_group.xhm(), c.vm, c.rf_z, c.tf_z, + 100 * c.clash_fraction, r(c.r_work), r(c.r_free), + std::isfinite(c.null_sigma) ? fmt::format("{:+.1f}", c.null_sigma) : std::string("-"), c.verdict); + } + return 0; +} From dcb01100e7aa13068257f1d8aa094022057d2868 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 12 Oct 2026 00:26:49 +0200 Subject: [PATCH 2/3] RefineRigidBody: a body refined beside fixed copies of it An optional `fixed` model: other atoms of the crystal held where they are - the other copies of a molecular-replacement search - whose Fcalc is added to the body's at every evaluation and whose atoms the bulk-solvent mask covers together with the body's (computed once per zone, with the body where the zone starts, and held). With them the origin is no longer free, so the gauge projection is skipped. On the CPU the fixed part is gridded once per zone; on the GPU the fixed atoms must be copies of the body, and each one's placement is read off by Kabsch superposition and composed on the device with the engine's own density, transform and mask (SetFixedCopies), so a fit beside three copies of a 90 kDa model costs what a fit alone does. Without `fixed` nothing changes (the composition kernel's extra term is null). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- .../structure_refinement/RigidBodyGPU.cpp | 42 +++++++++++- .../structure_refinement/RigidBodyGPU.cu | 64 +++++++++++++++++-- .../structure_refinement/RigidBodyGPUEngine.h | 5 ++ .../structure_refinement/RigidBodyRefine.cpp | 42 +++++++++++- .../structure_refinement/RigidBodyRefine.h | 20 +++++- 5 files changed, 163 insertions(+), 10 deletions(-) diff --git a/image_analysis/structure_refinement/RigidBodyGPU.cpp b/image_analysis/structure_refinement/RigidBodyGPU.cpp index 5316f580b..99509714b 100644 --- a/image_analysis/structure_refinement/RigidBodyGPU.cpp +++ b/image_analysis/structure_refinement/RigidBodyGPU.cpp @@ -280,6 +280,10 @@ void RigidBodyTargetGPU::Placement(const double q[6], double rotation[9], double } void RigidBodyTargetGPU::SetZone(const gemmi::AsuData> &fobs, double d_min) { + // The body where this zone starts, for the mask it shares with fixed atoms: the last evaluation's. + std::array q_start{}; + if (have_point_) + q_start = q_; d_min_ = d_min; solvent_fitted_ = false; host_scale_ = false; @@ -287,6 +291,42 @@ void RigidBodyTargetGPU::SetZone(const gemmi::AsuData> zone_ = pool_.ObservedZone(model_, cell_, sg_, fobs, d_min); try { engine_.SetZone(*zone_); + if (fixed != nullptr) { + // The fixed atoms are copies of the body (RefineRigidBody): each one's placement of the body's + // atoms, by the least-squares superposition of the two. + // Following Kabsch (1976) Acta Cryst. A32, 922-923 + const std::vector all = ModelPositions(*fixed); + const size_t n = base_.size(); + std::vector> copies; + for (size_t c = 0; c + n <= all.size(); c += n) { + Eigen::Vector3d mean = Eigen::Vector3d::Zero(); + for (size_t i = 0; i < n; i++) + mean += Eigen::Vector3d(all[c + i].x, all[c + i].y, all[c + i].z); + mean /= static_cast(n); + Eigen::Matrix3d h = Eigen::Matrix3d::Zero(); + for (size_t i = 0; i < n; i++) { + const Eigen::Vector3d p(base_[i].x - centre_.x, base_[i].y - centre_.y, base_[i].z - centre_.z); + const Eigen::Vector3d q(all[c + i].x - mean.x(), all[c + i].y - mean.y(), all[c + i].z - mean.z()); + h += p * q.transpose(); + } + const Eigen::JacobiSVD svd(h, Eigen::ComputeFullU | Eigen::ComputeFullV); + Eigen::Matrix3d d = Eigen::Matrix3d::Identity(); + d(2, 2) = (svd.matrixV() * svd.matrixU().transpose()).determinant() < 0 ? -1.0 : 1.0; + const Eigen::Matrix3d r = svd.matrixV() * d * svd.matrixU().transpose(); + // The body's relative coordinates have their mean at zero, so the copy's centroid is t. + std::array x{}; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + x[3 * i + j] = r(i, j); + x[9] = mean.x(); + x[10] = mean.y(); + x[11] = mean.z(); + copies.push_back(x); + } + std::array moving{}; + Placement(q_start.data(), moving.data(), moving.data() + 9); + engine_.SetFixedCopies(copies, moving); + } } catch (const JFJochException &e) { throw RigidBodyGPUFailure(e.what()); } @@ -331,7 +371,7 @@ bool RigidBodyTargetGPU::Residuals(const double q[6], double *residuals) { if (zone_->row_hkl.empty()) return false; engine_.Fcalc(rotation, translation); - if (!(hold_mask && have_point_)) + if (fixed == nullptr && !(hold_mask && have_point_)) engine_.Fmask(); if (zone_->point_obs.empty()) diff --git a/image_analysis/structure_refinement/RigidBodyGPU.cu b/image_analysis/structure_refinement/RigidBodyGPU.cu index b8e316280..b5a40a784 100644 --- a/image_analysis/structure_refinement/RigidBodyGPU.cu +++ b/image_analysis/structure_refinement/RigidBodyGPU.cu @@ -106,7 +106,7 @@ __global__ void place_kernel(const double *rel, int n, Placement p, int nu, int // dF/dt where `df_dt` is not null. `f1` is cuFFT's forward transform, which carries exp(-2 pi i h.x): // conjugated and scaled by V/N it is gemmi's F1. __global__ void compose_kernel(const float2 *f1, float norm, const RigidBodyGPUTerm *terms, const double *row_scale, - int rows, int ops, double2 *f, double2 *df_dt) { + int rows, int ops, const double2 *fixed, double2 *f, double2 *df_dt) { const int m = blockIdx.x * blockDim.x + threadIdx.x; if (m >= rows) return; @@ -126,7 +126,8 @@ __global__ void compose_kernel(const float2 *f1, float norm, const RigidBodyGPUT } } const double scale = row_scale[m]; - f[m] = make_double2(scale * sr, scale * si); + f[m] = fixed != nullptr ? make_double2(scale * sr + fixed[m].x, scale * si + fixed[m].y) + : make_double2(scale * sr, scale * si); if (df_dt != nullptr) { // times 2 pi i const double two_pi = 2 * 3.14159265358979323846; @@ -135,6 +136,12 @@ __global__ void compose_kernel(const float2 *f1, float norm, const RigidBodyGPUT } } +__global__ void add_kernel(const double2 *a, int n, double2 *sum) { + const int m = blockIdx.x * blockDim.x + threadIdx.x; + if (m < n) + sum[m] = make_double2(sum[m].x + a[m].x, sum[m].y + a[m].y); +} + // Fmask at each row, read at h directly (the mask is symmetric): gemmi's get_value_by_hkl() on the // transform of the mask, taking the Friedel mate where h is outside the stored half. __global__ void fmask_kernel(const float2 *fm, float norm, const int *row_hkl, int rows, int nu, int nv, int nw, @@ -324,6 +331,8 @@ struct RigidBodyGPUEngineImpl { CudaDevicePtr mask_slot; CudaDevicePtr mask_radius; CudaDevicePtr mask_atoms; // 4 per mask atom: fractional x, y, z in [0,1), radius + CudaDevicePtr copies_mask_atoms; // the same for the body and fixed copies of it, grown as needed + size_t copies_mask_capacity = 0; CudaDevicePtr grid; // 3 grids back to back CudaDevicePtr spectrum; // their 3 transforms @@ -332,7 +341,8 @@ struct RigidBodyGPUEngineImpl { CudaDevicePtr terms; CudaDevicePtr row_hkl; CudaDevicePtr row_scale, row_stol2; - CudaDevicePtr f, df_dt, shifted; + CudaDevicePtr f, df_dt, shifted, fixed; + bool has_fixed = false; // the zone's fixed Fcalc is in `fixed` CudaDevicePtr fmask; CudaDevicePtr obs_row, point_obs; @@ -375,6 +385,7 @@ struct RigidBodyGPUEngineImpl { f = CudaDevicePtr(nr, sync); df_dt = CudaDevicePtr(3 * nr, sync); shifted = CudaDevicePtr(3 * nr, sync); + fixed = CudaDevicePtr(nr, sync); fmask = CudaDevicePtr(nr, sync); obs_row = CudaDevicePtr(no, sync); point_obs = CudaDevicePtr(no, sync); @@ -443,7 +454,7 @@ struct RigidBodyGPUEngineImpl { void Compose(int slot, double2 *out, double2 *derivative) { compose_kernel<<>>(spectrum.get() + slot * ComplexPoints(), Norm(), terms, row_scale, static_cast(rows), static_cast(ops), - out, derivative); + has_fixed ? fixed.get() : nullptr, out, derivative); cuda_err(cudaGetLastError()); } @@ -500,7 +511,7 @@ size_t RigidBodyGPUEngine::DeviceBytes(const RigidBodyGPUCapacity &c) { bytes += 3 * c.grid_points * sizeof(float) + 3 * c.complex_points * sizeof(float2) + c.fft_work_bytes; bytes += ModelDensityGPU::DeviceBytes(c.atoms, c.pairs, c.bricks); bytes += c.terms * sizeof(RigidBodyGPUTerm); - bytes += nr * (3 * sizeof(int) + 2 * sizeof(double) + 7 * sizeof(double2) + sizeof(float2)); + bytes += nr * (3 * sizeof(int) + 2 * sizeof(double) + 8 * sizeof(double2) + sizeof(float2)); bytes += no * (2 * sizeof(int) + sizeof(float) + (1 + 6 + MAX_SCALE_PARAMS) * sizeof(double) + 2 * sizeof(float2)); bytes += ModelMaskGPU::DeviceBytes(c.grid_points) + ModelScaleGPU::DeviceBytes(no); return bytes; @@ -559,6 +570,7 @@ void RigidBodyGPUEngine::SetZone(const RigidBodyGPUZone &zone) { e.n_constraints = zone.constraints.size(); e.constraints = zone.constraints; e.f_mean = zone.f_mean; + e.has_fixed = false; std::vector slot(zone.atoms.size(), -1); std::vector radius(zone.atoms.size(), 0.0f); @@ -607,6 +619,48 @@ void RigidBodyGPUEngine::Fcalc(const double rotation[9], const double translatio e.Compose(0, e.f, e.df_dt); } +void RigidBodyGPUEngine::SetFixedCopies(const std::vector> &fixed, + const std::array &moving) { + RigidBodyGPUEngineImpl &e = *impl_; + cuda_err(cudaSetDevice(e.device)); + const size_t copies = fixed.size() + 1; + if (e.copies_mask_capacity < copies * e.n_mask_atoms) { + e.copies_mask_capacity = copies * e.n_mask_atoms; + e.copies_mask_atoms = CudaDevicePtr(4 * std::max(e.copies_mask_capacity, 1), CudaAlloc::Synchronous); + } + cuda_err(cudaMemsetAsync(e.fixed, 0, e.rows * sizeof(double2), e.stream)); + for (size_t c = 0; c < copies; c++) { + const std::array &x = c < fixed.size() ? fixed[c] : moving; + Placement p{}; + std::copy(x.begin(), x.begin() + 9, p.r); + std::copy(x.begin() + 9, x.end(), p.t); + std::copy(e.frac, e.frac + 9, p.frac); + place_kernel<<>>(e.rel, e.n_atoms, p, e.geom.nu, e.geom.nv, e.geom.nw, + e.mask_slot, e.mask_radius, e.pos, + e.copies_mask_atoms.get() + 4 * c * e.n_mask_atoms); + cuda_err(cudaGetLastError()); + if (c < fixed.size()) { + e.density->Compute(e.pos, e.grid.get()); + cufft_err(cufftExecR2C(e.Plan(1), e.grid, reinterpret_cast(e.spectrum.get()))); + e.has_fixed = false; + e.Compose(0, e.shifted, nullptr); + add_kernel<<>>(e.shifted, static_cast(e.rows), e.fixed); + cuda_err(cudaGetLastError()); + } + } + float *mask_grid = e.grid.get() + e.GridPoints(); + cufftComplex *mask_spectrum = reinterpret_cast(e.spectrum.get() + e.ComplexPoints()); + e.mask->Compute(reinterpret_cast(e.copies_mask_atoms.get()), + static_cast(copies * e.n_mask_atoms), mask_grid); + cufft_err(cufftExecR2C(e.Plan(1), mask_grid, mask_spectrum)); + fmask_kernel<<>>(reinterpret_cast(mask_spectrum), e.Norm(), + e.row_hkl, static_cast(e.rows), e.geom.nu, e.geom.nv, + e.geom.nw, e.fmask); + cuda_err(cudaGetLastError()); + cuda_err(cudaStreamSynchronize(e.stream)); + e.has_fixed = true; +} + void RigidBodyGPUEngine::Fmask() { RigidBodyGPUEngineImpl &e = *impl_; cuda_err(cudaSetDevice(e.device)); diff --git a/image_analysis/structure_refinement/RigidBodyGPUEngine.h b/image_analysis/structure_refinement/RigidBodyGPUEngine.h index a390fe781..e68ad8e51 100644 --- a/image_analysis/structure_refinement/RigidBodyGPUEngine.h +++ b/image_analysis/structure_refinement/RigidBodyGPUEngine.h @@ -109,6 +109,11 @@ public: // solvent mask of the same placement (ModelMaskGPU) and Fmask. void Fcalc(const double rotation[9], const double translation[3]); void Fmask(); + // Per zone, after SetZone(): fixed atoms beside the body that are copies of it, each placed at (R, t) - + // rotation row-major, then translation, as Fcalc() takes them. Their Fcalc, composed one copy at a time, + // is added to the body's in every composition from now on, and the bulk-solvent mask of them and of the + // body at `moving` is held at the rows in place of Fmask()'s. + void SetFixedCopies(const std::vector> &fixed, const std::array &moving); // The scale at this evaluation's Fcalc and Fmask (ModelScaleGPU): FitModelScale()'s solvent grid, and // the overall scale and anisotropic B at a fixed solvent. FitSolvent() throws ModelScaleGPUTooFewReflections where // the grid's fits would chain (too few reflections for the isotropic start), which the host then does. diff --git a/image_analysis/structure_refinement/RigidBodyRefine.cpp b/image_analysis/structure_refinement/RigidBodyRefine.cpp index 51e3c9a24..bc60a94ca 100644 --- a/image_analysis/structure_refinement/RigidBodyRefine.cpp +++ b/image_analysis/structure_refinement/RigidBodyRefine.cpp @@ -274,6 +274,31 @@ void RigidBodyTarget::SetZone(const gemmi::AsuData> &fo for (const auto &hv : fobs_.v) hkl.push_back(hv.hkl); composition_.emplace(grid, d_min, hkl); + if (fixed != nullptr) + fixed_part_ = RigidBodyFixed(model_, *fixed, cell_, sg_, d_min_, *composition_, nthreads_); +} + +RigidBodyFixedPart RigidBodyFixed(const gemmi::Model &moving, const gemmi::Model &fixed, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup &sg, double d_min, const SymmetryComposition &composition, + size_t nthreads) { + RigidBodyFixedPart out; + gemmi::DensityCalculator dc; + dc.d_min = d_min; + dc.rate = GRID_RATE; + dc.grid.unit_cell = cell; + dc.grid.spacegroup = &sg; + dc.set_refmac_compatible_blur(fixed); + PutModelDensityOnGrid(dc, fixed, {}, nthreads); + composition.Compose(MapToFPhi(dc.grid), dc.blur, out.fcalc, nullptr, nthreads); + + gemmi::Model all = moving; + all.chains.insert(all.chains.end(), fixed.chains.begin(), fixed.chains.end()); + gemmi::Grid mask_grid = RigidBodyZoneGrid(cell, sg, d_min); + PutMaskOnGrid(mask_grid, all, OrbitLeaders(mask_grid, nthreads), nthreads); + const gemmi::FPhiGrid fm = MapToFPhi(mask_grid); + for (const gemmi::Miller &h : composition.Hkl()) + out.fmask.push_back(fm.get_value_by_hkl(h)); + return out; } // The model's Fcalc at the zone's composed indices, from one copy of it on the zone's grid. @@ -300,7 +325,11 @@ bool RigidBodyTarget::Residuals(const double q[6], double *residuals) { const std::vector &hkl = composition_->Hkl(); std::vector> fmask; - if (hold_mask && have_point_) { + if (fixed != nullptr) { + for (size_t m = 0; m < fc.size(); m++) + fc[m] += fixed_part_.fcalc[m]; + fmask = fixed_part_.fmask; + } else if (hold_mask && have_point_) { fmask = fmask_; } else { // The mask is symmetric as gridded, so it is read at h directly. @@ -400,6 +429,9 @@ bool RigidBodyTarget::Jacobian(const double q[6], double *jacobian) { qj[j] += step; Place(qj, column_models_[j]); CopyFcalc(column_models_[j], shifted[j], nullptr); + if (fixed != nullptr) + for (size_t m = 0; m < shifted[j].size(); m++) + shifted[j][m] += fixed_part_.fcalc[m]; })); for (std::future &f : running) f.get(); @@ -498,7 +530,8 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, double d_min, Logger &logger, size_t nthreads, - RigidBodyGPUPool *gpu) { + RigidBodyGPUPool *gpu, + const gemmi::Model *fixed) { const auto t0 = std::chrono::steady_clock::now(); RigidBodyRefineResult result; if (fobs.v.empty()) @@ -513,9 +546,12 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, RigidBodyTargetBase &target = *target_backend; if (!target.Usable()) return result; + target.fixed = fixed; const std::vector ladder = RigidBodyLadder(d_min); - const gemmi::Mat33 gauge = RigidBodyGaugeProjector(sg, cell); + // With fixed atoms beside the body the origin is theirs, not free. + const gemmi::Mat33 gauge = fixed != nullptr ? gemmi::Mat33(0, 0, 0, 0, 0, 0, 0, 0, 0) + : RigidBodyGaugeProjector(sg, cell); double q[6] = {0, 0, 0, 0, 0, 0}; bool any_zone_solved = false; for (double zone : ladder) { diff --git a/image_analysis/structure_refinement/RigidBodyRefine.h b/image_analysis/structure_refinement/RigidBodyRefine.h index d334164a4..91cbd2437 100644 --- a/image_analysis/structure_refinement/RigidBodyRefine.h +++ b/image_analysis/structure_refinement/RigidBodyRefine.h @@ -53,6 +53,11 @@ double RigidBodyReachDeg(const gemmi::Model &model, double d_min); // // `gpu`: evaluate the target on an engine of this pool (RigidBodyGPU.h, CUDA builds only) rather than // on the CPU. The same target either way; null is the CPU. +// +// `fixed`: other atoms of the crystal that do not move - the other copies of a molecular-replacement +// search - whose Fcalc is added to the body's at every evaluation, and which the bulk-solvent mask covers +// as well; that mask is computed once per zone, with the body where the zone starts, and held. With them +// the origin is no longer free, and the translation is refined in all three directions. class RigidBodyGPUPool; RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, const gemmi::UnitCell &cell, @@ -61,7 +66,8 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, double d_min, Logger &logger, size_t nthreads = 1, - RigidBodyGPUPool *gpu = nullptr); + RigidBodyGPUPool *gpu = nullptr, + const gemmi::Model *fixed = nullptr); // --- The pieces RefineRigidBody is built from, exposed for the tests and the GPU backend --- @@ -123,6 +129,16 @@ private: std::vector terms_; }; +// What RefineRigidBody's fixed atoms add to a zone: their Fcalc at the zone's composed indices, and the +// bulk-solvent mask of them and the body together (`moving`, where the zone starts) at the same indices. +struct RigidBodyFixedPart { + std::vector> fcalc; + std::vector> fmask; +}; +RigidBodyFixedPart RigidBodyFixed(const gemmi::Model &moving, const gemmi::Model &fixed, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup &sg, double d_min, const SymmetryComposition &composition, + size_t nthreads); + // The target RefineRigidBody minimises over one resolution zone: the amplitude residuals of the model // placed at q (the angle-axis rotation about the centroid times the rms radius, then the translation, // all in angstroms), with the scale re-fitted at every evaluation. The placement and the counters are @@ -146,6 +162,7 @@ public: virtual bool Jacobian(const double q[6], double *jacobian) = 0; bool hold_mask = false; // for the tests: keep the last evaluation's bulk-solvent mask + const gemmi::Model *fixed = nullptr; // RefineRigidBody's fixed atoms; set before the first zone int evaluations = 0; int jacobians = 0; int unmatched = 0; // zone observations with no calculated amplitude to compare against @@ -186,6 +203,7 @@ private: double f_mean_ = 1; std::vector orbit_leaders_; std::optional composition_; // the zone's + RigidBodyFixedPart fixed_part_; // the zone's, with fixed atoms bool solvent_fitted_ = false; // The last evaluation, which the Jacobian is taken at. From 63e0a854670f1d1079cef71346f22c232e3ab226 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 12 Oct 2026 00:27:17 +0200 Subject: [PATCH 3/3] rugnux --af-model: several copies in the asymmetric unit, tNCS, a search-matched null, faster Copy number. The hard "Vm >= 4: MULTI_COPY, no decision" rule is gone. Copies are searched one after another up to the largest copy number with a Matthews probability of 1% (a Gumbel distribution with Kantardjieff & Rupp's 2003 mode and mean, no resolution dependence); COPIES_LIKELY reports the most probable. Sequential search. Copy n's translation function carries copies 1..n-1 as a known partial structure: the cross term 2 Re(conj(F_known) F(t)) goes on the same FFT grid (CPU FFTW and GPU cuFFT; KnownStructureFactors, AssemblyCorrelation). Its peaks, past those that clash, are ranked by the correlation of the whole assembly with E^2 (the TF itself, linear in |F|^2, rewards copying a placed copy onto an origin shift), the first that stands out (TF Z >= 8), packs and lowers R-free by 0.01 is kept; each copy is refined beside the others held fixed (RefineRigidBody's fixed copies), and the first copy again beside the second, the earlier ones once more at the end. Packing: ClashFraction counts C-alpha clashes with other copies and all their mates, binned (7mzt: 51 s of clash checks -> 1.9 s). Search quality: rotation-function peaks hill-climbed off the 6 deg sampling (6ze4: RF Z 10.4 -> 16.0, 3.3 deg off the deposited orientation -> 0.2); the first copy goes down the TF list until three placements that pack are refined, so a group is CLASH only where no offered orientation packs; groups are searched side by side, and every group the best one leads is searched again with the best one's copy orientations tried first, the better search kept. tNCS: the native Patterson of E^2 (5-10 A) off the lattice; placed copies' orientations join every further copy's search; and where the tNCS vector turns one group's screw axes into another's (TncsRelatesGroups), a group whose search stopped at fewer copies does not lose to the other - UNDECIDED. Measured on an in-house tNCS crystal (vector c/3 along the 3-fold): both hands fit with three copies to within 0.01 in R-free, and where one hand's search failed the other was taken as ENANTIOMORPH; now UNDECIDED. Null: the replicates are placed by the translation function in their own orientation, as the search placed the real one (a null at the real centroid reads high), and the gate is 4.51 = t(8) sqrt(1 + 1/9), 3 sd of a nine-draw null, the gate --model's MODEL_FIT is now read at. Speed: placements scored on the GPU structure-factor engines; further copies' trials score with the first copy's bulk solvent held, the assembly refitted at the end. Survey of 193 sets (jfjoch_mr_replay on stored merges, RTX 5080 shared with other jobs): MR median 1.15 s, max 24 s; demos 0.8-1.1 s. Report: COPIES and COPIES_LIKELY columns, MODEL_PRIOR_COPIES. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- docs/ACKNOWLEDGEMENT.md | 10 +- docs/RUGNUX_REPORT.md | 37 +- .../MolecularReplacement.cpp | 929 ++++++++++++++++-- .../MolecularReplacement.h | 99 +- .../MolecularReplacementGPU.cu | 55 +- .../MolecularReplacementGPU.h | 6 +- rugnux/ResultReport.cpp | 43 +- tests/MolecularReplacementTest.cpp | 121 ++- tools/jfjoch_mr_replay.cpp | 9 +- 9 files changed, 1152 insertions(+), 157 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index bc92f289d..b42dccda8 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -465,7 +465,15 @@ crystallographic software" (2007), J. Appl. Cryst. 40, 658-674 [doi:10.1107/S0021889807021206](https://doi.org/10.1107/S0021889807021206). Whether one copy fits the cell is read off the Matthews coefficient: B. W. Matthews, "Solvent content of protein crystals" (1968), J. Mol. Biol. 33, 491-497 -[doi:10.1016/0022-2836(68)90205-2](https://doi.org/10.1016/0022-2836(68)90205-2). +[doi:10.1016/0022-2836(68)90205-2](https://doi.org/10.1016/0022-2836(68)90205-2). How many copies +to search for is weighed by a distribution of the Matthews coefficient over protein crystals made from +the mode and mean of Kantardjieff and Rupp's: K. A. Kantardjieff and B. Rupp, "Matthews coefficient +probabilities: Improved estimates for unit cell contents of proteins, DNA, and protein-nucleic acid +complex crystals" (2003), Protein Sci. 12, 1865-1871 +[doi:10.1110/ps.0350503](https://doi.org/10.1110/ps.0350503). A copy refined beside the others has the +others' placements read off by least-squares superposition: W. Kabsch, "A solution for the best rotation +to relate two sets of vectors" (1976), Acta Cryst. A32, 922-923 +[doi:10.1107/S0567739476001873](https://doi.org/10.1107/S0567739476001873). ## Software and computing methods diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 96c20d134..9c6ffdb57 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -964,22 +964,35 @@ nothing else. `--af-model ` takes a predicted model - an AlphaFold mmCIF, with the pLDDT in the B-factor column - and places it by molecular replacement: residues below pLDDT 70 removed, B-factors from the -pLDDT, then a rotation function over the point group (6 deg sampling, 4-12 A), a translation function -for the 50 best orientations in every space group the merge cannot tell apart, and the best three -placements per group through the rigid body, bulk solvent and R-free (to `SCORE_D_MIN`, 3 A or the -data's limit). One copy in the asymmetric unit only, for now. +pLDDT, then a rotation function over the point group (6 deg sampling, 4-12 A, each of the 50 best +peaks refined off the sampling), a translation function for those orientations in every space group the +merge cannot tell apart, and the best three placements per group that pack through the rigid body, bulk +solvent and R-free (to `SCORE_D_MIN`, 3 A or the data's limit). Further copies follow one at a time, up +to the largest copy number with a Matthews probability of 1% (`COPIES_LIKELY`: the most probable): the +translation function again with the copies already placed as a known partial structure, its peaks +ranked by the correlation of all the copies with the data, each new copy refined beside the others held +fixed, and kept where its `TF_Z` reaches 8, it packs and R-free drops by 0.01 (`COPIES`, and +`MODEL_PRIOR_COPIES=` for the group chosen). The groups are searched side by side; every group the best +one leads is searched once more with the best one's copies' orientations tried first, and the better +search kept. -The table has one row per group and merge. A row is `PLACED` where its translation-function Z-score -(`TF_Z`) reaches 8, at most 5% of its C-alpha atoms clash with symmetry mates (`CLASH`), one copy fits -the cell (`MATTHEWS`: Vm below 1.5 A^3/Da; `MULTI_COPY`: Vm of 4 or more, which is reported but not -supported yet). The null - the same model refitted from random orientations, whose R-work the -placement's has to beat by 3 sd (`NULL_SIGMA`) - is built for the candidate the choice falls on; one -that does not beat it is marked `NULL` and the choice falls on the next. A placement clashing with its -symmetry mates over twice the limit is not refined. The rotation and translation functions run on the +The table has one row per group and merge. A row is `PLACED` where one of its copies' translation- +function Z-scores (`TF_Z`, the highest) reaches 8, at most 5% of the C-alpha atoms clash with another +copy or a symmetry mate (`CLASH`), and one copy fits the cell (`MATTHEWS`: Vm below 1.5 A^3/Da). The +null - the same model refitted from random orientations, each placed by the translation function as the +search placed the real one (all copies turned together about their centroid where there are several) - +is built for the candidate the choice falls on; its R-work has to beat the null's by 4.51 sd +(`NULL_SIGMA`: 3 sd allowing for the null being nine draws, as `MODEL_FIT_SIGMA`), and one that does not +is marked `NULL` and the choice falls on the next. A placement clashing over twice the limit is not +refined. Translational NCS is looked for in the native Patterson (a peak off the lattice of a fifth of +the origin's): with it, a group whose screw axes the tNCS translation turns into the chosen group's and +whose search placed fewer copies leaves the choice `UNDECIDED`, since one copy beside a translated twin +of itself fits the data no better than at random. The rotation and translation functions, the rigid +body and the structure factors the R-factors are taken on run on the GPU where there is one (`MODEL_PRIOR_SEARCH_DEVICE=`). `MODEL_PRIOR_DECISION=` is then: - `AGREES` - the placed group with the lowest R-free is the answer's, ahead of every other group of - the merge that packs by 0.03 in R-free; + the merge that packs by 0.03 in R-free (with tNCS, see above); - `ENANTIOMORPH` - it is the answer's enantiomorph, and the answer's files take that label (a label only, as with `--model`; no reflection moves). It is taken only by the rule a `--model`'s enantiomorph is taken by: the model has to fit against its null (the choice already required that) diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp index dc2556006..6d1be24f6 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.cpp +++ b/image_analysis/structure_refinement/MolecularReplacement.cpp @@ -9,7 +9,12 @@ #include #include #include +#include #include +#include +#include +#include +#include #include #include @@ -30,6 +35,8 @@ #ifdef JFJOCH_USE_CUDA #include "MolecularReplacementGPU.h" #include "RigidBodyGPU.h" +#include "RigidBodyGPUEngine.h" +#include "ModelStructureFactorsGPU.h" #include "../../common/CUDAWrapper.h" #include "../../common/JFJochException.h" #endif @@ -49,6 +56,8 @@ constexpr double ROTATION_STEP_DEG = 6.0; // 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. @@ -61,9 +70,15 @@ constexpr double SCORE_D_MIN = 3.0; 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. +// Below VM_MIN one copy does not fit the cell. constexpr double VM_MIN = 1.5; -constexpr double VM_TWO_COPIES = 4.0; +// 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; @@ -72,7 +87,10 @@ constexpr double R_FREE_MARGIN = 0.03; constexpr int NULL_REPLICATES = 9; constexpr unsigned NULL_SEED = 20261011; constexpr double NULL_D_MIN = 3.5; -constexpr double NULL_SIGMA_MIN = 3.0; +// 3 sd of the null, one-sided, allowing for its mean and sd being estimated from NULL_REPLICATES draws: +// a new draw lies t(n-1) sqrt(1 + 1/n) sample sds above their mean at that probability, 4.28 x 1.054 for +// n = 9 - the gate model validation's fit is read at (ModelValidation.cpp, MODEL_FIT_SIGMA). +constexpr double NULL_SIGMA_MIN = 4.51; double seconds_since(std::chrono::steady_clock::time_point t0) { return std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); @@ -90,6 +108,16 @@ gemmi::Mat33 random_rotation(std::mt19937 &rng) { 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]; } @@ -142,11 +170,9 @@ gemmi::AsuData> working_amplitudes(const std::vector &merged, double d_min, size_t nthreads) { +// 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> &fmodel, + gemmi::AsuData> &fmask, size_t nthreads) { gemmi::DensityCalculator dc; dc.d_min = d_min; dc.rate = 1.5; @@ -158,17 +184,14 @@ PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGro mask.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up); const std::vector leaders = OrbitLeaders(mask, nthreads); PutMaskOnGrid(mask, st.models[0], leaders, nthreads); - gemmi::AsuData> fmask = MapToFPhi(mask).prepare_asu_data(d_min, 0); + fmask = MapToFPhi(mask).prepare_asu_data(d_min, 0); PutModelDensityOnGrid(dc, st.models[0], leaders, nthreads); - gemmi::AsuData> fmodel = MapToFPhi(dc.grid).prepare_asu_data(d_min, dc.blur, false, false, false); - - gemmi::AsuData> fobs = working_amplitudes(merged, st.cell, sg, d_min); - gemmi::Scaling 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); + 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> &fmodel, const gemmi::Scaling &scaling, + const gemmi::SpaceGroup &sg, const std::vector &merged, double d_min) { const gemmi::GroupOps gops = sg.operations(); const gemmi::ReciprocalAsu asu(&sg); std::unordered_map> obs; @@ -202,6 +225,83 @@ PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGro return out; } +// Model validation's structure-factor engines on the GPU, one per group, copy count and resolution, made +// where the first placement of that kind is scored - or none, where there is no GPU, the cell is too small +// for its gridding or one would take over a quarter of the card. Decided on the card's total memory, so +// the same input always takes the same path. +#ifdef JFJOCH_USE_CUDA +class ScoreEngines { +public: + ModelStructureFactorsGPU *get(const gemmi::Structure &st, const gemmi::SpaceGroup &sg, double d_min) { + if (get_gpu_count() == 0) + return nullptr; + size_t atoms = 0; + for (const gemmi::Chain &ch : st.models[0].chains) + for (const gemmi::Residue &r : ch.residues) + atoms += r.atoms.size(); + std::lock_guard lock(m_); + auto &e = engines_[{sg.xhm(), atoms, d_min}]; + if (!e.tried) { + e.tried = true; + try { + auto sf = std::make_unique(RigidBodyGPUEngine::CurrentDevice(), st.models[0], + st.cell, sg, d_min); + if (sf->Supported() && sf->DeviceBytes() <= sf->DeviceTotalMemory() / 4) { + sf->Reserve(); + e.sf = std::move(sf); + } + } catch (const JFJochException &) { + cuda_clear_error(); + } + } + return e.sf.get(); + } + +private: + struct Entry { + bool tried = false; + std::unique_ptr sf; + }; + std::map, Entry> engines_; + std::mutex m_; +}; +#else +class ScoreEngines {}; +#endif + +// The model's Fcalc and bulk-solvent mask to d_min - on the GPU where `engines` has one - scaled to the +// working amplitudes as model validation scales them (ModelValidation.cpp, fit_model), and the R-factors +// over the reflections to d_min. `solvent`: k_sol and b_sol held at its, not fitted. +PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGroup &sg, + const std::vector &merged, double d_min, size_t nthreads, + ScoreEngines *engines = nullptr, const PlacementScore *solvent = nullptr) { + gemmi::AsuData> fmask, fmodel; +#ifdef JFJOCH_USE_CUDA + if (ModelStructureFactorsGPU *gpu = engines != nullptr ? engines->get(st, sg, d_min) : nullptr) + gpu->Compute(st.models[0], fmodel, fmask); + else +#endif + model_factors(st, d_min, fmodel, fmask, nthreads); + + gemmi::AsuData> fobs = working_amplitudes(merged, st.cell, sg, d_min); + gemmi::Scaling scaling(st.cell, &sg); + scaling.use_solvent = true; + scaling.prepare_points(fmodel, fobs, &fmask); + if (solvent != nullptr) { + // The bulk solvent of an earlier fit of the same crystal, held, as the rigid body holds it. + scaling.k_sol = solvent->k_sol; + scaling.b_sol = solvent->b_sol; + scaling.fix_k_sol = true; + scaling.fix_b_sol = true; + scaling.fit_isotropic_b_approximately(); + scaling.fit_parameters(); + } else { + FitModelScale(scaling, {.b_hi = SolventBMax(st.models[0])}, nthreads); + } + scaling.scale_data(fmodel, &fmask); + return r_factors(fmodel, scaling, sg, merged, d_min); +} + // The rigid-body engines for one model in one group, where there is a GPU. struct RigidBodyEngines { #ifdef JFJOCH_USE_CUDA @@ -236,32 +336,56 @@ RigidBodyEngines make_engines(const gemmi::Structure &st, const gemmi::SpaceGrou // 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 &merged, double d_min, RigidBodyGPUPool *pool, - size_t nthreads) { +// `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 &merged, + double d_min, RigidBodyGPUPool *pool, size_t nthreads, gemmi::Mat33 *moved = nullptr, + const gemmi::Model *fixed = nullptr) { const std::vector 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); + nthreads, pool, fixed); if (!rb.converged) SetModelPositions(st.models[0], start); - return score_placement(st, sg, merged, d_min, nthreads); + if (moved != nullptr) + *moved = rb.converged ? rb.rotation : gemmi::Mat33(); +} + +PlacementScore refine_placement(gemmi::Structure &st, const gemmi::SpaceGroup &sg, + const std::vector &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 fractional position of the model's centroid. +gemmi::Fractional centroid(const gemmi::Structure &st) { + const std::vector 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(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 &merged, double d_min, RigidBodyGPUPool *pool, - size_t nthreads) { + size_t nthreads, ScoreEngines *scorers, + const std::function(const std::vector &)> &translate) { const double null_d_min = std::max(d_min, NULL_D_MIN); std::vector 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 double real = score_placement(placed, sg, coarse, null_d_min, nthreads, scorers).r_work; const std::vector pg = CartesianPointGroup(sg, placed.cell); const double reach = RigidBodyReachDeg(placed.models[0], null_d_min); @@ -272,6 +396,13 @@ double null_sigma(const gemmi::Structure &placed, const gemmi::Mat33 &orientatio if (RotationDistanceDeg(r.multiply(orientation), orientation, pg) >= reach) draws.push_back(r); } + std::vector where; + if (translate) { + std::vector turned; + for (const gemmi::Mat33 &r : draws) + turned.push_back(r.multiply(orientation)); + where = translate(turned); + } std::vector r_work(NULL_REPLICATES); std::atomic next{0}; std::vector> workers; @@ -284,10 +415,11 @@ double null_sigma(const gemmi::Structure &placed, const gemmi::Mat33 &orientatio for (const gemmi::Position &p : pos) c += p; c *= 1.0 / static_cast(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) + c); + 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).r_work; + r_work[i] = refine_placement(rep, sg, coarse, null_d_min, pool, nthreads, scorers).r_work; } })); for (auto &w : workers) @@ -622,7 +754,7 @@ std::vector RotationFunction(const MolecularTransform &ft, const SearchR // 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::SpaceGroup &sg, bool known_part) { const gemmi::GroupOps gops = sg.operations(); const int nops = static_cast(gops.sym_ops.size()); const gemmi::Mat33 frac_t = cell.frac.mat.transpose(); @@ -648,6 +780,8 @@ TranslationSetup MakeTranslationSetup(const SearchReflections &refl, const gemmi 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}; @@ -656,6 +790,10 @@ TranslationSetup MakeTranslationSetup(const SearchReflections &refl, const gemmi 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); @@ -667,15 +805,84 @@ TranslationSetup MakeTranslationSetup(const SearchReflections &refl, const gemmi 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(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> KnownStructureFactors(const MolecularTransform &ft, const TranslationSetup &setup, + const gemmi::UnitCell &cell, + const std::vector &rotations, + const std::vector &positions) { + const size_t terms = setup.Terms(), nops = setup.ops; + std::vector> out(terms, 0.0); + for (size_t c = 0; c < rotations.size(); c++) { + const gemmi::Mat33 rt = rotations[c].transpose(); + const gemmi::Position x = cell.orthogonalize(positions[c]); + for (size_t t = 0; t < terms; t++) + for (size_t j = 0; j < nops; j++) { + const gemmi::Vec3 &sj = setup.s[t * nops + j]; + // (h R_j).t = s_j.x for the Cartesian s_j and x + out[t] += ft.At(rt.multiply(sj)) * setup.phase[t * nops + j] * std::polar(1.0, 2.0 * PI * sj.dot(x)); + } + } + return out; +} + +double AssemblyCorrelation(const MolecularTransform &ft, const TranslationSetup &setup, const gemmi::UnitCell &cell, + const std::vector> &known, const gemmi::Mat33 &rotation, + const gemmi::Fractional &t) { + const size_t terms = setup.Terms(), nops = setup.ops; + const gemmi::Mat33 rt = rotation.transpose(); + const gemmi::Position x = cell.orthogonalize(t); + std::vector ic(terms), bin_sum(setup.bins, 0.0); + std::vector bin_n(setup.bins, 0); + for (size_t i = 0; i < terms; i++) { + std::complex f = known[i]; + for (size_t j = 0; j < nops; j++) { + const gemmi::Vec3 &sj = setup.s[i * nops + j]; + f += ft.At(rt.multiply(sj)) * setup.phase[i * nops + j] * std::polar(1.0, 2.0 * PI * sj.dot(x)); + } + ic[i] = std::norm(f); + bin_sum[setup.bin[i]] += ic[i]; + bin_n[setup.bin[i]]++; + } + double w = 0, sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; + for (size_t i = 0; i < terms; i++) { + const double m = setup.multiplicity[i], a = ic[i] / (bin_sum[setup.bin[i]] / bin_n[setup.bin[i]]), + b = setup.e2[i]; + w += m; + sx += m * a; + sy += m * b; + sxx += m * a * a; + syy += m * b * b; + sxy += m * a * b; + } + const double vx = sxx - sx * sx / w, vy = syy - sy * sy / w; + return vx > 0 && vy > 0 ? (sxy - sx * sy / w) / std::sqrt(vx * vy) : 0.0; +} + TranslationPeak TranslationFunction(const MolecularTransform &ft, const TranslationSetup &setup, - const gemmi::Mat33 &rotation) { + const gemmi::Mat33 &rotation, const std::vector> *known) { 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 + 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> g(terms * nops); std::vector bin_sum(setup.bins, 0.0); @@ -701,6 +908,14 @@ TranslationPeak TranslationFunction(const MolecularTransform &ft, const Translat += 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 v = std::conj((*known)[t]) * g[t * nops + setup.cross_op[e]]; + c[setup.cross_u[e] + static_cast(nu) * (setup.cross_v[e] + static_cast(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 hkl; @@ -729,6 +944,114 @@ TranslationPeak TranslationFunction(const MolecularTransform &ft, const Translat return out; } +// --- translational NCS --- + +TranslationalNcs DetectTranslationalNcs(const std::vector &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> 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 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(std::ceil(3 * cell.a / D_MIN)))); + const int nv = good_even_size(std::max(2 * max_m[1] + 2, static_cast(std::ceil(3 * cell.b / D_MIN)))); + const int nw = good_even_size(std::max(2 * max_m[2] + 2, static_cast(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(nu) * (((k[1] % nv + nv) % nv) + + static_cast(nv) * k[2]); + hkl.data[idx] = std::complex(static_cast(e2), 0.0f); + } + const gemmi::Grid 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(map.nu) * (v + static_cast(map.nv) * w)]; + if (!(value > out.height * origin)) + continue; + const gemmi::Fractional f(static_cast(u) / map.nu, static_cast(v) / map.nv, + static_cast(w) / map.nw); + bool near_lattice = false; + for (const gemmi::Op::Tran &c : gops.cen_ops) { + gemmi::Fractional d(f.x - static_cast(c[0]) / gemmi::Op::DEN, + f.y - static_cast(c[1]) / gemmi::Op::DEN, + f.z - static_cast(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(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, @@ -745,57 +1068,141 @@ gemmi::Structure PlaceSearchModel(const SearchModel &model, const gemmi::Mat33 & return st; } -double ClashFraction(const gemmi::Structure &placed) { +double ClashFraction(const gemmi::Structure &placed, size_t copies) { const gemmi::UnitCell &cell = placed.cell; const gemmi::SpaceGroup *sg = placed.find_spacegroup(); - std::vector ca; - for (const gemmi::Chain &ch : placed.models[0].chains) - for (const gemmi::Residue &r : ch.residues) + const size_t chains = placed.models[0].chains.size(); + copies = std::max(1, std::min(copies, chains)); + // Each copy's C-alpha atoms, fractional, its centroid and its radius. + struct Copy { + std::vector ca; + gemmi::Fractional c0; + double radius = 0; + }; + std::vector 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") - ca.push_back(cell.fractionalize(a.pos)); - if (ca.empty() || sg == nullptr) + 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(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; - gemmi::Fractional c0; - for (const gemmi::Fractional &f : ca) - c0 = c0 + f; - c0 = gemmi::Fractional(c0 * (1.0 / static_cast(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 clash(ca.size(), 0); - for (const gemmi::Op &op : sg->operations().all_ops_sorted()) { - std::vector mate; - for (const gemmi::Fractional &f : ca) { - const std::array x = op.apply_to_xyz({{f.x, f.y, f.z}}); - mate.push_back(gemmi::Fractional(x[0], x[1], x[2])); + // 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(std::floor(x.x / CLASH_A)) + dx + 4096) * 67108864L + + (static_cast(std::floor(x.y / CLASH_A)) + dy + 4096) * 8192L + + (static_cast(std::floor(x.z / CLASH_A)) + dz + 4096); + }; + std::vector images; + for (size_t b = 0; b < copies; b++) + for (const gemmi::Op &op : sg->operations().all_ops_sorted()) { + std::vector mate; + for (const gemmi::Fractional &f : body[b].ca) { + const std::array 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(mate.size()))); + // The lattice shifts that bring this image within reach of some copy. + std::vector 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(b), unmoved}); + } + } } - gemmi::Fractional cm; - for (const gemmi::Fractional &f : mate) - cm = cm + f; - cm = gemmi::Fractional(cm * (1.0 / static_cast(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; + 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(a)) && x.dist_sq(it->x) < CLASH_A * CLASH_A) { + clash = true; break; } } - } + clashing += clash; + } + return static_cast(clashing) / static_cast(total); +} + +gemmi::Structure AssembleCopies(const std::vector &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; +} + +// 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 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 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))); } - return static_cast(std::count(clash.begin(), clash.end(), 1)) / static_cast(ca.size()); + const double sum = std::accumulate(p.begin(), p.end(), 0.0); + for (double &x : p) + x /= sum; + return p; } std::vector SameRotationGroups(const gemmi::SpaceGroup &sg) { @@ -860,7 +1267,12 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path 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 sphere = UniformRotations(ROTATION_STEP_DEG); + ScoreEngines scorers; std::vector 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 @@ -902,17 +1314,58 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path if (static_cast(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 peak_rotations; + std::vector 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 active(peaks.size(), 1); + for (int iter = 0; iter < PEAK_REFINE_ITERATIONS + && std::find(active.begin(), active.end(), 1) != active.end(); iter++) { + std::vector 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 v = gpu ? gpu->RotationFunction(refl, cell, merge.space_group, trial) + : RotationFunction(ft, refl, cell, merge.space_group, trial, nthreads); +#else + const std::vector 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 ? (rf[peaks.front()] - rf_mean) / rf_sd : 0.0); + 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 groups = SameRotationGroups(merge.space_group); - std::vector peak_rotations; - for (size_t p : peaks) - peak_rotations.push_back(rotations[p]); std::vector tf_all; for (const gemmi::SpaceGroup &sg : groups) { const TranslationSetup setup = MakeTranslationSetup(refl, cell, sg); @@ -929,75 +1382,317 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path } result.t_translation += seconds_since(t0); - for (size_t g = 0; g < groups.size(); g++) { + // 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 copy_rotation; // of every copy placed + Logger log = Logger::Buffered(); + }; + const auto search_group = [&](size_t g, const std::vector &seeds) { + GroupResult out; + Logger &glog = out.log; const gemmi::SpaceGroup &sg = groups[g]; MRCandidate cand; cand.merge = static_cast(m); cand.space_group = sg; const gemmi::GroupOps gops = sg.operations(); - cand.vm = cell.volume / static_cast(gops.order()) / model.mass_da; - cand.copies_likely = cand.vm >= VM_TWO_COPIES ? 2 : 1; - const TranslationPeak *tf = &tf_all[g * peaks.size()]; + const double asu_volume = cell.volume / static_cast(gops.order()); + cand.vm = asu_volume / model.mass_da; + const std::vector 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(n) + 1; + if (pm[n] >= MATTHEWS_P_MIN) + cand.copies_max = static_cast(n) + 1; + } + // The orientations searched: the rotation function's peaks and the seeds, which are tried first. + std::vector rots = peak_rotations; + std::vector rfs = peak_rf; + std::vector tfs(tf_all.begin() + g * peaks.size(), tf_all.begin() + (g + 1) * peaks.size()); std::vector 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; }); + 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 tf_seed(seeds.size()); +#ifdef JFJOCH_USE_CUDA + const std::vector 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 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 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 best few through packing, the rigid body and R-free --- - t0 = std::chrono::steady_clock::now(); + // --- the first copy: the best few through packing, the rigid body and R-free --- 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 + // Where not even one copy fits the cell, only the 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(one_copy ? RESCORED_PER_GROUP : 1, static_cast(by_tf.size())); k++) { + const bool fits = cand.vm >= VM_MIN; + // Each copy placed: its orientation and its centroid, where the rigid body left them. + std::vector copies; + std::vector copy_rotation; + std::vector copy_position; + // Down the list until RESCORED_PER_GROUP placements that pack have been refined: a placement + // clashing with its own mates over twice the limit is not refined - the rigid body moves a + // search placement by a few degrees and a fraction of an angstrom, and the clash fraction moves + // by a percent or two with it - and does not count, so that a group is CLASH only where no + // orientation the rotation function offers packs in it. + int refined = 0; + for (size_t k = 0; k < by_tf.size() && refined < (fits ? RESCORED_PER_GROUP : 1); k++) { const size_t i = by_tf[k]; auto placed = std::make_shared( - 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. + PlaceSearchModel(model, rots[i], tf[i].t, cell, sg)); double clash = ClashFraction(*placed); PlacementScore score{NAN, NAN, NAN, NAN}; - if (one_copy && clash <= 2 * CLASH_MAX) { + gemmi::Mat33 moved; + if (fits && clash <= 2 * CLASH_MAX) { if (engines.get() == nullptr) engines = make_engines(*placed, sg, *merge.merged, d_min, 1, logger); - score = refine_placement(*placed, sg, *merge.merged, d_min, engines.get(), nthreads); + score = refine_placement(*placed, sg, *merge.merged, d_min, engines.get(), nthreads, &scorers, + &moved); clash = ClashFraction(*placed); + refined++; + } else if (fits && have) { + continue; + } else if (!fits) { + refined++; } // The refined placement 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.rf_z = rf_sd > 0 ? (rfs[i] - rf_mean) / rf_sd : 0.0; cand.tf_z = tf[i].z; + cand.copy_tf_z = {tf[i].z}; cand.clash_fraction = clash; cand.r_work = score.r_work; cand.r_free = score.r_free; cand.k_sol = score.k_sol; cand.b_sol = score.b_sol; cand.placed = placed; - orientation = rotations[peaks[i]]; + cand.copies = 1; + orientation = rots[i]; + copies = {*placed}; + copy_rotation = {moved.multiply(rots[i])}; + copy_position = {centroid(*placed)}; } - 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, + glog.Info("Model prior: {:<12} RF Z {:5.1f} TF Z {:5.1f} clash {:4.1f}% R-work {:.3f} R-free {:.3f}", + sg.xhm(), rf_sd > 0 ? (rfs[i] - rf_mean) / rf_sd : 0.0, tf[i].z, 100 * clash, score.r_work, score.r_free); } - result.t_rescore += seconds_since(t0); + // --- 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 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 rest; + for (size_t o = 0; o < copies.size(); o++) + if (static_cast(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]); + } + 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(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 rots2 = rots; + rots2.insert(rots2.end(), copy_rotation.begin(), copy_rotation.end()); + if (!setup2) + setup2 = MakeTranslationSetup(refl, cell, sg, true); + const std::vector> known = KnownStructureFactors(ft, *setup2, cell, copy_rotation, + copy_position); + std::vector tf2(rots2.size()); +#ifdef JFJOCH_USE_CUDA + if (gpu) + tf2 = gpu->TranslationFunction(*setup2, rots2, &known); + else +#endif + ParallelFor(static_cast(rots2.size()), nthreads, [&](int i) { + tf2[i] = TranslationFunction(ft, *setup2, rots2[i], &known); + }); + std::vector 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 cc = 0; + PlacementScore score{NAN, NAN, NAN, NAN}; + gemmi::Structure copy, assembly; + gemmi::Mat33 rotation; + }; + std::optional 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 correlation of the assembly with E^2, and the best few of those through the rigid body. + std::vector trials; + for (size_t k = 0; k < by_tf2.size() && trials.size() < COPY_TRIALS; k++) { + Trial trial; + trial.i = by_tf2[k]; + trial.copy = PlaceSearchModel(model, rots2[trial.i], tf2[trial.i].t, cell, sg); + std::vector 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)); + } + ParallelFor(static_cast(trials.size()), nthreads, [&](int k) { + trials[k].cc = AssemblyCorrelation(ft, *setup2, cell, known, rots2[trials[k].i], + tf2[trials[k].i].t); + }); + std::stable_sort(trials.begin(), trials.end(), [](const Trial &a, const Trial &b) { return a.cc > b.cc; }); + if (trials.size() > RESCORED_PER_GROUP) + trials.resize(RESCORED_PER_GROUP); + const gemmi::Structure others = AssembleCopies(copies); + for (Trial &trial : trials) { + std::vector 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]); + glog.Info("Model prior: {:<12} copy {} TF Z {:5.1f} CC {:.3f} clash {:4.1f}% R-work {:.3f} R-free {:.3f}", + sg.xhm(), n, tf2[trial.i].z, trial.cc, 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.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)); + 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.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(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.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)); + 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> running; + for (size_t g = 0; g < groups.size(); g++) + running.push_back(std::async(std::launch::async, search_group, g, std::vector())); + std::vector found; + for (auto &r : running) + found.push_back(r.get()); + // Every group the best one leads 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(g); + if (lead >= 0) { + std::vector> again(found.size()); + for (size_t g = 0; g < found.size(); g++) + if (static_cast(g) != lead && found[g].cand.vm >= VM_MIN) + 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)); } }; @@ -1023,8 +1718,32 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path const auto t_null = std::chrono::steady_clock::now(); RigidBodyEngines engines = make_engines(*c.placed, c.space_group, *merges[merge].merged, result.d_min_score, NULL_REPLICATES, logger); + // One copy: its replicates placed by the translation function, as it was. + std::function(const std::vector &)> translate; + if (c.copies == 1) + translate = [&](const std::vector &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 tf(turned.size()); +#ifdef JFJOCH_USE_CUDA + if (gpu) + tf = gpu->TranslationFunction(setup, turned); + else +#endif + ParallelFor(static_cast(turned.size()), nthreads, [&](int i) { + tf[i] = TranslationFunction(ft, setup, turned[i]); + }); + std::vector out; + for (const TranslationPeak &p : tf) + out.push_back(p.t); + return out; + }; c.null_sigma = null_sigma(*c.placed, orientations[best], c.space_group, *merges[merge].merged, - result.d_min_score, engines.get(), nthreads); + result.d_min_score, engines.get(), nthreads, &scorers, translate); result.t_null += seconds_since(t_null); if (c.null_sigma >= NULL_SIGMA_MIN) break; @@ -1034,6 +1753,14 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path const MRCandidate &c = result.candidates[i]; if (static_cast(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. + if (result.tncs.present && c.copies < result.candidates[best].copies + && TncsRelatesGroups(c.space_group, result.candidates[best].space_group, result.tncs.vector)) { + result.decision = fmt::format("UNDECIDED: {} placed fewer copies than {} beside tNCS", + c.space_group.xhm(), result.candidates[best].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(), diff --git a/image_analysis/structure_refinement/MolecularReplacement.h b/image_analysis/structure_refinement/MolecularReplacement.h index c20a6e44d..30db4f53e 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.h +++ b/image_analysis/structure_refinement/MolecularReplacement.h @@ -3,7 +3,7 @@ #pragma once -// Molecular replacement of ONE copy of a predicted model (rugnux --af-model, experimental): a global +// Molecular replacement of a predicted model (rugnux --af-model, experimental): a global // search for the model's orientation and position in the cell, so that a model with no placement of // its own - an AlphaFold prediction - can be scored against the data, and the space groups the data // cannot tell apart can be scored against it. @@ -13,11 +13,13 @@ // correlation of the observed E^2 with the model's own symmetry-summed |F|^2 - read off the // model's molecular transform, computed once by FFT and interpolated per orientation. // - Translation: the full-symmetry translation function on one FFT per orientation and space group. +// - Further copies, one at a time, as many as the cell can hold: the translation function again with +// the copies already placed as a known partial structure, each new copy refined against the others. // - Rescoring: the best placements through the rigid body, bulk solvent, R-free, a packing check // against the symmetry mates and a null of the same model in random orientations. // -// Everything up to the rescoring runs on the CPU, in both builds; the rigid body runs on the GPU -// where there is one. +// The rotation and translation functions, the rigid body and the structure factors the R-factors are +// taken on run on the GPU where there is one, and on the CPU otherwise; the same functions either way. #include #include @@ -114,14 +116,32 @@ struct TranslationSetup { std::vector s; // terms x ops std::vector> phase; // terms x ops std::vector weight; // terms + std::vector e2, multiplicity; // terms std::vector bin; // terms int bins = 0; std::vector pair_j, pair_k; // the ordered pairs std::vector target_u, target_v, target_l; // terms x pairs + // With a known partial structure (MakeTranslationSetup's `known_part`): where the cross terms + // conj(F_known) g_j go, at h R_j - or, as their conjugate F_known conj(g_j), at -h R_j where that + // has l >= 0; both where l = 0. One entry per coefficient: its term, operator, index and whether it + // is the conjugate. + std::vector cross_term, cross_op, cross_u, cross_v, cross_l; + std::vector cross_conj; size_t Terms() const { return weight.size(); } }; +// `known_part`: the grid is sized for the cross terms with a known partial structure as well, which +// reach out to h R_j where the copies' own terms reach h R_j - h R_k (in a polar group no further than +// the plane the polar axis is free in). TranslationSetup MakeTranslationSetup(const SearchReflections &refl, const gemmi::UnitCell &cell, - const gemmi::SpaceGroup &sg); + const gemmi::SpaceGroup &sg, bool known_part = false); + +// F of copies of the model already placed - copy c turned by rotations[c] about its centroid, the +// centroid at fractional positions[c] - at each term of `setup`, summed over the copies and their +// symmetry mates, in the units the translation function's g_j are in. +std::vector> KnownStructureFactors(const MolecularTransform &ft, const TranslationSetup &setup, + const gemmi::UnitCell &cell, + const std::vector &rotations, + const std::vector &positions); // The translation function of the model in orientation `rotation`: its highest point as a Z-score over // every point of the function, and where it is (the fractional position of the model's centroid). @@ -129,20 +149,60 @@ struct TranslationPeak { double z = 0; gemmi::Fractional t; }; +// With `known` (KnownStructureFactors on a setup made with known_part), the function of one more copy +// beside the known ones: the overlap of E^2 - 1 with |F_known + F(t)|^2, whose cross term +// 2 Re(conj(F_known) F(t)) is a Fourier series in t too and goes on the same grid. TranslationPeak TranslationFunction(const MolecularTransform &ft, const TranslationSetup &setup, - const gemmi::Mat33 &rotation); + const gemmi::Mat33 &rotation, + const std::vector> *known = nullptr); + +// The correlation of E^2 with |F|^2 of the known copies (KnownStructureFactors) and one more copy turned by +// `rotation` with its centroid at `t`, over the terms of `setup`, each weighted by its multiplicity and each +// bin's |F|^2 over its mean: the score the translation function's peaks for a further copy are ranked by. +// Unlike the translation function, which is linear in |F|^2, it does not reward a copy that only repeats +// the known ones - their own copy moved by one of the group's origin shifts doubles F on half the +// reflections and cancels it on the other half, which the overlap with E^2 counts as a gain. +double AssemblyCorrelation(const MolecularTransform &ft, const TranslationSetup &setup, const gemmi::UnitCell &cell, + const std::vector> &known, const gemmi::Mat33 &rotation, + const gemmi::Fractional &t); // The model placed: rotated about its centroid by `rotation`, the centroid at `t`, in `cell` and `sg`. gemmi::Structure PlaceSearchModel(const SearchModel &model, const gemmi::Mat33 &rotation, const gemmi::Fractional &t, const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg); -// The fraction of the model's C-alpha atoms within 3 A of a C-alpha of a symmetry mate. -double ClashFraction(const gemmi::Structure &placed); +// The fraction of the C-alpha atoms within 3 A of a C-alpha of another copy: a symmetry mate, or one of +// the other `copies` the structure's chains are split into evenly (in order), or their mates. +double ClashFraction(const gemmi::Structure &placed, size_t copies = 1); + +// Copies of the model as one structure: their chains in order, renamed so that no two share a name. +gemmi::Structure AssembleCopies(const std::vector &copies); + +// The probabilities of 1, 2, ... copies of a model of `mass_da` in the asymmetric unit of a cell of +// `asu_volume` A^3, from the distribution of the Matthews coefficient over protein crystals; one entry +// per copy number whose Vm is at least 1.5 A^3/Da. +std::vector MatthewsProbabilities(double asu_volume, double mass_da); // The space groups with the rotations and the centring of `sg`, in its setting: those merged // intensities cannot tell apart, apart from systematic absences. `sg` itself first. std::vector SameRotationGroups(const gemmi::SpaceGroup &sg); +// Translational non-crystallographic symmetry: copies in the same orientation, a translation apart, make +// the native Patterson peak off its origin where no lattice vector is. The highest peak of the Patterson +// of E^2 (5-10 A), at least 15 A from the origin and from every centring vector, as a fraction of the +// origin's height, and where it is; tNCS where that is a fifth or more. +struct TranslationalNcs { + bool present = false; + double height = 0; // over the origin's + gemmi::Fractional vector; +}; +TranslationalNcs DetectTranslationalNcs(const std::vector &merged, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup &sg); + +// Whether groups `a` and `b` (the same rotations) differ only in screw translations that the tNCS vector +// `v` makes up: then copies a translation v apart put one group's operators onto the other's, and neither +// is told from the other by a search that has not placed all of them. +bool TncsRelatesGroups(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b, const gemmi::Fractional &v); + // One merge the model is searched against: the answer's, or the twin hypothesis' in a subgroup. struct MRMerge { std::string name; // "MERGE" or "TWIN_SUBGROUP_MERGE" @@ -151,20 +211,23 @@ struct MRMerge { gemmi::SpaceGroup space_group; // the group it was merged in }; -// One space group's best single-copy placement. +// One space group's best placement: its first copy and every further one the search added. struct MRCandidate { int merge = 0; // index into the merges searched gemmi::SpaceGroup space_group; double vm = 0; // Matthews coefficient of one copy - int copies_likely = 1; - double rf_z = 0; // of the orientation, over the rotation function - double tf_z = 0; - double clash_fraction = 0; + int copies_likely = 1; // the most probable copy number, by Matthews + int copies_max = 1; // the most the search tried (Matthews probability >= 1%) + int copies = 0; // placed + double rf_z = 0; // of the first copy's orientation, over the rotation function + double tf_z = 0; // the highest of the copies' + std::vector copy_tf_z; // each copy's, the later ones beside the earlier + double clash_fraction = 0; // of all the copies double r_work = NAN, r_free = NAN; // after the rigid body, to d_min_score double k_sol = NAN, b_sol = NAN; // the bulk solvent R was computed with double null_sigma = NAN; // R-work above its random-orientation null, in sd std::string verdict; // PLACED, or why not - std::shared_ptr placed; + std::shared_ptr placed; // every copy }; struct MolecularReplacementResult { @@ -180,6 +243,7 @@ struct MolecularReplacementResult { int chosen = -1; std::string decision; bool gpu_search = false; // the rotation and translation functions ran on the GPU + TranslationalNcs tncs; // of the answer's merge // An ENANTIOMORPH decision is taken into the answer's label only as model validation takes a model's // enantiomorph (ValidateAgainstModel): the placed model has to fit against its null and the anomalous // map at its atoms must not be inverted. Filled by the caller, which runs that validation. @@ -199,10 +263,11 @@ struct MolecularReplacementResult { double t_model = 0, t_transform = 0, t_rotation = 0, t_translation = 0, t_rescore = 0, t_null = 0; }; -// Search `merges` (the answer's first) for one copy of the predicted model at `model_path`, in every -// space group with the rotations of each merge's group, and decide. A candidate is PLACED where its -// translation function stands out (TF Z >= 8), it packs, one copy fits the cell, and its R-work beats -// its own null. The model chooses among the answer merge's groups where one placed candidate leads +// Search `merges` (the answer's first) for the predicted model at `model_path`, in every space group with +// the rotations of each merge's group - one copy, then further copies one at a time while each new one +// stands out, packs and lowers R-free - and decide. A candidate is PLACED where a copy's translation +// function stands out (TF Z >= 8), the copies pack, one copy fits the cell, and its R-work beats its own +// null. The model chooses among the answer merge's groups where one placed candidate leads // every other by an R-free margin; it chooses the twin hypothesis' merge only where none of the // answer's groups can hold the model at all (Matthews, packing) and a group of the subgroup can. MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector &merges, diff --git a/image_analysis/structure_refinement/MolecularReplacementGPU.cu b/image_analysis/structure_refinement/MolecularReplacementGPU.cu index 89cf87171..8a062fee8 100644 --- a/image_analysis/structure_refinement/MolecularReplacementGPU.cu +++ b/image_analysis/structure_refinement/MolecularReplacementGPU.cu @@ -194,6 +194,28 @@ __global__ void translation_gather(const int *offset, const int *entry, const fl coef[c] = sum; } +// The cross terms with a known partial structure, added to the coefficients: w conj(F_known) g_j, or its +// conjugate, for every entry landing on a grid point, in the order they were listed. +__global__ void translation_gather_cross(const int *offset, const int *entry, const float2 *copies, + const float *scaled, const float2 *known, const int *cross_term, + const int *cross_op, const char *cross_conj, int nops, size_t cells, + float2 *coef) { + const size_t c = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (c >= cells) + return; + float2 sum = coef[c]; + for (int e = offset[c]; e < offset[c + 1]; e++) { + const int x = entry[e], t = cross_term[x]; + const float2 k = known[t]; + const float2 g = copies[static_cast(t) * nops + cross_op[x]]; + const float w = scaled[t]; + const float re = k.x * g.x + k.y * g.y, im = k.x * g.y - k.y * g.x; // conj(k) g + sum.x += w * re; + sum.y += w * (cross_conj[x] ? -im : im); + } + coef[c] = sum; +} + // Per block: the sum and the sum of squares of the function, and its highest point (the lowest index // on a tie). __global__ void translation_stats(const float *map, size_t n, double *sum, double *sum2, float *best, @@ -329,12 +351,13 @@ std::vector MolecularReplacementGPU::RotationFunction(const SearchReflec return out; } -std::vector MolecularReplacementGPU::TranslationFunction(const TranslationSetup &setup, - const std::vector &rotations) { +std::vector MolecularReplacementGPU::TranslationFunction( + const TranslationSetup &setup, const std::vector &rotations, + const std::vector> *known) { std::vector out(rotations.size()); const int terms = static_cast(setup.Terms()), nops = setup.ops; const int npairs = static_cast(setup.pair_j.size()); - if (nops < 2 || terms == 0 || rotations.empty()) + if (terms == 0 || rotations.empty() || (nops < 2 && known == nullptr)) return out; cudaStream_t stream = state_->stream; const int nu = setup.nu, nv = setup.nv, nw = setup.nw, nwh = nw / 2 + 1; @@ -356,6 +379,28 @@ std::vector MolecularReplacementGPU::TranslationFunction(const entry[fill[(static_cast(setup.target_u[e]) * nv + setup.target_v[e]) * nwh + setup.target_l[e]]++] = static_cast(e); } + // The same for the cross terms. + std::vector cross_offset(cells + 1, 0), cross_entry; + std::vector known_f; + if (known != nullptr) { + for (size_t e = 0; e < setup.cross_term.size(); e++) + cross_offset[(static_cast(setup.cross_u[e]) * nv + setup.cross_v[e]) * nwh + setup.cross_l[e] + 1]++; + for (size_t c = 0; c < cells; c++) + cross_offset[c + 1] += cross_offset[c]; + cross_entry.resize(cross_offset[cells]); + std::vector fill(cross_offset.begin(), cross_offset.end() - 1); + for (size_t e = 0; e < setup.cross_term.size(); e++) + cross_entry[fill[(static_cast(setup.cross_u[e]) * nv + setup.cross_v[e]) * nwh + + setup.cross_l[e]]++] = static_cast(e); + for (const std::complex &k : *known) + known_f.push_back(make_float2(static_cast(k.real()), static_cast(k.imag()))); + } + auto d_cross_offset = upload(cross_offset, stream); + auto d_cross_entry = upload(cross_entry, stream); + auto d_known = upload(known_f, stream); + auto d_cross_term = upload(setup.cross_term, stream); + auto d_cross_op = upload(setup.cross_op, stream); + auto d_cross_conj = upload(setup.cross_conj, stream); std::vector phase(setup.phase.size()); for (size_t i = 0; i < phase.size(); i++) phase[i] = make_float2(static_cast(setup.phase[i].real()), static_cast(setup.phase[i].imag())); @@ -391,6 +436,10 @@ std::vector MolecularReplacementGPU::TranslationFunction(const scaled); translation_gather<<((cells + THREADS - 1) / THREADS), THREADS, 0, stream>>>( d_offset, d_entry, copies, scaled, d_pair_j, d_pair_k, nops, npairs, cells, coef); + if (known != nullptr) + translation_gather_cross<<((cells + THREADS - 1) / THREADS), THREADS, 0, stream>>>( + d_cross_offset, d_cross_entry, copies, scaled, d_known, d_cross_term, d_cross_op, d_cross_conj, + nops, cells, coef); cufft_err(cufftExecC2R(plan, reinterpret_cast(coef.get()), map.get())); translation_stats<<>>( map, points, d_sum.get() + r * stat_blocks, d_sum2.get() + r * stat_blocks, diff --git a/image_analysis/structure_refinement/MolecularReplacementGPU.h b/image_analysis/structure_refinement/MolecularReplacementGPU.h index 38f0ebdbf..4efb7f170 100644 --- a/image_analysis/structure_refinement/MolecularReplacementGPU.h +++ b/image_analysis/structure_refinement/MolecularReplacementGPU.h @@ -28,9 +28,11 @@ public: std::vector RotationFunction(const SearchReflections &refl, const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, const std::vector &rotations); - // The translation function of one group at each of `rotations`. + // The translation function of one group at each of `rotations`; with `known`, beside a known partial + // structure (KnownStructureFactors). std::vector TranslationFunction(const TranslationSetup &setup, - const std::vector &rotations); + const std::vector &rotations, + const std::vector> *known = nullptr); private: std::unique_ptr state_; diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 353a46c35..645127f20 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -2579,7 +2579,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyReal("SCORE_D_MIN", mp.d_min_score, "{:.2f}")); Add(s, KeyText("MODEL_PRIOR_DECISION", mp.decision)); if (mp.chosen >= 0) + { Add(s, KeyText("MODEL_PRIOR_SPACE_GROUP", mp.candidates[mp.chosen].space_group.xhm())); + Add(s, KeyInt("MODEL_PRIOR_COPIES", mp.candidates[mp.chosen].copies)); + } if (mp.hand_check.performed) { const auto &h = mp.hand_check; Add(s, KeyEnum("MODEL_PRIOR_HAND_CHECK", h.adopted ? "ADOPTED" : "VETOED", {"ADOPTED", "VETOED"})); @@ -2601,22 +2604,25 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, Blank()); ReportEntry e; e.kind = ReportEntry::Kind::Table; - e.table.columns = {"MERGE", "SPACE_GROUP", "VM", "RF_Z", "TF_Z", "CLASH", "R_WORK", "R_FREE", - "NULL_SIGMA", "VERDICT"}; - e.table.text_header = " MERGE SPACE_GROUP VM RF_Z TF_Z CLASH R_WORK R_FREE NULL_SIGMA VERDICT\n" - " ------------------- ----------- ----- ----- ----- ----- ------ ------ ---------- ----------"; + e.table.columns = {"MERGE", "SPACE_GROUP", "VM", "COPIES", "COPIES_LIKELY", "RF_Z", "TF_Z", "CLASH", + "R_WORK", "R_FREE", "NULL_SIGMA", "VERDICT"}; + e.table.text_header = " MERGE SPACE_GROUP VM COPIES LIKELY RF_Z TF_Z CLASH R_WORK R_FREE NULL_SIGMA VERDICT\n" + " ------------------- ----------- ----- ------ ------ ----- ----- ----- ------ ------ ---------- ----------"; for (const auto &c : mp.candidates) { const std::string merge = c.merge == 0 ? "MERGE" : "TWIN_SUBGROUP_MERGE"; const auto r = [](double v) { return std::isfinite(v) ? fmt::format("{:.3f}", v) : std::string("-"); }; e.table.text_rows.push_back(fmt::format( - " {:<19s} {:<11s} {:5.2f} {:5.1f} {:5.1f} {:4.1f}% {:>6s} {:>6s} {:>10s} {}", - merge, c.space_group.xhm(), c.vm, c.rf_z, c.tf_z, 100 * c.clash_fraction, r(c.r_work), + " {:<19s} {:<11s} {:5.2f} {:6d} {:6d} {:5.1f} {:5.1f} {:4.1f}% {:>6s} {:>6s} {:>10s} {}", + merge, c.space_group.xhm(), c.vm, c.copies, c.copies_likely, c.rf_z, c.tf_z, + 100 * c.clash_fraction, r(c.r_work), r(c.r_free), std::isfinite(c.null_sigma) ? fmt::format("{:+.1f}", c.null_sigma) : std::string("-"), c.verdict)); std::vector row; row.push_back(KeyText("", merge).value); row.push_back(KeyText("", c.space_group.xhm()).value); row.push_back(KeyReal("", c.vm, "{:.2f}").value); + row.push_back(KeyInt("", c.copies).value); + row.push_back(KeyInt("", c.copies_likely).value); row.push_back(KeyReal("", c.rf_z, "{:.1f}").value); row.push_back(KeyReal("", c.tf_z, "{:.1f}").value); row.push_back(KeyReal("", c.clash_fraction, "{:.3f}").value); @@ -2629,17 +2635,20 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, std::move(e)); Add(s, Blank()); Add(s, Prose(fmt::format( - " One copy of the predicted model, trimmed to pLDDT >= 70, was searched for in every space group\n" - " each merge cannot tell apart (rotation function on a 6 deg sampling, translation function,\n" - " {:.1f}-{:.1f} A), and the best placements taken through the rigid body; R_WORK and R_FREE are\n" - " to SCORE_D_MIN. PLACED needs TF_Z >= 8,\n" - " at most 5% of the C-alpha atoms clashing with symmetry mates, one copy fitting the cell\n" - " (MATTHEWS: Vm below 1.5; MULTI_COPY: Vm of 4 or more, not supported yet) and, for the\n" - " candidate the choice falls on, an R-work at least 3 sd better than the same model refitted\n" - " from random orientations (NULL where it is not; the next candidate is then taken). The model chooses\n" - " within a merge where its placed group leads every other that packs by 0.03 in R-free, and\n" - " the twin hypothesis' merge - searched only then - where no group of the answer's merge can\n" - " hold it at all. A group one copy cannot fill is checked for packing only, not refined.\n" + " The predicted model, trimmed to pLDDT >= 70, was searched for in every space group each merge\n" + " cannot tell apart (rotation function on a 6 deg sampling, translation function, {:.1f}-{:.1f} A),\n" + " and the best placements taken through the rigid body; then further copies, one at a time, up to\n" + " the largest copy number with a Matthews probability of 1% (COPIES_LIKELY: the most probable),\n" + " each with the translation function beside the copies already placed, kept where its TF_Z is 8\n" + " or more, it packs and R-free drops by 0.01. R_WORK and R_FREE are of all COPIES, to\n" + " SCORE_D_MIN; TF_Z is the highest of the copies'. PLACED needs TF_Z >= 8, at most 5% of the\n" + " C-alpha atoms clashing with another copy or a symmetry mate, one copy fitting the cell\n" + " (MATTHEWS: Vm below 1.5) and, for the candidate the choice falls on, an R-work at least 3 sd\n" + " better than the same copies refitted from random orientations (NULL where it is not; the next\n" + " candidate is then taken). The model chooses within a merge where its placed group leads every\n" + " other that packs by 0.03 in R-free, and the twin hypothesis' merge - searched only then - where\n" + " no group of the answer's merge can hold it at all. A group one copy cannot fit is checked for\n" + " packing only, not refined.\n" " The answer's files take the enantiomorph from it - only where model validation of the placed\n" " model accepts it as it accepts a --model's: a fit over its null and an anomalous map at its\n" " atoms that is not inverted (ENANTIOMORPH_VETOED otherwise) - and nothing else; another group it\n" diff --git a/tests/MolecularReplacementTest.cpp b/tests/MolecularReplacementTest.cpp index be04c8f21..cc401a164 100644 --- a/tests/MolecularReplacementTest.cpp +++ b/tests/MolecularReplacementTest.cpp @@ -11,9 +11,12 @@ #include "../image_analysis/structure_refinement/ModelFFT.h" #include "../image_analysis/structure_refinement/MolecularReplacement.h" +#include "../image_analysis/structure_refinement/RigidBodyRefine.h" +#include "../common/Logger.h" #ifdef JFJOCH_USE_CUDA #include "../common/CUDAWrapper.h" #include "../image_analysis/structure_refinement/MolecularReplacementGPU.h" +#include "../image_analysis/structure_refinement/RigidBodyGPU.h" #endif namespace { @@ -39,7 +42,7 @@ namespace { c *= 1.0 / pos.size(); for (const auto &p : pos) { gemmi::Atom a; - a.name = "C"; + a.name = "CA"; // every atom counts in the clash check a.element = gemmi::El::C; a.pos = gemmi::Position(p - c); a.b_iso = 15.0f; @@ -173,6 +176,114 @@ TEST_CASE("MolecularReplacement_SyntheticPlacement", "[MolecularReplacement]") { } } +// Two copies in the asymmetric unit: with the first one known, the translation function of the second +// at its true orientation peaks at its true position - not merely at one of the group's origins, which +// only the cross term with the known copy tells apart. +TEST_CASE("MolecularReplacement_SecondCopyBesideKnownFirst", "[MolecularReplacement]") { + const SearchModel model = synthetic_model(); + const gemmi::UnitCell cell(52, 60, 70, 90, 90, 90); + const gemmi::SpaceGroup &sg = gemmi::get_spacegroup_by_name("P 21 21 21"); + const gemmi::Mat33 r1 = axis_angle({1, 2, 3}, 50), r2 = axis_angle({-2, 1, 1}, 110); + const gemmi::Fractional t1(0.17, 0.31, 0.09), t2(0.41, 0.12, 0.33); + const gemmi::Structure crystal = AssembleCopies({PlaceSearchModel(model, r1, t1, cell, sg), + PlaceSearchModel(model, r2, t2, cell, sg)}); + const auto merged = observed(crystal); + const SearchReflections refl = MakeSearchReflections(merged, cell, sg, 4.0, 12.0); + const MolecularTransform ft(model, 4.0, 4); + const TranslationSetup setup = MakeTranslationSetup(refl, cell, sg, true); + const auto known = KnownStructureFactors(ft, setup, cell, {r1}, {t1}); + const TranslationPeak peak = TranslationFunction(ft, setup, r2, &known); + CHECK(peak.z > 10.0); + for (double d : {peak.t.x - t2.x, peak.t.y - t2.y, peak.t.z - t2.z}) + CHECK(std::fabs(d - std::round(d)) < 0.04); + // Two copies of the same body in one cell do not clash where they are placed apart. + CHECK(ClashFraction(crystal, 2) < 0.05); + CHECK(ClashFraction(AssembleCopies({PlaceSearchModel(model, r1, t1, cell, sg), + PlaceSearchModel(model, r2, t1, cell, sg)}), 2) > 0.5); +} + +TEST_CASE("MolecularReplacement_MatthewsProbabilities", "[MolecularReplacement]") { + // Vm 2.5 A^3/Da for one copy: one copy is by far the most probable, four impossible (Vm 0.6). + const auto p = MatthewsProbabilities(2.5 * 20000, 20000); + REQUIRE(p.size() == 1); + CHECK(p[0] == Catch::Approx(1.0)); + // Vm 9 for one copy: three or four copies (Vm 3, 2.25) are the likely ones, one copy is not. + const auto q = MatthewsProbabilities(9.0 * 20000, 20000); + REQUIRE(q.size() == 6); + CHECK(q[0] < 0.01); + CHECK(std::max_element(q.begin(), q.end()) - q.begin() == 3); +} + +TEST_CASE("MolecularReplacement_TranslationalNcs", "[MolecularReplacement]") { + const SearchModel model = synthetic_model(); + const gemmi::UnitCell cell(52, 60, 90, 90, 90, 90); + const gemmi::SpaceGroup &sg = gemmi::get_spacegroup_by_name("P 21 21 21"); + const gemmi::Mat33 r = axis_angle({1, 2, 3}, 50); + const gemmi::Fractional t(0.17, 0.31, 0.09), v(0.0, 0.0, 0.5); + const auto one = observed(PlaceSearchModel(model, r, t, cell, sg)); + CHECK_FALSE(DetectTranslationalNcs(one, cell, sg).present); + const auto two = observed(AssembleCopies({PlaceSearchModel(model, r, t, cell, sg), + PlaceSearchModel(model, r, gemmi::Fractional(t.x, t.y, t.z + v.z), cell, sg)})); + const TranslationalNcs tncs = DetectTranslationalNcs(two, cell, sg); + REQUIRE(tncs.present); + CHECK(std::fabs(tncs.vector.x - std::round(tncs.vector.x)) < 0.03); + CHECK(std::fabs(tncs.vector.y - std::round(tncs.vector.y)) < 0.03); + CHECK(std::fabs(tncs.vector.z - 0.5) < 0.03); +} + +TEST_CASE("MolecularReplacement_TncsRelatesGroups", "[MolecularReplacement]") { + const auto &p3121 = gemmi::get_spacegroup_by_name("P 31 2 1"), &p3221 = gemmi::get_spacegroup_by_name("P 32 2 1"); + // A translation of c/3 turns a 3(1) axis into a 3(2) one, wherever it lies in the plane. + CHECK(TncsRelatesGroups(p3121, p3221, gemmi::Fractional(0.3, 0.1, 1.0 / 3))); + CHECK_FALSE(TncsRelatesGroups(p3121, p3221, gemmi::Fractional(0.3, 0.1, 0.1))); + // A translation of a/2 turns 2(1) along a into 2, and nothing else. + const auto &p212121 = gemmi::get_spacegroup_by_name("P 21 21 21"); + CHECK(TncsRelatesGroups(gemmi::get_spacegroup_by_name("P 2 21 21"), p212121, gemmi::Fractional(0.5, 0, 0))); + CHECK_FALSE(TncsRelatesGroups(gemmi::get_spacegroup_by_name("P 21 2 2"), p212121, gemmi::Fractional(0.5, 0, 0))); +} + +// A copy refined beside another one held fixed goes back where it belongs, on the CPU and on the GPU alike. +TEST_CASE("MolecularReplacement_RigidBodyBesideFixedCopy", "[MolecularReplacement]") { + const SearchModel model = synthetic_model(); + const gemmi::UnitCell cell(52, 60, 70, 90, 90, 90); + const gemmi::SpaceGroup &sg = gemmi::get_spacegroup_by_name("P 21 21 21"); + const gemmi::Mat33 r1 = axis_angle({1, 2, 3}, 50), r2 = axis_angle({-2, 1, 1}, 110); + const gemmi::Fractional t1(0.17, 0.31, 0.09), t2(0.41, 0.12, 0.33); + const gemmi::Structure first = PlaceSearchModel(model, r1, t1, cell, sg); + const gemmi::Structure truth = PlaceSearchModel(model, r2, t2, cell, sg); + const auto merged = observed(AssembleCopies({first, truth})); + gemmi::AsuData> fobs; + fobs.unit_cell_ = cell; + fobs.spacegroup_ = &sg; + for (const MergedReflection &m : merged) + fobs.v.push_back({{{m.h, m.k, m.l}}, {m.F, 1.0f}}); + fobs.ensure_asu(); + fobs.ensure_sorted(); + // 3 deg and 0.5 A out. + const gemmi::Structure start = PlaceSearchModel(model, axis_angle({0, 0, 1}, 3).multiply(r2), + gemmi::Fractional(t2.x + 0.5 / 52, t2.y, t2.z), cell, sg); + const auto worst = [&](const gemmi::Structure &st) { + const auto a = ModelPositions(st.models[0]), b = ModelPositions(truth.models[0]); + double w = 0; + for (size_t i = 0; i < a.size(); i++) + w = std::max(w, a[i].dist(b[i])); + return w; + }; + Logger logger("MolecularReplacement_RigidBodyBesideFixedCopy"); + gemmi::Structure cpu = start; + RefineRigidBody(cpu.models[0], cell, sg, fobs, 3.0, logger, 4, nullptr, &first.models[0]); + CHECK(worst(start) > 1.0); + CHECK(worst(cpu) < 0.1); +#ifdef JFJOCH_USE_CUDA + if (get_gpu_count() > 0) { + gemmi::Structure gpu = start; + auto pool = RigidBodyGPUPool::Create(gpu.models[0], cell, sg, 3.0, fobs.v.size(), 1, logger); + REQUIRE(pool); + RefineRigidBody(gpu.models[0], cell, sg, fobs, 3.0, logger, 4, pool.get(), &first.models[0]); + CHECK(worst(gpu) < 0.1); + } +#endif +} #ifdef JFJOCH_USE_CUDA TEST_CASE("MolecularReplacement_GPUMatchesCPU", "[MolecularReplacement]") { @@ -207,5 +318,13 @@ TEST_CASE("MolecularReplacement_GPUMatchesCPU", "[MolecularReplacement]") { CHECK(tf_gpu[i].t.z == Catch::Approx(cpu.t.z)); } } + // Beside a known copy. + const TranslationSetup setup2 = MakeTranslationSetup(refl, cell, sg, true); + const auto known = KnownStructureFactors(ft, setup2, cell, {rotations[5]}, {gemmi::Fractional(0.3, 0.1, 0.4)}); + const auto tf2_gpu = gpu.TranslationFunction(setup2, orientations, &known); + for (size_t i = 0; i < orientations.size(); i++) { + const TranslationPeak cpu = TranslationFunction(ft, setup2, orientations[i], &known); + CHECK(std::fabs(tf2_gpu[i].z - cpu.z) < 1e-3); + } } #endif diff --git a/tools/jfjoch_mr_replay.cpp b/tools/jfjoch_mr_replay.cpp index 2d751dcea..addf45aea 100644 --- a/tools/jfjoch_mr_replay.cpp +++ b/tools/jfjoch_mr_replay.cpp @@ -90,8 +90,10 @@ int main(int argc, char **argv) { return 0; } std::cout << "MODEL_PRIOR_DECISION= " << mr.decision << "\n"; - if (mr.chosen >= 0) + if (mr.chosen >= 0) { std::cout << "MODEL_PRIOR_SPACE_GROUP= " << mr.candidates[mr.chosen].space_group.xhm() << "\n"; + std::cout << "MODEL_PRIOR_COPIES= " << mr.candidates[mr.chosen].copies << "\n"; + } if (mr.hand_check.performed) std::cout << "MODEL_PRIOR_HAND_CHECK= " << (mr.hand_check.adopted ? "ADOPTED" : "VETOED") << "\n"; std::cout << fmt::format("MODEL_PRIOR_SECONDS= {:.2f}\n", mr.t_model + mr.t_transform + mr.t_rotation @@ -100,8 +102,9 @@ int main(int argc, char **argv) { for (const MRCandidate &c : mr.candidates) { const auto r = [](double v) { return std::isfinite(v) ? fmt::format("{:.3f}", v) : std::string("-"); }; std::cout << fmt::format( - "# {:<19s} {:<11s} {:5.2f} {:5.1f} {:5.1f} {:4.1f}% {:>6s} {:>6s} {:>10s} {}\n", - c.merge == 0 ? "MERGE" : "TWIN_SUBGROUP_MERGE", c.space_group.xhm(), c.vm, c.rf_z, c.tf_z, + "# {:<19s} {:<11s} {:5.2f} {:6d} {:6d} {:5.1f} {:5.1f} {:4.1f}% {:>6s} {:>6s} {:>10s} {}\n", + c.merge == 0 ? "MERGE" : "TWIN_SUBGROUP_MERGE", c.space_group.xhm(), c.vm, c.copies, c.copies_likely, + c.rf_z, c.tf_z, 100 * c.clash_fraction, r(c.r_work), r(c.r_free), std::isfinite(c.null_sigma) ? fmt::format("{:+.1f}", c.null_sigma) : std::string("-"), c.verdict); }