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 326ac6e9e..eff93826a 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -974,22 +974,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 72ed74c6c..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 @@ -25,10 +30,13 @@ #include "ModelFFT.h" #include "ModelGrid.h" #include "ModelScaling.h" +#include "ModelValidation.h" #include "RigidBodyRefine.h" #ifdef JFJOCH_USE_CUDA #include "MolecularReplacementGPU.h" #include "RigidBodyGPU.h" +#include "RigidBodyGPUEngine.h" +#include "ModelStructureFactorsGPU.h" #include "../../common/CUDAWrapper.h" #include "../../common/JFJochException.h" #endif @@ -48,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. @@ -60,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; @@ -71,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(); @@ -89,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]; } @@ -141,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; @@ -157,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; @@ -201,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 @@ -235,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); @@ -271,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; @@ -283,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) @@ -621,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(); @@ -647,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}; @@ -655,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); @@ -666,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); @@ -700,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; @@ -728,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, @@ -744,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) { @@ -859,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 @@ -901,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); @@ -928,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)); } }; @@ -1022,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; @@ -1033,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(), @@ -1091,3 +1819,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..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,11 +263,20 @@ 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, 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/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/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. diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index f57c05fab..a71e3a743 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -2591,7 +2591,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"})); @@ -2613,22 +2616,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); @@ -2641,17 +2647,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/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/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/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..addf45aea --- /dev/null +++ b/tools/jfjoch_mr_replay.cpp @@ -0,0 +1,112 @@ +// 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"; + 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 + + 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} {: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); + } + return 0; +}