diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index f5ca5d0a..935a7930 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -130,7 +130,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & dc.set_grid_cell_and_spacegroup(st); dc.set_refmac_compatible_blur(st.models[0]); dc.put_model_density_on_grid(st.models[0]); - gemmi::AsuData> fmodel = + gemmi::AsuData> fcalc = gemmi::transform_map_to_f_phi(dc.grid, true).prepare_asu_data(dc.d_min, dc.blur, false, false, false); // --- flat bulk-solvent mask -> Fmask --- @@ -144,107 +144,146 @@ ModelValidationResult ValidateAgainstModel(const std::vector & gemmi::AsuData> fmask = gemmi::transform_map_to_f_phi(mask_grid, true).prepare_asu_data(dc.d_min, 0); - // --- observed amplitudes into the model ASU, keyed by hkl (also remember free flag) --- gemmi::GroupOps gops = sg->operations(); gemmi::ReciprocalAsu asu(sg); - gemmi::AsuData> fobs; - fobs.unit_cell_ = ucell; - fobs.spacegroup_ = sg; - std::unordered_map> obs_by_hkl; // hkl -> (Fobs, is_free) - // Observed amplitudes are the French-Wilson |F| already computed at the end of the merge - // (MergedReflection.F), so the model R-free / maps use exactly the same amplitudes as the - // written reflection file. - for (const MergedReflection &r : obs) { - if (std::isnan(r.F)) continue; - gemmi::Miller h{{r.h, r.k, r.l}}; - if (!asu.is_in(h)) h = asu.to_asu(h, gops).first; - fobs.v.push_back({h, {r.F, 1.0f}}); - obs_by_hkl[hkl_key(h)] = {r.F, r.rfree_flag}; - } - fobs.ensure_asu(); - fobs.ensure_sorted(); + // --- fit the (scaled, solvent-corrected) model to one observed set and score it --- + // Factored into a lambda so we can probe indexing (merohedral) ambiguities: run the same scale + + // R computation on each reindexing of the observed reflections and keep the lowest-R-free one. + struct Fit { + gemmi::AsuData> fmodel, map2fofc, mapfofc; + std::unordered_map> obs_by_hkl; // hkl -> (Fobs, is_free) + double r_work = 1, r_free = 1, k_sol = 0, b_sol = 0, k_overall = 0; + int n_w = 0, n_f = 0, n_shells = 0; + }; + auto fit_model = [&](const std::vector &obs_in) -> Fit { + Fit out; + out.fmodel = fcalc; // copy the atomic structure factors; scaling mutates them in place - // --- scale Fmodel(+solvent) to Fobs: k_overall, anisotropic B, k_sol, b_sol --- - gemmi::Scaling scaling(ucell, sg); - scaling.use_solvent = true; - scaling.prepare_points(fmodel, fobs, &fmask); - scaling.fit_isotropic_b_approximately(); - scaling.fit_parameters(); - scaling.scale_data(fmodel, &fmask); // fmodel now holds the scaled, solvent-corrected Fmodel - - // --- per-resolution-shell scaling on top of the overall scaling --- - // Fit a smooth scale K(1/d^2) least-squares per shell (sum(Fo*Fc)/sum(Fc^2)) and apply it to - // every reflection. This mops up the residual radial Fobs/Fmodel mismatch that a single overall - // B leaves behind. The shells are fit on the work set only, so R-free stays untouched by the fit. - int n_shells = 0; - { - struct ShellPoint { double s2, fo, fc; }; - std::vector pts; - pts.reserve(fmodel.v.size()); - for (const auto &hv : fmodel.v) { - auto it = obs_by_hkl.find(hkl_key(hv.hkl)); - if (it == obs_by_hkl.end() || it->second.second) continue; // skip missing + free set - if (it->second.first <= 0) continue; - pts.push_back({ucell.calculate_1_d2(hv.hkl), it->second.first, std::abs(hv.value)}); + // --- observed amplitudes into the model ASU, keyed by hkl (also remember free flag) --- + // Observed amplitudes are the French-Wilson |F| already computed at the end of the merge + // (MergedReflection.F), so the model R-free / maps use exactly the same amplitudes as the + // written reflection file. + gemmi::AsuData> fobs; + fobs.unit_cell_ = ucell; + fobs.spacegroup_ = sg; + for (const MergedReflection &r : obs_in) { + if (std::isnan(r.F)) continue; + gemmi::Miller h{{r.h, r.k, r.l}}; + if (!asu.is_in(h)) h = asu.to_asu(h, gops).first; + fobs.v.push_back({h, {r.F, 1.0f}}); + out.obs_by_hkl[hkl_key(h)] = {r.F, r.rfree_flag}; } - std::sort(pts.begin(), pts.end(), [](const ShellPoint &a, const ShellPoint &b) { return a.s2 < b.s2; }); - int nb = std::max(6, std::min(40, static_cast(pts.size() / 300))); // ~300 refl/shell - std::vector shell_s2, shell_k; - for (int b = 0; b < nb; ++b) { - size_t i0 = pts.size() * b / nb, i1 = pts.size() * (b + 1) / nb; - double sum_fofc = 0, sum_fc2 = 0, sum_s2 = 0; - int n = 0; - for (size_t i = i0; i < i1; ++i) { - sum_fofc += pts[i].fo * pts[i].fc; - sum_fc2 += pts[i].fc * pts[i].fc; - sum_s2 += pts[i].s2; - ++n; + fobs.ensure_asu(); + fobs.ensure_sorted(); + + // --- scale Fmodel(+solvent) to Fobs: k_overall, anisotropic B, k_sol, b_sol --- + gemmi::Scaling scaling(ucell, sg); + scaling.use_solvent = true; + scaling.prepare_points(out.fmodel, fobs, &fmask); + scaling.fit_isotropic_b_approximately(); + scaling.fit_parameters(); + scaling.scale_data(out.fmodel, &fmask); // out.fmodel now holds the scaled, solvent-corrected Fmodel + out.k_sol = scaling.k_sol; + out.b_sol = scaling.b_sol; + out.k_overall = scaling.k_overall; + + // --- per-resolution-shell scaling on top of the overall scaling --- + // Fit a smooth scale K(1/d^2) least-squares per shell (sum(Fo*Fc)/sum(Fc^2)) and apply it to + // every reflection. This mops up the residual radial Fobs/Fmodel mismatch that a single overall + // B leaves behind. The shells are fit on the work set only, so R-free stays untouched by the fit. + { + struct ShellPoint { double s2, fo, fc; }; + std::vector pts; + pts.reserve(out.fmodel.v.size()); + for (const auto &hv : out.fmodel.v) { + auto it = out.obs_by_hkl.find(hkl_key(hv.hkl)); + if (it == out.obs_by_hkl.end() || it->second.second) continue; // skip missing + free set + if (it->second.first <= 0) continue; + pts.push_back({ucell.calculate_1_d2(hv.hkl), it->second.first, std::abs(hv.value)}); } - if (sum_fc2 > 0 && n > 0) { shell_s2.push_back(sum_s2 / n); shell_k.push_back(sum_fofc / sum_fc2); } - } - auto scale_at = [&](double s2) -> double { - if (shell_s2.empty()) return 1.0; - if (s2 <= shell_s2.front()) return shell_k.front(); - if (s2 >= shell_s2.back()) return shell_k.back(); - for (size_t i = 1; i < shell_s2.size(); ++i) - if (s2 <= shell_s2[i]) { - double t = (s2 - shell_s2[i - 1]) / (shell_s2[i] - shell_s2[i - 1]); - return shell_k[i - 1] * (1 - t) + shell_k[i] * t; + std::sort(pts.begin(), pts.end(), [](const ShellPoint &a, const ShellPoint &b) { return a.s2 < b.s2; }); + int nb = std::max(6, std::min(40, static_cast(pts.size() / 300))); // ~300 refl/shell + std::vector shell_s2, shell_k; + for (int b = 0; b < nb; ++b) { + size_t i0 = pts.size() * b / nb, i1 = pts.size() * (b + 1) / nb; + double sum_fofc = 0, sum_fc2 = 0, sum_s2 = 0; + int n = 0; + for (size_t i = i0; i < i1; ++i) { + sum_fofc += pts[i].fo * pts[i].fc; + sum_fc2 += pts[i].fc * pts[i].fc; + sum_s2 += pts[i].s2; + ++n; } - return shell_k.back(); - }; - for (auto &hv : fmodel.v) - hv.value *= static_cast(scale_at(ucell.calculate_1_d2(hv.hkl))); - n_shells = static_cast(shell_k.size()); - } + if (sum_fc2 > 0 && n > 0) { shell_s2.push_back(sum_s2 / n); shell_k.push_back(sum_fofc / sum_fc2); } + } + auto scale_at = [&](double s2) -> double { + if (shell_s2.empty()) return 1.0; + if (s2 <= shell_s2.front()) return shell_k.front(); + if (s2 >= shell_s2.back()) return shell_k.back(); + for (size_t i = 1; i < shell_s2.size(); ++i) + if (s2 <= shell_s2[i]) { + double t = (s2 - shell_s2[i - 1]) / (shell_s2[i] - shell_s2[i - 1]); + return shell_k[i - 1] * (1 - t) + shell_k[i] * t; + } + return shell_k.back(); + }; + for (auto &hv : out.fmodel.v) + hv.value *= static_cast(scale_at(ucell.calculate_1_d2(hv.hkl))); + out.n_shells = static_cast(shell_k.size()); + } - // --- R-work / R-free and map coefficients (2Fo-Fc and Fo-Fc, model phases) --- - gemmi::AsuData> map2fofc, mapfofc; - map2fofc.unit_cell_ = ucell; map2fofc.spacegroup_ = sg; - mapfofc.unit_cell_ = ucell; mapfofc.spacegroup_ = sg; - double num_w = 0, den_w = 0, num_f = 0, den_f = 0; - int n_w = 0, n_f = 0; - for (const auto &hv : fmodel.v) { - auto it = obs_by_hkl.find(hkl_key(hv.hkl)); - if (it == obs_by_hkl.end()) continue; - double Fo = it->second.first; - double Fc = std::abs(hv.value); - double phi = std::arg(hv.value); - if (it->second.second) { num_f += std::fabs(Fo - Fc); den_f += Fo; ++n_f; } - else { num_w += std::fabs(Fo - Fc); den_w += Fo; ++n_w; } - std::complex ph = std::polar(1.0f, static_cast(phi)); - map2fofc.v.push_back({hv.hkl, static_cast(2 * Fo - Fc) * ph}); - mapfofc.v.push_back({hv.hkl, static_cast(Fo - Fc) * ph}); - } + // --- R-work / R-free and map coefficients (2Fo-Fc and Fo-Fc, model phases) --- + out.map2fofc.unit_cell_ = ucell; out.map2fofc.spacegroup_ = sg; + out.mapfofc.unit_cell_ = ucell; out.mapfofc.spacegroup_ = sg; + double num_w = 0, den_w = 0, num_f = 0, den_f = 0; + for (const auto &hv : out.fmodel.v) { + auto it = out.obs_by_hkl.find(hkl_key(hv.hkl)); + if (it == out.obs_by_hkl.end()) continue; + double Fo = it->second.first; + double Fc = std::abs(hv.value); + double phi = std::arg(hv.value); + if (it->second.second) { num_f += std::fabs(Fo - Fc); den_f += Fo; ++out.n_f; } + else { num_w += std::fabs(Fo - Fc); den_w += Fo; ++out.n_w; } + std::complex ph = std::polar(1.0f, static_cast(phi)); + out.map2fofc.v.push_back({hv.hkl, static_cast(2 * Fo - Fc) * ph}); + out.mapfofc.v.push_back({hv.hkl, static_cast(Fo - Fc) * ph}); + } + out.r_work = den_w > 0 ? num_w / den_w : 1; + out.r_free = den_f > 0 ? num_f / den_f : 1; + return out; + }; - result.r_work = den_w > 0 ? num_w / den_w : 0; - result.r_free = den_f > 0 ? num_f / den_f : 0; - result.n_work = n_w; - result.n_free = n_f; - result.k_sol = scaling.k_sol; - result.b_sol = scaling.b_sol; - result.k_overall = scaling.k_overall; + // --- probe indexing (merohedral) ambiguities: keep the reindexing with the lowest R-free --- + // No-op for a holohedral crystal (no twin laws), e.g. lysozyme, where the only candidate is the + // identity. The enantiomorph/screw ambiguity is NOT probed here: |Fcalc| is the same for both + // hands, so R-free cannot distinguish them - that is resolved from the model hand above. + const auto reindex_ops = ReindexAmbiguityOperators(cell, sg->number); + Fit best = fit_model(obs); + bool did_reindex = false; + for (const auto &op : reindex_ops) { + Fit cand = fit_model(ReindexReflections(obs, op)); + if (cand.r_free < best.r_free) { best = std::move(cand); did_reindex = true; } + } + if (!reindex_ops.empty()) + logger.Info("Model validation: probed {} indexing solution(s); {} (R-free {:.4f})", + reindex_ops.size() + 1, + did_reindex ? "reindexed to the lower-R-free solution" : "kept the current indexing", + best.r_free); + + gemmi::AsuData> &fmodel = best.fmodel; + gemmi::AsuData> &map2fofc = best.map2fofc; + gemmi::AsuData> &mapfofc = best.mapfofc; + std::unordered_map> &obs_by_hkl = best.obs_by_hkl; + int n_shells = best.n_shells; + + result.r_work = best.r_work; + result.r_free = best.r_free; + result.n_work = best.n_w; + result.n_free = best.n_f; + result.k_sol = best.k_sol; + result.b_sol = best.b_sol; + result.k_overall = best.k_overall; // --- write the maps and score the 2Fo-Fc map at atom centres (a real map peaks there) --- const std::string p2 = output_prefix + "_2fofc.ccp4";