mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases 72 Wiki Activity
Files
b0234a1b4b402594faaff91081b71ae3f2d4164a
Jungfraujoch/rugnux/ModelValidation.cpp
T
leonarski_fandClaude Opus 5 b0234a1b4b rugnux: --model names the strongest anomalous scatterers, not just R-free 
Where the merge kept the Bijvoet split - which the default Friedel-averaged
merge does, so this needs no -A - --model now computes an anomalous difference
Fourier, coefficients F(+)-F(-) on the model phase turned back by 90 degrees
over the acentric reflections, writes it as <prefix>_anom.ccp4, and reads it at
the model's own atom centres. The ten highest are logged and land in the results
report as ANOMALOUS_SITE_01..10, each named by the atom, residue and chain it
sits on, with its height in map sigma. Reading at the sites rather than
searching the map for blobs is what ANODE does, and it says what carries the
signal instead of leaving coordinates to look up; a scatterer the model does not
contain is by construction absent from the list, which is what the map file is
for.

The list is always ten long, so its information is in the heights. On an
S-SAD dataset at 5 keV the six cysteine sulfurs take ranks 1-6 (10.1 down to
4.2 sigma) and rank 7 falls to 2.2; the same crystal at 13 keV, where sulfur has
no anomalous signal to speak of, tops out at 2.7 sigma with a water and a main
chain carbon inside the top five. So the same list reads as a signal / no-signal
gate, not only as a ranking.

The reading is cubic. Grid::interpolate_value defaults to trilinear, and on the
d_min/3 map grid that under-reads a peak this sharp by up to a quarter, unevenly
enough to reorder the sites; cubic on the same grid is within about 2% of a
sample-rate-10 one. The map r.m.s. is unchanged by the sample rate, so the
normalisation was never at fault. The 2Fo-Fc reading behind
MEAN_ATOM_DENSITY_SIGMA is deliberately left on the default so a number already
in existing reports does not move.

An anomalous merge keeps each Bijvoet mate as a row of its own and both rows
carry the same F(+)/F(-) pair, attached to the + index of the Friedel ASU of the
frame the merge was made in, so the differences are read from the merged
reflections as they came in - carried into the model's frame here, on the +
rows only. Taking the - rows as well would give one reflection both signs of its
difference and let the last row written decide.

write_ccp4 splits into an FFT and a write so the anomalous map can be both
written and sampled from one transform; the 2Fo-Fc atom-centre score, which used
to repeat that FFT, now reuses the grid.

Cross-checked against ANODE run on this output: identical identification - its
three peaks are the three disulfides, each maximum 0.64-0.70 A off a sulfur we
name - with our heights about 1.35x lower. Sharpening, grid and sigma
normalisation are excluded; the residual is SHELXC's FA preparation, which fed
ANODE 2166 reflections against our 2534.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_014zTy4Bpi4pPHw4bybf7q2R
2026-08-27 14:51:06 +02:00

