mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 67 Wiki Activity
Files
386f345758a6e80652e2c19618ac39148b6ff2dc
Jungfraujoch/rugnux/ModelValidation.cpp
T
leonarski_fandClaude Opus 4.8 386f345758 rugnux: probe indexing ambiguities against the model in validation 
ModelValidation now fits the scaled model to each reindexing of the observed
reflections (identity plus the twin laws from ReindexAmbiguityOperators) and
keeps the one with the lowest R-free. The scale + per-shell + R computation is
factored into a per-candidate lambda so every candidate is scored identically.

No-op for a holohedral crystal (no twin laws), e.g. lysozyme, where the only
candidate is the identity. The enantiomorph/screw ambiguity is still resolved
from the model hand, not probed - |Fcalc| is the same for both hands so R-free
cannot distinguish them.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
2026-07-13 08:20:27 +02:00

357 lines
17 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "ModelValidation.h"
#include <cmath>
#include <complex>
#include <array>
#include <algorithm>
#include <vector>
#include <unordered_map>
#include <gemmi/pdb.hpp> // read_pdb
#include <gemmi/gz.hpp> // MaybeGzipped
#include <gemmi/it92.hpp> // IT92 x-ray form factors
#include <gemmi/dencalc.hpp> // DensityCalculator
#include <gemmi/fourier.hpp> // transform_map_to_f_phi, get_f_phi_on_grid, transform_f_phi_grid_to_map
#include <gemmi/solmask.hpp> // SolventMasker
#include <gemmi/scaling.hpp> // Scaling (bulk solvent + anisotropic B)
#include <gemmi/ccp4.hpp> // Ccp4 map I/O
#include <gemmi/mtz.hpp> // Mtz (map-coefficient output)
#include "../common/Logger.h"
#include "../image_analysis/scale_merge/ReindexAmbiguity.h" // ReindexReflections
namespace {
using Table = gemmi::IT92<float>;
// Stable key for a Miller index reduced into the ASU (indices are small, well within +/-512).
long hkl_key(const gemmi::Miller &h) {
return (h[0] + 512L) * 1048576 + (h[1] + 512L) * 1024 + (h[2] + 512L);
}
// Write an FFT of ASU map coefficients as a CCP4 map; return its RMS (for reporting / sigma units).
double write_ccp4(gemmi::AsuData<std::complex<float>> &coef, const std::string &path) {
coef.ensure_sorted();
std::array<int, 3> size = gemmi::get_size_for_hkl(coef, {{0, 0, 0}}, 3.0);
gemmi::Grid<float> map =
gemmi::transform_f_phi_grid_to_map(gemmi::get_f_phi_on_grid<float>(coef, size, true));
gemmi::Ccp4<float> ccp4;
ccp4.grid = map;
ccp4.update_ccp4_header(2);
ccp4.write_ccp4_map(path);
return ccp4.hstats.rms;
}
} // namespace
ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &merged,
const UnitCell &cell,
const std::string &model_path,
const std::string &output_prefix,
Logger &logger,
std::optional<int> data_space_group_number) {
ModelValidationResult result;
// --- read the atomic model ---
gemmi::Structure st;
try {
st = gemmi::read_pdb(gemmi::MaybeGzipped(model_path));
} catch (const std::exception &e) {
logger.Error("Model validation: cannot read model {}: {}", model_path, e.what());
return result;
}
if (st.models.empty() || !st.cell.is_crystal()) {
logger.Error("Model validation: model {} has no atoms or no unit cell", model_path);
return result;
}
const gemmi::SpaceGroup *sg = st.find_spacegroup();
if (!sg) {
logger.Error("Model validation: model {} has no usable space group", model_path);
return result;
}
// If the data was indexed in the enantiomorph of the model's space group (e.g. data P4(1)2(1)2,
// model P4(3)2(1)2 - the merged intensities cannot tell them apart), reindex the observed
// reflections into the model's hand so the two settings agree. This does not change the
// R-factors (which use |F|), but keeps the observed data consistent with the model.
std::vector<MergedReflection> reindexed;
const std::vector<MergedReflection> *obs_ptr = &merged;
if (data_space_group_number && *data_space_group_number != sg->number) {
const gemmi::SpaceGroup *dsg = gemmi::find_spacegroup_by_number(*data_space_group_number);
if (dsg && dsg->is_enantiomorphic() && sg->is_enantiomorphic()) {
gemmi::GroupOps eops = dsg->operations();
eops.change_basis_forward(dsg->change_of_hand_op());
const gemmi::SpaceGroup *enant = gemmi::find_spacegroup_by_ops(eops);
if (enant && enant->number == sg->number) {
reindexed = ReindexReflections(merged, dsg->change_of_hand_op());
obs_ptr = &reindexed;
logger.Info("Model validation: data space group {} is the enantiomorph of the model {}; "
"reindexed the observed reflections into the model's hand",
dsg->short_name(), sg->hm);
}
}
}
const std::vector<MergedReflection> &obs = *obs_ptr;
// Resolution limit from the data (the merged set is already resolution-trimmed).
double d_min = 0.0;
for (const MergedReflection &r : obs)
if (r.d > 0 && (d_min == 0.0 || r.d < d_min))
d_min = r.d;
if (d_min <= 0.0) {
logger.Error("Model validation: merged reflections carry no resolution");
return result;
}
// Re-fractionalize the model into the data cell (rigid cell adjustment; no refinement).
const gemmi::UnitCell data_cell = cell; // UnitCell -> gemmi::UnitCell
if (data_cell.is_crystal()) {
gemmi::UnitCell old = st.cell;
for (gemmi::Model &m : st.models)
for (gemmi::Chain &ch : m.chains)
for (gemmi::Residue &r : ch.residues)
for (gemmi::Atom &a : r.atoms)
a.pos = data_cell.orthogonalize(old.fractionalize(a.pos));
st.cell = data_cell;
}
st.setup_cell_images();
const gemmi::UnitCell &ucell = st.cell;
logger.Info("Model validation: {} atoms, cell a={:.2f} b={:.2f} c={:.2f}, sg {}, to {:.2f} A",
gemmi::count_atom_sites(st.models[0]), ucell.a, ucell.b, ucell.c, sg->hm, d_min);
// --- Fcalc (atomic) via electron density on a grid + FFT ---
gemmi::DensityCalculator<Table, float> 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]);
dc.put_model_density_on_grid(st.models[0]);
gemmi::AsuData<std::complex<float>> 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 ---
// Refmac radii give a slightly lower R than the Cctbx set on our test cases, at no cost.
gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac);
gemmi::Grid<float> 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);
masker.put_mask_on_grid(mask_grid, st.models[0]);
gemmi::AsuData<std::complex<float>> fmask =
gemmi::transform_map_to_f_phi(mask_grid, true).prepare_asu_data(dc.d_min, 0);
gemmi::GroupOps gops = sg->operations();
gemmi::ReciprocalAsu asu(sg);
// --- 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<std::complex<float>> fmodel, map2fofc, mapfofc;
std::unordered_map<long, std::pair<double, bool>> 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<MergedReflection> &obs_in) -> Fit {
Fit out;
out.fmodel = fcalc; // copy the atomic structure factors; scaling mutates them in place
// --- 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<gemmi::ValueSigma<float>> 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};
}
fobs.ensure_asu();
fobs.ensure_sorted();
// --- scale Fmodel(+solvent) to Fobs: k_overall, anisotropic B, k_sol, b_sol ---
gemmi::Scaling<float> 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<ShellPoint> 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)});
}
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<int>(pts.size() / 300))); // ~300 refl/shell
std::vector<double> 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;
}
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<float>(scale_at(ucell.calculate_1_d2(hv.hkl)));
out.n_shells = static_cast<int>(shell_k.size());
}
// --- 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<float> ph = std::polar(1.0f, static_cast<float>(phi));
out.map2fofc.v.push_back({hv.hkl, static_cast<float>(2 * Fo - Fc) * ph});
out.mapfofc.v.push_back({hv.hkl, static_cast<float>(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;
};
// --- 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<std::complex<float>> &fmodel = best.fmodel;
gemmi::AsuData<std::complex<float>> &map2fofc = best.map2fofc;
gemmi::AsuData<std::complex<float>> &mapfofc = best.mapfofc;
std::unordered_map<long, std::pair<double, bool>> &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";
const std::string pd = output_prefix + "_fofc.ccp4";
double rms2 = write_ccp4(map2fofc, p2);
write_ccp4(mapfofc, pd);
// Recompute the 2Fo-Fc real-space grid once more to sample it at atom positions.
{
map2fofc.ensure_sorted();
std::array<int, 3> size = gemmi::get_size_for_hkl(map2fofc, {{0, 0, 0}}, 3.0);
gemmi::Grid<float> grid =
gemmi::transform_f_phi_grid_to_map(gemmi::get_f_phi_on_grid<float>(map2fofc, size, true));
double s = 0; int n = 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) { s += grid.interpolate_value(a.pos); ++n; }
result.mean_atom_density_sigma = (n > 0 && rms2 > 0) ? (s / n) / rms2 : 0;
}
// --- MTZ of map coefficients so the maps can be re-opened / rebuilt in Coot etc. ---
try {
gemmi::Mtz mtz(true);
mtz.spacegroup = sg;
mtz.set_cell_for_all(ucell);
mtz.add_dataset("model_validation");
mtz.add_column("FP", 'F', -1, -1, false);
mtz.add_column("FC", 'F', -1, -1, false);
mtz.add_column("PHIC", 'P', -1, -1, false);
mtz.add_column("FWT", 'F', -1, -1, false);
mtz.add_column("PHWT", 'P', -1, -1, false);
mtz.add_column("DELFWT", 'F', -1, -1, false);
mtz.add_column("PHDELWT", 'P', -1, -1, false);
mtz.add_column("FREE", 'I', -1, -1, false);
std::vector<float> data;
int nref = 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_deg = gemmi::phase_in_angles(hv.value);
data.insert(data.end(), {static_cast<float>(hv.hkl[0]), static_cast<float>(hv.hkl[1]),
static_cast<float>(hv.hkl[2]),
static_cast<float>(Fo), static_cast<float>(Fc),
static_cast<float>(phi_deg),
static_cast<float>(2 * Fo - Fc), static_cast<float>(phi_deg),
static_cast<float>(Fo - Fc), static_cast<float>(phi_deg),
it->second.second ? 0.0f : 1.0f});
++nref;
}
mtz.nreflections = nref;
mtz.data = std::move(data);
mtz.write_to_file(output_prefix + "_maps.mtz");
} catch (const std::exception &e) {
logger.Warning("Model validation: could not write map MTZ: {}", e.what());
}
result.ok = true;
result.maps_prefix = output_prefix;
logger.Info("Model validation: R-work={:.4f} ({} refl) R-free={:.4f} ({} refl) [overall + {}-shell scaling]",
result.r_work, result.n_work, result.r_free, result.n_free, n_shells);
logger.Info("Model validation: bulk solvent k_sol={:.3f} b_sol={:.1f}, k_overall={:.3f}",
result.k_sol, result.b_sol, result.k_overall);
logger.Info("Model validation: mean 2Fo-Fc density at atom centres = {:.2f} sigma", result.mean_atom_density_sigma);
logger.Info("Model validation: wrote {}_2fofc.ccp4, {}_fofc.ccp4, {}_maps.mtz",
output_prefix, output_prefix, output_prefix);
return result;
}
Reference in New Issue View Git Blame Copy Permalink