mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 67 Wiki Activity
Files
d6924aa90130234472ee4ef21708b5fd6522b9b0
Jungfraujoch/rugnux/ModelValidation.cpp
T
leonarski_fandClaude Opus 4.8 d6924aa901 rugnux: add --model for R-free and 2Fo-Fc / Fo-Fc maps from an atomic model 
New offline capability: `rugnux --model model.pdb` validates the merged intensities
against a PDB model. It computes structure factors from the model on a grid (IT92 x-ray
form factors, FFT), adds an optimized flat bulk solvent (SolventMasker + gemmi Scaling:
k_sol, b_sol, overall k and anisotropic B) and a per-resolution-shell scale, then reports
R-work / R-free and writes 2Fo-Fc and Fo-Fc CCP4 maps plus a map-coefficient MTZ
(FP, FC, PHIC, FWT/PHWT, DELFWT/PHDELWT, FREE). No refinement of the structure itself;
the model is only re-fractionalized into the data cell (a rigid cell adjustment).

The observed amplitudes are the French-Wilson |F| already on the merged reflections
(MergedReflection.F, computed at end of merge), so the R-free/maps use exactly the same
amplitudes as the written reflection file.

ModelValidation.{h,cpp} holds the logic; wired via --model (ProcessConfig.model_path) into
both the full pipeline (Rugnux::Run write site) and the --scale re-merge path. Links the
unified `gemmi` static target for the model / structure-factor / map headers.

Verified on lyso_ref: R-free 0.346, mean 2Fo-Fc at atom centres +2.16 sigma, maps written.

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

293 lines
13 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"
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) {
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;
}
// Resolution limit from the data (the merged set is already resolution-trimmed).
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 (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>> fmodel =
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);
// --- observed amplitudes into the model ASU, keyed by hkl (also remember free flag) ---
gemmi::GroupOps gops = sg->operations();
gemmi::ReciprocalAsu asu(sg);
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
fobs.unit_cell_ = ucell;
fobs.spacegroup_ = sg;
std::unordered_map<long, std::pair<double, bool>> 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 : merged) {
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();
// --- 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(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<ShellPoint> 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)});
}
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 : fmodel.v)
hv.value *= static_cast<float>(scale_at(ucell.calculate_1_d2(hv.hkl)));
n_shells = static_cast<int>(shell_k.size());
}
// --- R-work / R-free and map coefficients (2Fo-Fc and Fo-Fc, model phases) ---
gemmi::AsuData<std::complex<float>> 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<float> ph = std::polar(1.0f, static_cast<float>(phi));
map2fofc.v.push_back({hv.hkl, static_cast<float>(2 * Fo - Fc) * ph});
mapfofc.v.push_back({hv.hkl, static_cast<float>(Fo - Fc) * ph});
}
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;
// --- 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