569 lines
29 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 <algorithm>
#include <cmath>
#include <complex>
#include <array>
#include <string>
#include <vector>
#include <unordered_map>
#include <gemmi/mmread_gz.hpp> // read_structure_gz
#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);
}
// FFT ASU map coefficients into a real-space map.
gemmi::Grid<float> map_from_coefficients(gemmi::AsuData<std::complex<float>> &coef) {
coef.ensure_sorted();
std::array<int, 3> size = gemmi::get_size_for_hkl(coef, {{0, 0, 0}}, 3.0);
return gemmi::transform_f_phi_grid_to_map(gemmi::get_f_phi_on_grid<float>(coef, size, true));
}
// Write a map as CCP4; return its RMS (the sigma the map is read in).
double write_ccp4(const gemmi::Grid<float> &map, const std::string &path) {
gemmi::Ccp4<float> ccp4;
ccp4.grid = map;
ccp4.update_ccp4_header(2);
ccp4.write_ccp4_map(path);
return ccp4.hstats.rms;
}
// How many anomalous sites the report names. The strongest few are what says whether the anomalous
// signal is there and what carries it; a full site list is what the map file is for.
constexpr size_t MAX_ANOMALOUS_SITES = 10;
} // 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,
bool probe_indexing_ambiguity) {
ModelValidationResult result;
result.model_path = model_path;
// --- read the atomic model ---
gemmi::Structure st;
try {
// Detect, not the default: without it GEMMI picks the format from the extension and only
// falls back to the content when it does not recognise one. A model arrives named however
// whoever produced it named it, so the file itself is the better authority.
st = gemmi::read_structure_gz(model_path, gemmi::CoorFormat::Detect);
} catch (const std::exception &e) {
result.failure_reason = fmt::format("cannot read model {}: {}", model_path, e.what());
logger.Error("Model validation: {}", result.failure_reason);
return result;
}
if (st.models.empty() || !st.cell.is_crystal()) {
result.failure_reason = fmt::format("model {} has no atoms or no unit cell", model_path);
logger.Error("Model validation: {}", result.failure_reason);
return result;
}
const gemmi::SpaceGroup *sg = st.find_spacegroup();
if (!sg) {
result.failure_reason = fmt::format("model {} has no usable space group", model_path);
logger.Error("Model validation: {}", result.failure_reason);
return result;
}
result.model_space_group_number = sg->number;
// 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;
result.hand_op = dsg->change_of_hand_op();
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;
};
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;
// The model is scaled to the data with an overall scale, an anisotropic B and a flat bulk
// solvent only - the standard, few-parameter model that refinement programs use. A dataset-
// specific free-form per-resolution-shell rescale would lower this dataset's R a little, but
// it reshapes each map's radial amplitude profile differently, so a batch of maps would no
// longer be directly comparable. For a fragment-screening / PanDDA campaign, comparable maps
// across datasets matter more than the last bit of per-dataset R, so it is deliberately omitted.
// --- 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;
};
// --- indexing (merohedral) ambiguity ---
// When a reference MTZ was supplied, the data were already reindexed to agree with the reference
// intensities (at the merge stage for rotation data, per image in stills scaling), and that
// choice is authoritative - we keep it. Only with a model and NO reference do we resolve the
// ambiguity here, as a fallback, by fitting each candidate reindexing and keeping the lowest
// R-free. A no-op either way for a holohedral crystal (no twin laws). The
// enantiomorph/screw ambiguity is never probed by R-free: |Fcalc| is the same for both hands, so
// it cannot distinguish them - that is taken from the model hand above.
Fit best = fit_model(obs);
if (probe_indexing_ambiguity) {
const auto reindex_ops = ReindexAmbiguityOperators(cell, sg->number);
std::vector<double> candidate_r_free{best.r_free}; // identity first, then the twin laws
for (const auto &op : reindex_ops) {
Fit cand = fit_model(ReindexReflections(obs, op));
candidate_r_free.push_back(cand.r_free);
if (cand.r_free < best.r_free) { best = std::move(cand); result.indexing_op = op; }
}
if (!reindex_ops.empty()) {
// The runner-up as well as the winner: the margin between them is what says whether the
// choice was made on evidence, and on weak data the two can come out within noise.
std::sort(candidate_r_free.begin(), candidate_r_free.end());
logger.Info("Model validation: no reference - probed {} indexing solution(s) against the model; "
"{} (R-free {:.4f}, runner-up {:.4f})", candidate_r_free.size(),
result.indexing_op == gemmi::Op::identity() ? "kept the current indexing"
: "reindexed to the lower-R-free solution",
candidate_r_free[0], candidate_r_free[1]);
}
}
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;
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 gemmi::Grid<float> grid2fofc = map_from_coefficients(map2fofc);
const double rms2 = write_ccp4(grid2fofc, output_prefix + "_2fofc.ccp4");
write_ccp4(map_from_coefficients(mapfofc), output_prefix + "_fofc.ccp4");
{
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 += grid2fofc.interpolate_value(a.pos); ++n; }
result.mean_atom_density_sigma = (n > 0 && rms2 > 0) ? (s / n) / rms2 : 0;
}
// --- anomalous difference map, where the merge kept the Bijvoet split ---
// Coefficients F(+) - F(-) carried on the model phase turned back by 90 degrees. Its peaks sit on
// the anomalous scatterers, so reading the map at each of the model's own atoms names them,
// rather than leaving a list of coordinates for someone to look up.
// Following ANODE, Thorn & Sheldrick (2011) J. Appl. Cryst. 44, 1285-1287
{
// 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_number
? gemmi::find_spacegroup_by_number(*data_space_group_number) : nullptr;
if (data_sg == nullptr)
data_sg = sg;
const gemmi::ReciprocalAsu data_asu(data_sg);
const gemmi::GroupOps data_gops = data_sg->operations();
std::unordered_map<long, float> 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.hand_op == gemmi::Op::identity()))
h = result.hand_op.apply_to_hkl(h);
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<std::complex<float>> 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<float>(std::arg(hv.value) - M_PI / 2);
mapanom.v.push_back({hv.hkl, it->second * std::polar(1.0f, phi)});
}
result.anomalous_pairs = static_cast<int>(mapanom.v.size());
if (!mapanom.v.empty()) {
const gemmi::Grid<float> grid = map_from_coefficients(mapanom);
const double rms = write_ccp4(grid, output_prefix + "_anom.ccp4");
std::vector<ModelValidationResult::AnomalousSite> sites;
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;
// Cubic, not the default linear: the map is sampled every d_min/3, and a
// peak that sharp read by trilinear interpolation comes out up to a quarter
// low - unevenly enough to reorder the sites.
sites.push_back({fmt::format("{} {} {}{}", a.name, r.name, ch.name,
r.seqid.str()),
rms > 0 ? grid.interpolate_value(a.pos, 3) / rms : 0.0});
}
std::sort(sites.begin(), sites.end(),
[](const auto &x, const auto &y) { return x.sigma > y.sigma; });
if (sites.size() > MAX_ANOMALOUS_SITES)
sites.resize(MAX_ANOMALOUS_SITES);
result.anomalous_sites = std::move(sites);
}
}
// --- 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 + anisotropic B + bulk solvent]",
result.r_work, result.n_work, result.r_free, result.n_free);
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);
if (!result.anomalous_sites.empty()) {
std::string sites;
for (const auto &s : result.anomalous_sites)
sites += fmt::format("{}{} {:.1f}", sites.empty() ? "" : ", ", s.label, s.sigma);
logger.Info("Model validation: anomalous difference map from {} Bijvoet pairs; strongest "
"density at the model's atoms (sigma): {}", result.anomalous_pairs, sites);
}
logger.Info("Model validation: wrote {}_2fofc.ccp4, {}_fofc.ccp4{}, {}_maps.mtz",
output_prefix, output_prefix,
result.anomalous_sites.empty() ? "" : ", " + output_prefix + "_anom.ccp4",
output_prefix);
return result;
}
namespace {
// The reindexing operator as it reads on Miller indices ("k,h,-l" rather than "y,x,-z").
std::string hkl_triplet(const gemmi::Op &op) {
std::string t = op.triplet();
std::replace(t.begin(), t.end(), 'x', 'h');
std::replace(t.begin(), t.end(), 'y', 'k');
std::replace(t.begin(), t.end(), 'z', 'l');
return t;
}
} // namespace
int AdoptModelFrame(const ModelValidationResult &validation,
std::vector<MergedReflection> &merged,
int data_space_group_number,
bool merge_friedel,
Logger &logger) {
int space_group_number = data_space_group_number;
if (!validation.ok)
return space_group_number;
// The change of hand also changes the space group the reflections are in: P4(1)2(1)2 data written
// in the model's hand are P4(3)2(1)2 data. The merged intensities are the same measurements either
// way - the enantiomorph is exactly what they cannot distinguish - so the model's hand is the only
// evidence there is, and adopting it is what makes the file refinable against that model.
if (!(validation.hand_op == gemmi::Op::identity()) && validation.model_space_group_number > 0) {
merged = ReindexMergedIntoAsu(merged, validation.hand_op,
validation.model_space_group_number, merge_friedel);
space_group_number = validation.model_space_group_number;
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number);
logger.Info("Model validation: the written reflections take the model's enantiomorph, {} ({}), "
"reindexed by {}", sg ? sg->short_name() : "?", space_group_number,
hkl_triplet(validation.hand_op));
}
// The alternative indexing, by contrast, is metric- and group-preserving: only the labels move.
if (!(validation.indexing_op == gemmi::Op::identity())) {
merged = ReindexMergedIntoAsu(merged, validation.indexing_op, space_group_number, merge_friedel);
logger.Info("Model validation: the written reflections take the model's indexing, reindexed by {}",
hkl_triplet(validation.indexing_op));
}
return space_group_number;
}
std::vector<MergedReflection> ModelReferenceIntensities(const std::string &model_path,
const std::optional<UnitCell> &cell,
std::optional<int> space_group_number,
double d_min,
Logger &logger) {
std::vector<MergedReflection> out;
if (!(d_min > 0.0)) {
logger.Warning("Model reference: no resolution limit to compute the model intensities to");
return out;
}
gemmi::Structure st;
try {
// Detect, not the default: without it GEMMI picks the format from the extension and only
// falls back to the content when it does not recognise one. A model arrives named however
// whoever produced it named it, so the file itself is the better authority.
st = gemmi::read_structure_gz(model_path, gemmi::CoorFormat::Detect);
} catch (const std::exception &e) {
logger.Error("Model reference: cannot read model {}: {}", model_path, e.what());
return out;
}
if (st.models.empty() || !st.cell.is_crystal()) {
logger.Error("Model reference: model {} has no atoms or no unit cell", model_path);
return out;
}
// Put the model in the cell and group the run works in, where it knows them, so the reference is
// indexed the way the data are. The correlation that consumes this matches on hkl, so a small cell
// difference costs nothing; the space group is what has to agree.
if (cell.has_value()) {
const gemmi::UnitCell target = *cell;
if (target.is_crystal()) {
const 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 = target.orthogonalize(old.fractionalize(a.pos));
st.cell = target;
}
}
if (space_group_number.has_value())
if (const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(*space_group_number))
st.spacegroup_hm = sg->xhm();
const gemmi::SpaceGroup *sg = st.find_spacegroup();
if (!sg) {
logger.Error("Model reference: model {} has no usable space group", model_path);
return out;
}
st.setup_cell_images();
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 at the standard constants. Nothing here is fitted - there are no observations
// yet - but without it the few lowest-resolution reflections are the largest and the most wrong,
// and a correlation on raw intensities would be led by them.
constexpr double K_SOL = 0.35;
constexpr double B_SOL = 46.0;
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);
std::unordered_map<long, std::complex<float>> mask_by_hkl;
mask_by_hkl.reserve(fmask.v.size());
for (const auto &hv : fmask.v)
mask_by_hkl[hkl_key(hv.hkl)] = hv.value;
const gemmi::UnitCell &ucell = st.cell;
out.reserve(fcalc.v.size());
for (const auto &hv : fcalc.v) {
const double d = ucell.calculate_d(hv.hkl);
if (!(d > 0.0))
continue;
std::complex<float> f = hv.value;
const auto it = mask_by_hkl.find(hkl_key(hv.hkl));
if (it != mask_by_hkl.end())
f += static_cast<float>(K_SOL * std::exp(-B_SOL / (4.0 * d * d))) * it->second;
const double F = std::abs(f);
out.push_back(MergedReflection{.h = hv.hkl[0], .k = hv.hkl[1], .l = hv.hkl[2],
.I = static_cast<float>(F * F),
.d = static_cast<float>(d)});
}
logger.Info("Model reference: {} intensities computed from {} to {:.2f} A, space group {}",
out.size(), model_path, d_min, sg->short_name());
return out;
}
Reference in New Issue View Git Blame Copy Permalink