diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 714909a9c..96c20d134 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -981,9 +981,10 @@ GPU where there is one (`MODEL_PRIOR_SEARCH_DEVICE=`). `MODEL_PRIOR_DECISION=` i - `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; - `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 as a `--model`'s enantiomorph is: - model validation is run on the placed model, and the model has to fit against its null and the - anomalous map at its atoms must not be inverted (`MODEL_PRIOR_HAND_CHECK=`, + 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) + and model validation's anomalous difference map, read at the placed model's atoms, must not be + inverted (`MODEL_PRIOR_HAND_CHECK=`, `MODEL_PRIOR_HAND_FIT_SIGMA=`, `MODEL_PRIOR_ANOMALOUS_INVERTED=`). Where it is not, the decision reads `ENANTIOMORPH_VETOED` and the label stays; - `SPACE_GROUP` - another group with the answer's rotations; written as `_mr.mtz`/`.cif` from diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index ea927418c..e8ba4d186 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -31,6 +31,7 @@ #include // Scaling (bulk solvent + anisotropic B) #include // Ccp4 map I/O #include // Mtz (map-coefficient output) +#include // calculate_data_statistics #include "../../common/CorrelationCoefficient.h" #include "../../common/JFJochMath.h" // PI (M_PI is not standard, and MSVC does not define it) @@ -374,6 +375,93 @@ double frame_probe_r(const gemmi::Structure &st, const gemmi::SpaceGroup *sg, return scaling.calculate_r_factor(); } +// The anomalous difference map's coefficients on the reflections of `fmodel`: F(+) - F(-) carried on the +// model phase turned back by 90 degrees. +// +// Read from the merged reflections as they came in, and carry each one into the model's frame here: the +// hand each Bijvoet difference belongs to is a property of the frame the merge was made in, and both +// operators can change it. F(+) and F(-) are attached to the + index of the Friedel ASU of that frame, so +// an anomalous merge - which keeps each mate as a row of its own, both carrying the same pair - is read on +// its + rows only. Taking the - rows as well would give one reflection both signs of its difference, and +// the last row written would decide. +gemmi::AsuData> anomalous_coefficients(const std::vector &merged, + const gemmi::SpaceGroup *data_sg, + const gemmi::Op &indexing_op, + const gemmi::AsuData> &fmodel) { + const gemmi::ReciprocalAsu data_asu(data_sg); + const gemmi::GroupOps data_gops = data_sg->operations(); + const gemmi::ReciprocalAsu asu(fmodel.spacegroup_); + const gemmi::GroupOps gops = fmodel.spacegroup_->operations(); + + std::unordered_map danom_by_hkl; + for (const MergedReflection &r : merged) { + if (!std::isfinite(r.F_plus) || !std::isfinite(r.F_minus)) + continue; + gemmi::Op::Miller h{{r.h, r.k, r.l}}; + if (data_gops.is_reflection_centric(h)) // a centric reflection has no anomalous difference + continue; + if (!data_asu.to_asu_sign(h, data_gops).second) + continue; + if (!(indexing_op == gemmi::Op::identity())) + h = indexing_op.apply_to_hkl(h); + const auto [hasu, plus] = asu.to_asu_sign(h, gops); + danom_by_hkl[hkl_key(hasu)] = plus ? r.F_plus - r.F_minus : r.F_minus - r.F_plus; + } + + gemmi::AsuData> mapanom; + mapanom.unit_cell_ = fmodel.unit_cell_; + mapanom.spacegroup_ = fmodel.spacegroup_; + for (const auto &hv : fmodel.v) { + const auto it = danom_by_hkl.find(hkl_key(hv.hkl)); + if (it == danom_by_hkl.end()) + continue; + const auto phi = static_cast(std::arg(hv.value) - PI / 2); + mapanom.v.push_back({hv.hkl, it->second * std::polar(1.0f, phi)}); + } + return mapanom; +} + +// The anomalous map read at every atom of `st`, in units of its RMS: the strongest sites, the mean at the +// anomalous scatterers, and whether the map is inverted at the model's atoms. +void read_anomalous_map(const gemmi::Grid &grid, double rms, gemmi::Structure &st, + ModelValidationResult &result) { + std::vector sites; + double scatterer_sum = 0; + for (gemmi::Model &m : st.models) + for (gemmi::Chain &ch : m.chains) + for (gemmi::Residue &r : ch.residues) + for (gemmi::Atom &a : r.atoms) { + if (a.is_hydrogen()) // hydrogen scatters no anomalous signal + continue; + sites.push_back({fmt::format("{} {} {}{}", a.name, r.name, ch.name, r.seqid.str()), + rms > 0 ? grid.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER) / rms : 0.0}); + if (a.element.atomic_number() >= ANOMALOUS_SCATTERER_MIN_Z) { + scatterer_sum += sites.back().sigma; + result.anomalous_scatterers++; + } + } + if (result.anomalous_scatterers > 0) + result.anomalous_scatterer_mean_sigma = scatterer_sum / result.anomalous_scatterers; + std::sort(sites.begin(), sites.end(), [](const auto &x, const auto &y) { return x.sigma > y.sigma; }); + + // A model and a dataset in opposite hands turn every anomalous peak into a trough, so a map whose + // deepest hole at an atom is both deep and deeper than its highest peak says the two disagree about + // the hand. That is worth reporting: it is real evidence about the crystal, and the alternative - + // reindexing until the two agree - would erase it. + if (!sites.empty()) { + const auto &deepest = sites.back(); + if (deepest.sigma < -ANOMALOUS_INVERSION_SIGMA && -deepest.sigma > sites.front().sigma) { + result.anomalous_hands_disagree = true; + result.anomalous_deepest_site = deepest.label; + result.anomalous_deepest_sigma = deepest.sigma; + } + } + + if (sites.size() > MAX_ANOMALOUS_SITES) + sites.resize(MAX_ANOMALOUS_SITES); + result.anomalous_sites = std::move(sites); +} + #ifdef JFJOCH_USE_CUDA // The model's structure factors to the data's resolution, and the maps, on the GPU - or not, decided here, // once, before any of them is computed, and never revisited: nothing is moved to the CPU part way through. @@ -1500,84 +1588,13 @@ ModelValidationResult Validate(const std::vector &merged, // rather than leaving a list of coordinates for someone to look up. // Following ANODE, Thorn & Sheldrick (2011) J. Appl. Cryst. 44, 1285-1287 std::future anomalous_map = std::async(std::launch::async, [&] { - // Read from the merged reflections as they came in, and carry each one into the model's frame - // here: the hand each Bijvoet difference belongs to is a property of the frame the merge was - // made in, and both operators can change it. F(+) and F(-) are attached to the + index of the - // Friedel ASU of that frame, so an anomalous merge - which keeps each mate as a row of its own, - // both carrying the same pair - is read on its + rows only. Taking the - rows as well would - // give one reflection both signs of its difference, and the last row written would decide. - const gemmi::SpaceGroup *data_sg = data_space_group != nullptr ? data_space_group : sg; - const gemmi::ReciprocalAsu data_asu(data_sg); - const gemmi::GroupOps data_gops = data_sg->operations(); - - std::unordered_map danom_by_hkl; - for (const MergedReflection &r : merged) { - if (!std::isfinite(r.F_plus) || !std::isfinite(r.F_minus)) - continue; - gemmi::Op::Miller h{{r.h, r.k, r.l}}; - if (data_gops.is_reflection_centric(h)) // a centric reflection has no anomalous difference - continue; - if (!data_asu.to_asu_sign(h, data_gops).second) - continue; - if (!(result.indexing_op == gemmi::Op::identity())) - h = result.indexing_op.apply_to_hkl(h); - const auto [hasu, plus] = asu.to_asu_sign(h, gops); - danom_by_hkl[hkl_key(hasu)] = plus ? r.F_plus - r.F_minus : r.F_minus - r.F_plus; - } - - gemmi::AsuData> mapanom; - mapanom.unit_cell_ = ucell; - mapanom.spacegroup_ = sg; - for (const auto &hv : fmodel.v) { - const auto it = danom_by_hkl.find(hkl_key(hv.hkl)); - if (it == danom_by_hkl.end()) - continue; - const auto phi = static_cast(std::arg(hv.value) - PI / 2); - mapanom.v.push_back({hv.hkl, it->second * std::polar(1.0f, phi)}); - } + gemmi::AsuData> mapanom = + anomalous_coefficients(merged, data_space_group != nullptr ? data_space_group : sg, + result.indexing_op, fmodel); result.anomalous_pairs = static_cast(mapanom.v.size()); - if (!mapanom.v.empty()) { const gemmi::Grid grid = map_from_coefficients(mapanom, sf_gpu); - const double rms = write_ccp4(grid, output_prefix + "_anom.ccp4"); - std::vector sites; - double scatterer_sum = 0; - for (gemmi::Model &m : st.models) - for (gemmi::Chain &ch : m.chains) - for (gemmi::Residue &r : ch.residues) - for (gemmi::Atom &a : r.atoms) { - if (a.is_hydrogen()) // hydrogen scatters no anomalous signal - continue; - sites.push_back({fmt::format("{} {} {}{}", a.name, r.name, ch.name, - r.seqid.str()), - rms > 0 ? grid.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER) / rms - : 0.0}); - if (a.element.atomic_number() >= ANOMALOUS_SCATTERER_MIN_Z) { - scatterer_sum += sites.back().sigma; - result.anomalous_scatterers++; - } - } - if (result.anomalous_scatterers > 0) - result.anomalous_scatterer_mean_sigma = scatterer_sum / result.anomalous_scatterers; - std::sort(sites.begin(), sites.end(), - [](const auto &x, const auto &y) { return x.sigma > y.sigma; }); - - // A model and a dataset in opposite hands turn every anomalous peak into a trough, so a - // map whose deepest hole at an atom is both deep and deeper than its highest peak says - // the two disagree about the hand. That is worth reporting: it is real evidence about - // the crystal, and the alternative - reindexing until the two agree - would erase it. - if (!sites.empty()) { - const auto &deepest = sites.back(); - if (deepest.sigma < -ANOMALOUS_INVERSION_SIGMA && -deepest.sigma > sites.front().sigma) { - result.anomalous_hands_disagree = true; - result.anomalous_deepest_site = deepest.label; - result.anomalous_deepest_sigma = deepest.sigma; - } - } - - if (sites.size() > MAX_ANOMALOUS_SITES) - sites.resize(MAX_ANOMALOUS_SITES); - result.anomalous_sites = std::move(sites); + read_anomalous_map(grid, write_ccp4(grid, output_prefix + "_anom.ccp4"), st, result); } }); @@ -1794,6 +1811,79 @@ const gemmi::SpaceGroup *AdoptModelFrame(const ModelValidationResult &validation return space_group; } +ModelValidationResult ReadModelAnomalousMap(const std::vector &merged, + const gemmi::Structure &placed, const gemmi::SpaceGroup &data_space_group, + double k_sol, double b_sol, Logger &logger, size_t nthreads) { + ModelValidationResult result; + gemmi::Structure st = placed; + st.setup_cell_images(); + const gemmi::SpaceGroup *sg = st.find_spacegroup(); + double d_min = 0.0; + for (const MergedReflection &r : merged) + if (r.d > 0 && (d_min == 0.0 || r.d < d_min)) + d_min = r.d; + if (sg == nullptr || d_min <= 0.0) { + result.failure_reason = "no space group or no resolution"; + return result; + } + + // F_calc and F_mask to the data's resolution, as the validation computes them (compute_model_factors). + gemmi::AsuData> fcalc, fmask; + ModelStructureFactorsGPU *gpu = nullptr; +#ifdef JFJOCH_USE_CUDA + std::unique_ptr engine; + if (get_gpu_count() > 0) + engine = StructureFactorsGPU(RigidBodyGPUEngine::CurrentDevice(), st.models[0], st.cell, *sg, d_min, 0, logger); + gpu = engine.get(); + if (gpu != nullptr) { + try { + gpu->Compute(st.models[0], fcalc, fmask); + } catch (const JFJochException &e) { + throw RigidBodyGPUFailure(e.what()); + } + } +#else + (void) logger; +#endif + if (gpu == nullptr) { + gemmi::DensityCalculator dc; + dc.d_min = d_min; + dc.rate = 1.5; + dc.set_grid_cell_and_spacegroup(st); + dc.set_refmac_compatible_blur(st.models[0]); + gemmi::Grid mask_grid; + mask_grid.unit_cell = dc.grid.unit_cell; + mask_grid.spacegroup = dc.grid.spacegroup; + mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up); + const std::vector orbit_leaders = OrbitLeaders(mask_grid, nthreads); + PutMaskOnGrid(mask_grid, st.models[0], orbit_leaders, nthreads); + fmask = MapToFPhi(mask_grid).prepare_asu_data(d_min, 0); + PutModelDensityOnGrid(dc, st.models[0], orbit_leaders, nthreads); + fcalc = MapToFPhi(dc.grid).prepare_asu_data(dc.d_min, dc.blur, false, false, false); + } + + // With the bulk solvent the placement was scored with: the overall scale and the anisotropic B are + // real and positive and leave the phases where they are, so the solvent is all of fit_model's scaling + // that the map's phases see. + gemmi::AsuData> fmodel = fcalc; + gemmi::Scaling scaling(st.cell, sg); + scaling.use_solvent = true; + scaling.k_sol = k_sol; + scaling.b_sol = b_sol; + scaling.scale_data(fmodel, &fmask); + + gemmi::AsuData> mapanom = + anomalous_coefficients(merged, &data_space_group, gemmi::Op::identity(), fmodel); + result.anomalous_pairs = static_cast(mapanom.v.size()); + if (!mapanom.v.empty()) { + const gemmi::Grid grid = map_from_coefficients(mapanom, gpu); + // The RMS write_ccp4 reads the map in. + read_anomalous_map(grid, gemmi::calculate_data_statistics(grid.data).rms, st, result); + } + result.ok = true; + return result; +} + void KeepModelVerdict(ModelValidationResult &remade, const ModelValidationResult &decided) { remade.fit_tested = decided.fit_tested; remade.model_fits = decided.model_fits; diff --git a/image_analysis/structure_refinement/ModelValidation.h b/image_analysis/structure_refinement/ModelValidation.h index adc686bf2..ebd6b26cf 100644 --- a/image_analysis/structure_refinement/ModelValidation.h +++ b/image_analysis/structure_refinement/ModelValidation.h @@ -293,6 +293,17 @@ const gemmi::SpaceGroup *AdoptModelFrame(const ModelValidationResult &validation // because they are on the axes the files are written on. void KeepModelVerdict(ModelValidationResult &remade, const ModelValidationResult &decided); +// The anomalous difference map of ValidateAgainstModel, and nothing else, for a model that is already +// placed and scored in the data's cell (`placed` carries the cell and the model's group; k_sol and b_sol +// are the bulk solvent it was scored with): F_calc and F_mask to the data's resolution, the map +// F(+) - F(-) on the phases of F_calc + solvent, and its reading at the model's atoms - anomalous_pairs, +// the scatterers' mean, the sites and anomalous_hands_disagree, by the same code and the same rule as the +// validation's. No placement, no scale fit, no null, no file. For a caller that has already placed the +// model and judged its fit (the --af-model hand check). On the GPU where there is one. +ModelValidationResult ReadModelAnomalousMap(const std::vector &merged, + const gemmi::Structure &placed, const gemmi::SpaceGroup &data_space_group, + double k_sol, double b_sol, Logger &logger, size_t nthreads); + // Reference intensities computed from an atomic model: |F_model|^2, with a flat bulk-solvent // contribution at the standard constants, for every reflection to `d_min`. This is the same // information a reference MTZ carries, and it is what lets a model break the indexing ambiguity of a diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp index ed7719c3c..428cf31bf 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.cpp +++ b/image_analysis/structure_refinement/MolecularReplacement.cpp @@ -117,6 +117,7 @@ int good_even_size(int n) { struct PlacementScore { double r_work = 1, r_free = 1; + double k_sol = NAN, b_sol = NAN; }; gemmi::AsuData> working_amplitudes(const std::vector &merged, @@ -193,6 +194,8 @@ PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGro } } PlacementScore out; + out.k_sol = scaling.k_sol; + out.b_sol = scaling.b_sol; out.r_work = den_w > 0 ? num_w / den_w : 1; out.r_free = den_f > 0 ? num_f / den_f : 1; return out; @@ -936,7 +939,7 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path // body moves a search placement by a few degrees and a fraction of an angstrom, and the // clash fraction moves by a percent or two with it. double clash = ClashFraction(*placed); - PlacementScore score{NAN, NAN}; + PlacementScore score{NAN, NAN, NAN, NAN}; if (one_copy && clash <= 2 * CLASH_MAX) { if (engines.get() == nullptr) engines = make_engines(*placed, sg, *merge.merged, d_min, 1, logger); @@ -952,6 +955,8 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path 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]]; } diff --git a/image_analysis/structure_refinement/MolecularReplacement.h b/image_analysis/structure_refinement/MolecularReplacement.h index 3bc628534..fa5031d17 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.h +++ b/image_analysis/structure_refinement/MolecularReplacement.h @@ -161,6 +161,7 @@ struct MRCandidate { double tf_z = 0; double clash_fraction = 0; 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; diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index b822754c6..353a46c35 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -2648,7 +2648,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, result.model_prior_file_prefix.empty() ? std::string("_mr") : result.model_prior_file_prefix))); Add(s, Prose(fmt::format( " Seconds: model {:.2f}, transform {:.2f}, rotation {:.2f}, translation {:.2f}, rigid body and\n" - " R-free {:.2f}, null {:.2f}, hand check {:.2f}.", mp.t_model, mp.t_transform, mp.t_rotation, + " R-free {:.2f}, null {:.2f}, hand check {:.3f}.", mp.t_model, mp.t_transform, mp.t_rotation, mp.t_translation, mp.t_rescore, mp.t_null, mp.hand_check.seconds), true)); } doc.sections.push_back(std::move(s)); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index bb1f623be..a3391b14d 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -20,7 +20,6 @@ #include #include #include -#include #include #include #include @@ -4002,22 +4001,23 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { if (mr.chosen >= 0) { const MRCandidate &c = mr.candidates[mr.chosen]; if (mr.decision == "ENANTIOMORPH") { - // Taken into the answer as --model takes a model's enantiomorph: model validation runs on - // the placed model, in the group it chose, against the answer's merge, and its verdict - // decides - the model has to fit against its own null, and the anomalous map at its atoms - // must not come out inverted. The validation's files go to a scratch directory, removed after. + // 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(); - const std::filesystem::path scratch = config_.output_prefix + "_mr_hand"; - std::filesystem::create_directories(scratch); - const std::string tmp = (scratch / "m").string(); Logger quiet = Logger::Buffered(); - WritePlacedModel(*c.placed, *result.consensus_cell, c.space_group, tmp, quiet); - const auto data_sg = experiment_.GetGemmiSpaceGroup(); - const ModelValidationResult v = ValidateAgainstModel( - sm.merged, *result.consensus_cell, tmp + "_model.cif", tmp, quiet, - data_sg ? &*data_sg : nullptr, /*probe_indexing_ambiguity=*/false, - static_cast(config_.nthreads), experiment_.GetWavelength_A()); - std::filesystem::remove_all(scratch); + 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; @@ -4031,7 +4031,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { 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 - {} ({:.1f} s)", + 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)", @@ -4052,7 +4052,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { twin_file->space_group = *other; } logger.Info("Model prior: the answer is relabelled {}, the enantiomorph the placed model chose " - "(fit {:+.1f} sd over its null, anomalous map at its atoms not inverted, {:.1f} s)", + "(fit {:+.1f} sd over its null, anomalous map at its atoms not inverted, {:.3f} s)", c.space_group.xhm(), mr.hand_check.fit_sigma, mr.hand_check.seconds); } else if (mr.decision == "SPACE_GROUP") { ModelPriorFile f{sm.merged, *result.consensus_cell, c.space_group, std::nullopt, {}, c.placed};