leonarski_fandClaude Opus 5
06a8a9d143
The d_min <= 0 return was the one failure path that left failure_reason empty. The report's MODEL_VALIDATION= NOT_PERFORMED line and the WARNING that goes with it are both gated on that string, so this case printed a blank reason and no warning - the "indistinguishable from a run that was never given --model" outcome the field was added to prevent. All four failure returns now set it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_016L1qig74oYQzfUJJZbbxFh
602 lines
31 KiB
C++
602 lines
31 KiB
C++
// 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;
|
|
}
|
|
|
|
// Cubic, not the default linear, for reading a map at a point. The maps are 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 anomalous sites.
|
|
constexpr int MAP_INTERPOLATION_ORDER = 3;
|
|
|
|
// How deep a trough at an atom has to be before the anomalous map is called inverted. Well clear of
|
|
// the couple of sigma a map with no anomalous signal reaches at its noisiest atom.
|
|
constexpr double ANOMALOUS_INVERSION_SIGMA = 5.0;
|
|
|
|
// 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), adopt the model's group as
|
|
// a LABEL and leave the reflections exactly where they are.
|
|
//
|
|
// It is tempting to reindex by the change-of-hand operator instead, and that is wrong. The two
|
|
// groups of an enantiomorphic pair have the same rotation operations - only their translations
|
|
// differ - so they transform hkl identically, share a reciprocal ASU, and split into Bijvoet
|
|
// hands identically: the label carries no handedness at all, and nothing about it needs undoing.
|
|
// What does carry the hand is the indexing the data already have, from the diffraction geometry,
|
|
// and with it the sign of every anomalous difference. The change-of-hand operator is the
|
|
// inversion, so reindexing by it swaps I(+) with I(-) - it does not correct the hand, it flips
|
|
// it, on the strength of a label the space-group search itself reports as undetermined. Where
|
|
// the model really is the wrong enantiomorph for this crystal, that flip does not reveal the
|
|
// disagreement but manufactures agreement. The anomalous difference map below is the only honest
|
|
// arbiter, and it is used to report the disagreement rather than to bury it.
|
|
const std::vector<MergedReflection> &obs = 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) {
|
|
result.adopted_model_enantiomorph = true;
|
|
logger.Info("Model validation: data space group {} is the enantiomorph of the model {}; "
|
|
"taking the model's group as the label, with no reindexing - the two groups "
|
|
"index identically, and reindexing would flip the anomalous differences",
|
|
dsg->short_name(), sg->hm);
|
|
}
|
|
}
|
|
}
|
|
|
|
// 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) {
|
|
result.failure_reason = "the merged reflections carry no resolution";
|
|
logger.Error("Model validation: {}", result.failure_reason);
|
|
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, MAP_INTERPOLATION_ORDER); ++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.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;
|
|
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});
|
|
}
|
|
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);
|
|
}
|
|
}
|
|
|
|
// --- 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);
|
|
}
|
|
if (result.anomalous_hands_disagree)
|
|
logger.Warning("Model validation: the anomalous density at the model's atoms is inverted "
|
|
"({} reads {:.1f} sigma, deeper than the highest peak): the data and the model "
|
|
"are in opposite hands. The reflections have NOT been reindexed to make them "
|
|
"agree - either the model is the wrong enantiomorph for this crystal, or the "
|
|
"data were indexed in the wrong hand, and reindexing would hide which",
|
|
result.anomalous_deepest_site, result.anomalous_deepest_sigma);
|
|
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;
|
|
|
|
// Adopting the model's enantiomorph is a change of the space-group LABEL and nothing else. The
|
|
// two groups have the same rotation operations, so the same reflections, indexed the way they
|
|
// already are, are as good a description of one group as of the other; what the file gains is a
|
|
// group that agrees with the model it will be refined against. Reindexing here would swap the
|
|
// Bijvoet mates and so change the data - see the note in ValidateAgainstModel.
|
|
if (validation.adopted_model_enantiomorph && validation.model_space_group_number > 0) {
|
|
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, {} ({}), "
|
|
"as a label - no reflection moved", sg ? sg->short_name() : "?", space_group_number);
|
|
}
|
|
|
|
// 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;
|
|
}
|