diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index db415c610..7aa5a6692 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -41,6 +41,7 @@ #include "ModelFFT.h" #include "ModelGrid.h" // PutModelDensityOnGrid, PutMaskOnGrid #include "RigidBodyRefine.h" +#include "MolecularReplacement.h" // the warm-up of a prediction #include "SigmaA.h" #ifdef JFJOCH_USE_CUDA #include "ModelStructureFactorsGPU.h" @@ -1964,8 +1965,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & const ModelValidationSchedule &schedule, const gemmi::Op *twin_law, const gemmi::Structure *placed) { - if (const auto w = LastWarmUp(); w.valid()) - w.wait(); + WaitForModelWarmUp(); #ifdef JFJOCH_USE_CUDA // A CUDA failure ends the validation, not the run, and nothing is moved to the CPU part way through: // the validation reports why it did not finish, and a validation that did not finish decides nothing, @@ -1989,6 +1989,11 @@ ModelValidationResult ValidateAgainstModel(const std::vector & #endif } +void WaitForModelWarmUp() { + if (const auto w = LastWarmUp(); w.valid()) + w.wait(); +} + void WarmUpModelValidation(std::shared_ptr model, size_t nthreads, Logger logger) { #ifdef JFJOCH_USE_CUDA if (!model || get_gpu_count() == 0) @@ -2002,6 +2007,26 @@ void WarmUpModelValidation(std::shared_ptr model, size_t nthr // Held for as long as it runs: work beside a run waits for it where the card is too small for both. GPUWorkBeside beside; Logger quiet = Logger::Buffered(); + if (const SearchModel search = PrepareSearchModel(*model); search.failure_reason.empty() && search.predicted) { + std::string failure; + try { + const MolecularReplacementWarmUp w = WarmUpMolecularReplacement(*model, nthreads); + failure = w.failure_reason; + if (w.placed) { + ModelValidationSchedule placed_elsewhere; + placed_elsewhere.fit_test = false; + failure = Validate(w.merged, w.cell, *model, quiet, &w.space_group, false, nthreads, 0.0, {}, + placed_elsewhere, nullptr, true, w.placed.get()).failure_reason; + } + } catch (const std::exception &e) { + cuda_clear_error(); + failure = e.what(); + } + logger.Info("Molecular replacement warm-up on the GPU: {:.2f} s{}", + std::chrono::duration(std::chrono::steady_clock::now() - start).count(), + failure.empty() ? "" : " (did not finish: " + failure + ")"); + return; + } const double d_min = model->structure.resolution > 0.0 ? model->structure.resolution : 2.0; std::vector fake = ModelReferenceIntensities(*model, std::nullopt, nullptr, d_min, quiet); for (size_t i = 0; i < fake.size(); i++) { diff --git a/image_analysis/structure_refinement/ModelValidation.h b/image_analysis/structure_refinement/ModelValidation.h index 0eebb653b..9551ff0f9 100644 --- a/image_analysis/structure_refinement/ModelValidation.h +++ b/image_analysis/structure_refinement/ModelValidation.h @@ -416,4 +416,8 @@ std::vector ModelReferenceIntensities(const RugnuxModel &model // calls it when a model is uploaded, before the collection starts. ValidateAgainstModel waits for a // warm-up still running, so a real validation never shares the cards with one and computes exactly // what it would have without it. Its time goes to `logger` when it ends. No-op without a GPU. +// A prediction, which a run places by molecular replacement and validates only as placed there, is warmed up by a +// search instead (WarmUpMolecularReplacement) and the validation of what it placed. void WarmUpModelValidation(std::shared_ptr model, size_t nthreads, Logger logger); +// Returns once the warm-up last started has ended (at once where there is none). +void WaitForModelWarmUp(); diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp index 63b7689a0..816f0dc18 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.cpp +++ b/image_analysis/structure_refinement/MolecularReplacement.cpp @@ -3713,8 +3713,10 @@ std::vector SplitTranslationFunction(const std::vector &merges, - Logger &logger, size_t nthreads, size_t gpu_resident_bytes) { +namespace { + +MolecularReplacementResult run_molecular_replacement(const RugnuxModel &input, const std::vector &merges, + Logger &logger, size_t nthreads, size_t gpu_resident_bytes) { const auto t0 = std::chrono::steady_clock::now(); // What the search may hold on the card: the budget model validation keeps to. What the stream-ordered pool // holds with no buffer in it - freed by the steps before - goes back to the driver first, as before a @@ -3822,6 +3824,135 @@ MolecularReplacementResult RunMolecularReplacement(const RugnuxModel &input, con return result; } +} // namespace + +MolecularReplacementResult RunMolecularReplacement(const RugnuxModel &model, const std::vector &merges, + Logger &logger, size_t nthreads, size_t gpu_resident_bytes) { + WaitForModelWarmUp(); + return run_molecular_replacement(model, merges, logger, nthreads, gpu_resident_bytes); +} + +std::optional CarryMolecularReplacement(const MolecularReplacementResult &early, + const MRMerge &searched, const MRMerge &merge, + Logger &logger, size_t nthreads, + size_t gpu_resident_bytes) { + if (!early.ok || early.chosen < 0 || early.candidates[early.chosen].merge != 0) { + logger.Info("Molecular replacement: the search on the {} placed nothing ({}) - searched again", + searched.name, early.ok ? early.decision : early.failure_reason); + return std::nullopt; + } + if (searched.space_group.xhm() != merge.space_group.xhm() || !merge.cell.is_close(searched.cell, 0.01, 1.0)) { + logger.Info("Molecular replacement: the search on the {} was made in {} ({:.2f} {:.2f} {:.2f}), this merge is in " + "{} ({:.2f} {:.2f} {:.2f}) - searched again", searched.name, searched.space_group.xhm(), + searched.cell.a, searched.cell.b, searched.cell.c, merge.space_group.xhm(), merge.cell.a, + merge.cell.b, merge.cell.c); + return std::nullopt; + } + const auto t0 = std::chrono::steady_clock::now(); + const MRCandidate &c = early.candidates[early.chosen]; + const gemmi::UnitCell cell(merge.cell); + const gemmi::SpaceGroup &sg = c.space_group; // the group it was placed in, one of the merge's (SameRotationGroups) + double data_d_min = 1e9; + for (const MergedReflection &r : *merge.merged) + if (r.d > 0) + data_d_min = std::min(data_d_min, r.d); + const double d_min = std::max(SCORE_D_MIN, data_d_min); + // Every copy at its fractional position in this cell. + std::vector body = SplitCopies(*c.placed, c.copies); + for (gemmi::Structure &b : body) { + std::vector pos = ModelPositions(b.models[0]); + for (gemmi::Position &x : pos) + x = cell.orthogonalize(b.cell.fractionalize(x)); + SetModelPositions(b.models[0], pos); + b.cell = cell; + b.spacegroup_hm = sg.xhm(); + b.setup_cell_images(); + } + RigidBodyGPUPool *pool = nullptr; + ScoreEngines scorers; +#ifdef JFJOCH_USE_CUDA + // One rigid-body engine and the score engine, on this card, as a search's own (size_merge_gpu). + std::unique_ptr engines; + if (get_gpu_count() > 0) { + const int device = RigidBodyGPUEngine::CurrentDevice(); + const size_t budget = ModelGPUBudget(get_gpu_total_memory(device), gpu_resident_bytes); + size_t observations = 0; + for (const MergedReflection &r : *merge.merged) + if (!std::isnan(r.F) && r.d >= RigidBodyLadder(d_min).back()) + ++observations; + Logger quiet = Logger::Buffered(); + const gemmi::Structure all = AssembleCopies(body); + engines = RigidBodyGPUPool::Plan(body[0].models[0], cell, sg, d_min, observations, quiet); + if (engines && engines->Cover(all.models[0], sg, quiet) && engines->EngineBytes() <= budget / 2) { + engines->Reserve({1}); + pool = engines.get(); + } + scorers.SetRoom({{device, budget - (pool != nullptr ? engines->EngineBytes() : 0)}}, device, true); + } +#endif + for (size_t k = 0; k < body.size(); k++) { + std::vector rest; + for (size_t o = 0; o < body.size(); o++) + if (o != k) + rest.push_back(body[o]); + if (rest.empty()) { + refine_rigid_body(body[k], sg, *merge.merged, d_min, pool, nthreads); + } else { + const gemmi::Structure fixed = AssembleCopies(rest); + refine_rigid_body(body[k], sg, *merge.merged, d_min, pool, nthreads, nullptr, &fixed.models[0]); + } + } + auto placed = std::make_shared(AssembleCopies(body)); + const PlacementScore score = score_placement(*placed, sg, *merge.merged, d_min, nthreads, &scorers); + const bool holds = score.r_free <= c.r_free + R_FREE_MARGIN; + logger.Info("Molecular replacement: the search on the {} ({}): its placement in {} ({} copies, R-free {:.3f} there) refined against " + "this merge: R-work {:.3f} R-free {:.3f} ({:.2f} s) - {}", searched.name, + early.decision, sg.xhm(), c.copies, c.r_free, score.r_work, score.r_free, + seconds_since(t0), holds ? "it holds, not searched again" : "it does not hold, searched again"); + if (!holds) + return std::nullopt; + MolecularReplacementResult out = early; + MRCandidate &chosen = out.candidates[out.chosen]; + out.r_free_searched = chosen.r_free; + chosen.placed = placed; + chosen.r_work = score.r_work; + chosen.r_free = score.r_free; + chosen.k_sol = score.k_sol; + chosen.b_sol = score.b_sol; + out.d_min_score = d_min; + return out; +} + +MolecularReplacementWarmUp WarmUpMolecularReplacement(const RugnuxModel &model, size_t nthreads) { + MolecularReplacementWarmUp out; + const SearchModel search = PrepareSearchModel(model); + if (!search.failure_reason.empty()) { + out.failure_reason = search.failure_reason; + return out; + } + const double edge = 4.0 * search.radius_A; + out.cell = UnitCell{static_cast(edge), static_cast(edge + 2.0), static_cast(edge + 4.0), + 90.0f, 100.0f, 90.0f}; + out.space_group = *gemmi::find_spacegroup_by_name("P 1 21 1"); + RugnuxModel crystal; + crystal.name = model.name; + crystal.structure = PlaceSearchModel(search, gemmi::Mat33(), gemmi::Fractional(0.1, 0.2, 0.3), + gemmi::UnitCell(out.cell), out.space_group); + Logger quiet = Logger::Buffered(); + out.merged = ModelReferenceIntensities(crystal, std::nullopt, nullptr, 2.5, quiet); + for (size_t i = 0; i < out.merged.size(); i++) { + out.merged[i].F = std::sqrt(out.merged[i].I); + out.merged[i].rfree_flag = i % 20 == 0; + } + const MolecularReplacementResult mr = run_molecular_replacement( + model, {{"MERGE", &out.merged, out.cell, out.space_group}}, quiet, nthreads, 0); + if (mr.chosen >= 0) + out.placed = mr.candidates[mr.chosen].placed; + else + out.failure_reason = mr.ok ? mr.decision : mr.failure_reason; + return out; +} + void CheckEnantiomorphHand(MolecularReplacementResult &result, const std::vector &merged, const gemmi::SpaceGroup &data_space_group, Logger &logger, size_t nthreads, size_t gpu_resident_bytes) { diff --git a/image_analysis/structure_refinement/MolecularReplacement.h b/image_analysis/structure_refinement/MolecularReplacement.h index 8c90974a3..d56d402f0 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.h +++ b/image_analysis/structure_refinement/MolecularReplacement.h @@ -374,6 +374,9 @@ struct MolecularReplacementResult { // Seconds per step. double t_model = 0, t_transform = 0, t_rotation = 0, t_translation = 0, t_rescore = 0, t_null = 0; double t_domains = 0; // the PAE domains' own searches, wall clock (they run side by side) + // Searched on another merge of the crystal and carried over to this one (CarryMolecularReplacement): the + // chosen placement's R-free on the merge it was found on, NaN where this merge was searched itself. + double r_free_searched = NAN; }; // What a molecular-replacement search puts on its GPU, for one merge: the molecular transform the rotation and @@ -450,9 +453,34 @@ bool PackMolecularTransformGPU(const MolecularReplacementGPUNeed &need, size_t p // 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. // `gpu_resident_bytes`: what is held on the card meanwhile by what came before (a merge engine still at work). +// Waits for a model warm-up still running (WarmUpModelValidation), as a validation does. MolecularReplacementResult RunMolecularReplacement(const RugnuxModel &model, const std::vector &merges, Logger &logger, size_t nthreads, size_t gpu_resident_bytes = 0); +// A search made on another merge of the same crystal - `searched`, the geometry pre-pass's in a live run - taken +// over to `merge` where it holds there: the same space group, a cell within 1%, and the chosen placement - every +// copy refined by the rigid body beside the others against `merge`, and scored there - fitting it within +// R_FREE_MARGIN in R-free of how it fitted the merge it was found on. Its decision stands as it was taken; the +// placement and its R-factors are `merge`'s. Empty where it does not hold, and `merge` is searched itself. +std::optional CarryMolecularReplacement(const MolecularReplacementResult &early, + const MRMerge &searched, const MRMerge &merge, + Logger &logger, size_t nthreads, + size_t gpu_resident_bytes = 0); + +// A search for a prediction ahead of its run, which pays for what the GPU loads on first use - the kernels of the +// transform, the rotation and translation functions, the rigid body and the scores, on every card the fits are +// spread over. The crystal is made up from the model alone: its search model in one orientation in a P 1 21 1 cell +// four of its radii on each edge, and that crystal's own intensities to 2.5 A as the merge. The search is +// RunMolecularReplacement's, without the wait; `placed` is what it placed, null where it placed nothing. +struct MolecularReplacementWarmUp { + std::vector merged; + UnitCell cell; + gemmi::SpaceGroup space_group; + std::shared_ptr placed; + std::string failure_reason; +}; +MolecularReplacementWarmUp WarmUpMolecularReplacement(const RugnuxModel &model, 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 diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 03eb8bab2..776aa617c 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -1227,6 +1227,11 @@ class Rugnux { }; std::shared_ptr speculative_probe_; bool speculate_geometry_probe_ = false; + // fast_live, with a prediction (or --mr): its molecular replacement, started on the geometry pre-pass's merge + // and run beside the second pass; the answer's merge takes it over where it holds there + // (CarryMolecularReplacement) and searches itself where not. RugnuxScaleMerge.cpp. + struct EarlyMolecularReplacement; + std::shared_ptr early_mr_; // Whether this run is a probe pass running beside the main line under a GPUWorkBeside - the // speculative probe, the rotation-scale walk's indexing probes. Its indexer pools then wait for the // main line's image pass, never for that work, being that work themselves (IndexerThreadPool). diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 54951ad7f..661941763 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -347,6 +347,47 @@ namespace { } } +// The molecular replacement of a live run started on the geometry pre-pass's merge (Rugnux.h): the merge it +// searches - the pre-pass's P1 merge averaged into a group of its own search, on the axes the answer's files +// would be written on - and the search, its log held until the answer's merge reads it. +struct Rugnux::EarlyMolecularReplacement { + std::vector merged; + UnitCell cell; + gemmi::SpaceGroup space_group; + Logger log = Logger::Buffered(); + std::future run; +}; + +// A P1 merge averaged into the asymmetric unit of `sg`, each reflection weighted by 1/sigma^2; the +// reflections `sg` makes absent left out. +static std::vector AverageIntoGroup(const std::vector &p1, + const gemmi::SpaceGroup &sg, const UnitCell &cell) { + const HKLKeyGenerator key_gen(/*merge_friedel=*/true, sg); + std::map, std::pair> sums; // sum of w I, sum of w + for (const MergedReflection &r : p1) { + if (!std::isfinite(r.I) || !(r.sigma > 0.0f) || key_gen.IsSystematicallyAbsent(r.h, r.k, r.l)) + continue; + const HKLKey key = key_gen(r.h, r.k, r.l); + const double w = 1.0 / (static_cast(r.sigma) * r.sigma); + auto &[wi, ws] = sums[{key.h, key.k, key.l}]; + wi += w * r.I; + ws += w; + } + const gemmi::UnitCell uc(cell); + std::vector out; + for (const auto &[hkl, sum] : sums) { + MergedReflection m; + m.h = hkl[0]; + m.k = hkl[1]; + m.l = hkl[2]; + m.I = static_cast(sum.first / sum.second); + m.sigma = static_cast(1.0 / std::sqrt(sum.second)); + m.d = static_cast(uc.calculate_d({{m.h, m.k, m.l}})); + out.push_back(m); + } + return out; +} + // Proper rotations of a crystal system's holohedry - what a lattice of that class carries, and the // bound SearchSpaceGroup caps its point-group candidates at. static size_t HolohedryRotations(gemmi::CrystalSystem system) { @@ -1534,6 +1575,72 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { if (geometry_prepass && rsm && experiment_.IsRotationIndexing()) post_refine_geometry(); + // Live, the molecular replacement of a prediction (or of --mr) starts here, on this pass's merge, and runs + // beside the second pass on a card of its own where there is one, with a quarter of the threads. It is a + // hypothesis: the answer's merge takes its placement over only where it holds there (CarryMolecularReplacement) + // and searches itself otherwise. Its group is a quick one of its own, from this merge's search to the + // resolution the placements are scored at - the pre-pass decides no group for the run - so the merge searched is + // this one averaged into that group, on the axes the answer's files would be written on. + if (geometry_prepass && config_.fast_live && config_.model && result.consensus_cell + && config_.reference_data.empty() && (config_.force_mr || PrepareSearchModel(*config_.model).predicted)) { + early_mr_ = std::make_shared(); + early_mr_->run = std::async(std::launch::async, + [early = early_mr_.get(), p1 = sm.merged, cell = *result.consensus_cell, fixed_sg = initial_sg, + lattice = end_msg.rotation_lattice_type, isa = result.error_model_isa, model = config_.model, + rfree_fraction = experiment_.GetScalingSettings().GetRfreeFraction(), + nthreads = static_cast(std::max(1, config_.nthreads / 4))]() mutable { + GPUWorkBeside beside; + if (const int32_t cards = get_gpu_count(); cards > 0) + pin_gpu(1 % cards); + gemmi::SpaceGroup sg; + if (fixed_sg) { + sg = *fixed_sg; + } else { + SearchSpaceGroupOptions o; + o.nthreads = nthreads; + o.merge_friedel = true; + o.d_min_limit_A = 3.0; + if (lattice) { + o.lattice_system = lattice->crystal_system; + o.lattice_centring = lattice->centering; + } + o.cell = gemmi::UnitCell(cell); + o.enumerate_all_settings = true; + o.enumerate_all_rotation_sets = true; + o.merge_isa = isa; + const SearchSpaceGroupResult found = SearchSpaceGroup(p1, o); + if (found.best_space_group) + sg = *found.best_space_group; + else if (found.point_group_representative) + sg = *found.point_group_representative; + else { + MolecularReplacementResult none; + none.failure_reason = "no space group from the pre-pass merge"; + return none; + } + } + const gemmi::Op cob = ChooseOutputSetting(cell, sg); + if (const gemmi::SpaceGroup *moved = SpaceGroupInBasis(sg, cob); + !(cob == gemmi::Op::identity()) && moved != nullptr) { + p1 = ReindexReflections(p1, HklOperator(cob)); + cell = CellInBasis(cell, cob); + sg = *moved; + } + early->merged = AverageIntoGroup(p1, sg, cell); + early->cell = cell; + early->space_group = sg; + AssignRfreeFlags(early->merged, sg, rfree_fraction, 500, cell, nthreads); + FrenchWilsonOptions fw_opts; + fw_opts.num_threads = nthreads; + fw_opts.anisotropy_hkl = AnisotropyPriorHKL(early->merged, gemmi::UnitCell(cell), &sg); + ApplyFrenchWilson(early->merged, sg, fw_opts); + early->log.Info("Molecular replacement on the geometry pre-pass's merge, in {} from its own search, beside " + "the second pass", sg.xhm()); + return RunMolecularReplacement(*model, {{"geometry pre-pass's merge", &early->merged, cell, sg}}, + early->log, nthreads); + }); + } + // For a de-novo search, gate the P1 scaling+merge (the pass that feeds the space-group search) at // >= 1 over thin resolution shells: noise-dominated high-res shells otherwise mislead // the determination, where a single artifact (a zinger the box sum integrates whole, or a @@ -4444,8 +4551,27 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { if (model_twin) merges.push_back({"TWIN_SUBGROUP_MERGE", &model_twin->merged, model_twin->cell, model_twin->space_group, model_twin->law}); - MolecularReplacementResult mr = RunMolecularReplacement(*config_.model, merges, logger, - static_cast(config_.nthreads), merge_on_card); + // The search started on the geometry pre-pass's merge, where there is one, taken over where it holds on + // this one. Not where the twin hypothesis' merge is searched too: that search is this merge's own. + std::optional carried; + if (early_mr_ && !model_twin) { + TimingMark early_mark(logger, "early molecular replacement wait"); + try { + MolecularReplacementResult early = early_mr_->run.get(); + early_mark.Stop(); + early_mr_->log.ReplayInto(logger); + carried = CarryMolecularReplacement( + early, {"geometry pre-pass's merge", &early_mr_->merged, early_mr_->cell, early_mr_->space_group}, + merges[0], logger, static_cast(config_.nthreads), merge_on_card); + } catch (const std::exception &e) { + logger.Warning("Molecular replacement on the geometry pre-pass's merge failed ({}) - searched again", + e.what()); + } + } + MolecularReplacementResult mr = carried ? std::move(*carried) + : RunMolecularReplacement(*config_.model, merges, logger, + static_cast(config_.nthreads), + merge_on_card); if (mr.chosen >= 0) { const MRCandidate &c = mr.candidates[mr.chosen]; if (mr.decision == "ENANTIOMORPH") diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index ea281a597..2a1ff9f88 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -3186,6 +3186,9 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { if (simulate_live) logger.Info("Simulated live collection: image i readable (i+1) x {:.1f} ms after the run starts", simulated_frame_time.count() / 1000.0); + // As the broker does when a model is uploaded, before the collection it is for starts. + if (simulate_live && model) + WarmUpModelValidation(model, static_cast(config.nthreads), logger); simulated_start = std::chrono::steady_clock::now(); Rugnux process(reader, experiment, *dataset->pixel_mask, config); diff --git a/tests/MolecularReplacementTest.cpp b/tests/MolecularReplacementTest.cpp index 438144573..696e8fd06 100644 --- a/tests/MolecularReplacementTest.cpp +++ b/tests/MolecularReplacementTest.cpp @@ -646,6 +646,54 @@ TEST_CASE("MolecularReplacement_RigidBodyBesideFixedCopy", "[MolecularReplacemen #endif } +// A search made on one merge, carried over to another merge of the same crystal in a cell 0.5% longer: the +// placement, a degree and half an angstrom out where it was found, is refined against the new merge onto the +// truth and holds there; in another group it is not carried at all, and the merge is to be searched itself. +TEST_CASE("MolecularReplacement_CarryToAnotherMerge", "[MolecularReplacement]") { + const SearchModel model = synthetic_model(); + const gemmi::SpaceGroup &sg = gemmi::get_spacegroup_by_name("P 21 21 21"); + const gemmi::UnitCell found_cell(52, 60, 70, 90, 90, 90), cell(52.26, 60.3, 70.35, 90, 90, 90); + const gemmi::Mat33 r = axis_angle({1, 2, 3}, 50); + const gemmi::Fractional t(0.17, 0.31, 0.09); + const gemmi::Structure truth = PlaceSearchModel(model, r, t, cell, sg); + auto merged = observed(truth); + for (size_t i = 0; i < merged.size(); i++) + merged[i].rfree_flag = i % 10 == 0; + + MolecularReplacementResult early; + early.ok = true; + early.chosen = 0; + early.decision = "AGREES"; + MRCandidate c; + c.space_group = sg; + c.copies = 1; + c.r_free = 0.30; + c.placed = std::make_shared(PlaceSearchModel( + model, axis_angle({0, 0, 1}, 1).multiply(r), gemmi::Fractional(t.x + 0.5 / 52, t.y, t.z), found_cell, sg)); + early.candidates.push_back(c); + + const UnitCell found_uc{52, 60, 70, 90, 90, 90}, uc{52.26f, 60.3f, 70.35f, 90, 90, 90}; + const MRMerge searched{"first merge", &merged, found_uc, sg}; + Logger logger("MolecularReplacement_CarryToAnotherMerge"); + const auto carried = CarryMolecularReplacement(early, searched, {"MERGE", &merged, uc, sg}, logger, 4); + REQUIRE(carried); + CHECK(carried->decision == "AGREES"); + CHECK(carried->r_free_searched == Catch::Approx(0.30)); + CHECK(carried->candidates[0].placed->cell.a == Catch::Approx(52.26)); + 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; + }; + CHECK(worst(*c.placed) > 0.5); + CHECK(worst(*carried->candidates[0].placed) < 0.2); + + const gemmi::SpaceGroup &other = gemmi::get_spacegroup_by_name("P 2 21 21"); + CHECK_FALSE(CarryMolecularReplacement(early, searched, {"MERGE", &merged, uc, other}, logger, 4)); +} + #ifdef JFJOCH_USE_CUDA TEST_CASE("MolecularReplacement_GPUMatchesCPU", "[MolecularReplacement]") { if (get_gpu_count() == 0)