leonarski_fandClaude Opus 5
0903a8eee3
Two things wrong with <prefix>_maps.mtz beside <prefix>.mtz. It took the space group from the MODEL, unconditionally. That was invisible while a model's hand was always adopted; now that a model decides nothing until it beats its own null, a rejected model that asserts the other enantiomorph leaves the two files disagreeing - the reflections in the group the data were merged in, the map coefficients labelled with the model's. An enantiomorphic pair indexes identically, so the coefficients are the same numbers either way and only the label moves, but the two groups have different screw translations: a reader that expands symmetry out of the map file was doing it in the wrong group. It now takes the group the same way AdoptModelFrame does a few lines below - the model's only where the hand was adopted, the data's otherwise. And it carried no wavelength at all, DWAVEL 0.00000 on both datasets, because nothing in ValidateAgainstModel knew it; it is now passed in from the experiment the two callers already hold. Measured on one rotation data set: with the model accepted both files read P 43 21 2, with an unrelated model rejected both read P 41 21 2, and both carry the collection wavelength. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
959 lines
53 KiB
C++
959 lines
53 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 <numeric>
|
|
#include <random>
|
|
#include <string>
|
|
#include <utility>
|
|
#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/JFJochMath.h" // PI (M_PI is not standard, and MSVC does not define it)
|
|
#include "../common/Logger.h"
|
|
#include "../common/ParallelFor.h" // ParallelFor
|
|
#include "../image_analysis/scale_merge/ReindexAmbiguity.h" // ReindexReflections
|
|
#include "RigidBodyRefine.h"
|
|
#include "SigmaA.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);
|
|
}
|
|
|
|
// 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;
|
|
|
|
// How many random placements of the same model the real fit is compared against. The verdict is
|
|
// (real - mean)/sd of this sample, so what matters is not the mean but how well the SPREAD is
|
|
// pinned: the relative error on an sd from n draws is 1/sqrt(2(n-1)), and a sample that happens to
|
|
// come out narrow is what turns a model that does not fit into one that appears to. Measured on the
|
|
// case that sits closest to the gate - a model of an unrelated protein, 1.83 sigma against a
|
|
// threshold of 3 - the chance of it reading over the gate on a different seed is 31% at n=3, 17% at
|
|
// n=5, 6% at n=9.
|
|
//
|
|
// Nine rather than five because the replicates run concurrently: the null costs the SLOWEST of them
|
|
// rather than their sum, and the slowest of nine is barely above the slowest of five, so the extra
|
|
// four are close to free in wall clock (measured 2.5 s against 2.6 s) up to the thread count.
|
|
constexpr int NULL_REPLICATES = 9;
|
|
|
|
// Fixed, so the same data and the same model give the same verdict on every run.
|
|
constexpr unsigned NULL_SEED = 20260902;
|
|
|
|
// How far above its own null a fit has to sit before the model is allowed to decide anything. The cut
|
|
// is in sigma of that null and not in R: measured, the R a model that explains nothing reaches moves
|
|
// with the model's atom count and B-factors as much as with the data, so no value of R separates the
|
|
// two on its own. A model that belongs to the crystal measures around +4 sigma here and an unrelated
|
|
// one under +1.
|
|
constexpr double MODEL_FIT_SIGMA = 3.0;
|
|
|
|
// Mean and sample standard deviation of a small sample.
|
|
std::pair<double, double> mean_sd(const std::vector<double> &v) {
|
|
if (v.size() < 2)
|
|
return {v.empty() ? 0.0 : v.front(), 0.0};
|
|
const double mean = std::accumulate(v.begin(), v.end(), 0.0) / static_cast<double>(v.size());
|
|
double s2 = 0;
|
|
for (double x : v)
|
|
s2 += (x - mean) * (x - mean);
|
|
return {mean, std::sqrt(s2 / static_cast<double>(v.size() - 1))};
|
|
}
|
|
|
|
// A rotation drawn uniformly from SO(3), through a uniform random unit quaternion.
|
|
// Following Shoemake (1992) Graphics Gems III, 124-132
|
|
gemmi::Mat33 random_rotation(std::mt19937 &rng) {
|
|
std::uniform_real_distribution<double> u(0.0, 1.0);
|
|
const double u1 = u(rng), t2 = 2 * PI * u(rng), t3 = 2 * PI * u(rng);
|
|
const double r1 = std::sqrt(1 - u1), r2 = std::sqrt(u1);
|
|
const double x = r1 * std::sin(t2), y = r1 * std::cos(t2), z = r2 * std::sin(t3), w = r2 * std::cos(t3);
|
|
return {1 - 2 * (y * y + z * z), 2 * (x * y - z * w), 2 * (x * z + y * w),
|
|
2 * (x * y + z * w), 1 - 2 * (x * x + z * z), 2 * (y * z - x * w),
|
|
2 * (x * z - y * w), 2 * (y * z + x * w), 1 - 2 * (x * x + y * y)};
|
|
}
|
|
|
|
// Everything one fit moves: the model itself, and the structure factors that follow it. The real
|
|
// model has one of these and every null replicate gets a copy of its own, which is what lets the
|
|
// replicates run at the same time without stepping on each other.
|
|
struct ModelState {
|
|
gemmi::Structure st;
|
|
gemmi::AsuData<std::complex<float>> fcalc, fmask;
|
|
};
|
|
|
|
// Turn the model about its own centroid, so it keeps its place in the cell and loses its orientation.
|
|
void reorient_about_centroid(gemmi::Model &model, const gemmi::Mat33 &rot) {
|
|
std::vector<gemmi::Position> pos = ModelPositions(model);
|
|
if (pos.empty())
|
|
return;
|
|
gemmi::Vec3 centre;
|
|
for (const gemmi::Position &p : pos)
|
|
centre += p;
|
|
centre *= 1.0 / static_cast<double>(pos.size());
|
|
for (gemmi::Position &p : pos)
|
|
p = gemmi::Position(rot.multiply(gemmi::Vec3(p) - centre) + centre);
|
|
SetModelPositions(model, pos);
|
|
}
|
|
|
|
} // namespace
|
|
|
|
ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &merged,
|
|
const UnitCell &cell,
|
|
const std::string &model_path,
|
|
const std::string &output_prefix,
|
|
Logger &logger,
|
|
const gemmi::SpaceGroup *data_space_group,
|
|
bool probe_indexing_ambiguity,
|
|
size_t nthreads,
|
|
double wavelength_A) {
|
|
ModelValidationResult result;
|
|
result.model_path = model_path;
|
|
|
|
// --- read the atomic model ---
|
|
ModelState mdl;
|
|
gemmi::Structure &st = mdl.st; // the real model, the one the maps and the report describe
|
|
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), the model's group is a
|
|
// CANDIDATE for the label the reflections are written under. Only a candidate: this is arithmetic
|
|
// on two group numbers and says nothing about whether the model belongs to this crystal, and a
|
|
// model that does not is exactly as capable of rewriting the label as one that does. Whether it is
|
|
// taken up is settled at the end, on the fit against its own null and on the anomalous map.
|
|
//
|
|
// 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 && data_space_group->number != sg->number) {
|
|
if (data_space_group->is_enantiomorphic() && sg->is_enantiomorphic()) {
|
|
gemmi::GroupOps eops = data_space_group->operations();
|
|
eops.change_basis_forward(data_space_group->change_of_hand_op());
|
|
const gemmi::SpaceGroup *enant = gemmi::find_spacegroup_by_ops(eops);
|
|
if (enant && enant->number == sg->number) {
|
|
result.model_enantiomorph_candidate = true;
|
|
logger.Info("Model validation: data space group {} is the enantiomorph of the model {}; "
|
|
"the model's group is a candidate for the label, pending the fit - a change "
|
|
"of label only, since the two groups index identically and reindexing would "
|
|
"flip the anomalous differences",
|
|
data_space_group->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, plus a flat bulk-solvent mask -> Fmask.
|
|
// A lambda because the rigid-body step below moves the model and then needs both again, and it
|
|
// takes the state to work on so that a null replicate can run it on its own copy. ---
|
|
auto compute_model_factors = [&](ModelState &ms) {
|
|
gemmi::DensityCalculator<Table, float> dc;
|
|
dc.d_min = d_min;
|
|
dc.rate = 1.5;
|
|
dc.set_grid_cell_and_spacegroup(ms.st);
|
|
dc.set_refmac_compatible_blur(ms.st.models[0]);
|
|
dc.put_model_density_on_grid(ms.st.models[0]);
|
|
ms.fcalc = gemmi::transform_map_to_f_phi(dc.grid, true)
|
|
.prepare_asu_data(dc.d_min, dc.blur, false, false, false);
|
|
|
|
// 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, ms.st.models[0]);
|
|
ms.fmask = gemmi::transform_map_to_f_phi(mask_grid, true).prepare_asu_data(dc.d_min, 0);
|
|
};
|
|
compute_model_factors(mdl);
|
|
|
|
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;
|
|
gemmi::AsuData<gemmi::ValueSigma<float>> fobs; // what it was fitted to
|
|
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 = [&](ModelState &ms, const std::vector<MergedReflection> &obs_in) -> Fit {
|
|
Fit out;
|
|
out.fmodel = ms.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 = out.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, &ms.fmask);
|
|
scaling.fit_isotropic_b_approximately();
|
|
scaling.fit_parameters();
|
|
scaling.scale_data(out.fmodel, &ms.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.
|
|
// (The sigma_A weighting further down is a different thing and does not reopen this: it never
|
|
// rescales Fobs, and it leaves the R-factors below untouched. It does weight the map
|
|
// coefficients per shell and per dataset - see the note where it is computed.)
|
|
|
|
// --- R-work / R-free ---
|
|
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);
|
|
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; }
|
|
}
|
|
out.r_work = den_w > 0 ? num_w / den_w : 1;
|
|
out.r_free = den_f > 0 ? num_f / den_f : 1;
|
|
return out;
|
|
};
|
|
|
|
// --- the model placed against one observed set: fit, then a rigid-body placement ---
|
|
// The real model and every null replicate below go through this same lambda, so the two are
|
|
// comparable - a null that was not placed would be the null of a weaker procedure than the one it
|
|
// is there to judge. It leaves the model where it ends up, and the caller puts it back.
|
|
struct Placement {
|
|
Fit fit;
|
|
bool rb_applied = false;
|
|
double rb_angle_deg = 0, rb_shift_A = 0, r_free_before_rb = 0;
|
|
};
|
|
auto place_and_fit = [&](ModelState &ms, const std::vector<MergedReflection> &obs_in) -> Placement {
|
|
Placement out;
|
|
out.fit = fit_model(ms, obs_in);
|
|
|
|
// --- rigid-body placement of the model in the data cell ---
|
|
// Re-fractionalizing the model above puts it in the right box but not in the right place: a
|
|
// non-isomorphous cell squeezes the box without moving the body inside it, and the body's own
|
|
// position differs from crystal to crystal anyway. Six parameters recover that - a rotation about
|
|
// the model's centroid and a translation - which is all a fragment-screening model needs, since it
|
|
// arrives already solved. The step is committed only if the FREE reflections, which the refinement
|
|
// never saw, agree that it helped; on data the model cannot be placed against, the model stays
|
|
// exactly where it was read.
|
|
gemmi::AsuData<gemmi::ValueSigma<float>> fobs_work;
|
|
fobs_work.unit_cell_ = ucell;
|
|
fobs_work.spacegroup_ = sg;
|
|
for (const auto &hv : out.fit.fobs.v) {
|
|
const auto it = out.fit.obs_by_hkl.find(hkl_key(hv.hkl));
|
|
if (it != out.fit.obs_by_hkl.end() && !it->second.second)
|
|
fobs_work.v.push_back(hv);
|
|
}
|
|
const std::vector<gemmi::Position> before = ModelPositions(ms.st.models[0]);
|
|
const RigidBodyRefineResult rb =
|
|
RefineRigidBody(ms.st.models[0], ucell, *sg, fobs_work, d_min, logger);
|
|
if (rb.converged) {
|
|
compute_model_factors(ms);
|
|
Fit moved = fit_model(ms, obs_in);
|
|
const bool commit = moved.r_free < out.fit.r_free;
|
|
logger.Info("Model validation: rigid body held-out R-free {:.4f} -> {:.4f} => {}",
|
|
out.fit.r_free, moved.r_free, commit ? "committed" : "rejected, model put back");
|
|
if (commit) {
|
|
out.rb_applied = true;
|
|
out.rb_angle_deg = rb.angle_deg;
|
|
out.rb_shift_A = rb.shift_A;
|
|
out.r_free_before_rb = out.fit.r_free;
|
|
out.fit = std::move(moved);
|
|
} else {
|
|
SetModelPositions(ms.st.models[0], before);
|
|
compute_model_factors(ms);
|
|
}
|
|
}
|
|
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.
|
|
//
|
|
// The candidates are a property of the DATA's lattice and point group: it is the observed
|
|
// intensities that are relabelled. Asking the model's group instead enumerates nothing at all
|
|
// wherever the two groups differ - measured, a model in P4(3)2(1)2 against data merged in P4(3)
|
|
// probes zero operators - which is exactly the case the probe exists for.
|
|
const gemmi::SpaceGroup &ambiguity_sg = data_space_group != nullptr ? *data_space_group : *sg;
|
|
const std::vector<gemmi::Op> reindex_ops =
|
|
probe_indexing_ambiguity ? ReindexAmbiguityOperators(cell, ambiguity_sg)
|
|
: std::vector<gemmi::Op>{};
|
|
|
|
struct IndexingProbe {
|
|
gemmi::Op op = gemmi::Op::identity(); // the lowest-R-free relabelling
|
|
double margin = 0; // by how much in R-free it leads the runner-up
|
|
};
|
|
auto probe_indexing = [&](ModelState &ms) {
|
|
IndexingProbe out;
|
|
std::vector<double> r_free{fit_model(ms, obs).r_free}; // identity first, then the twin laws
|
|
double best_r_free = r_free.front();
|
|
for (const auto &op : reindex_ops) {
|
|
const double cand = fit_model(ms, ReindexReflections(obs, op)).r_free;
|
|
r_free.push_back(cand);
|
|
if (cand < best_r_free) {
|
|
best_r_free = cand;
|
|
out.op = op;
|
|
}
|
|
}
|
|
// 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(r_free.begin(), r_free.end());
|
|
out.margin = r_free[1] - r_free[0];
|
|
return out;
|
|
};
|
|
|
|
const std::vector<gemmi::Position> as_read = ModelPositions(st.models[0]);
|
|
|
|
IndexingProbe indexing;
|
|
std::vector<MergedReflection> reindexed;
|
|
const std::vector<MergedReflection> *obs_model = &obs;
|
|
if (!reindex_ops.empty()) {
|
|
indexing = probe_indexing(mdl);
|
|
if (!(indexing.op == gemmi::Op::identity())) {
|
|
reindexed = ReindexReflections(obs, indexing.op);
|
|
obs_model = &reindexed;
|
|
}
|
|
}
|
|
Placement real = place_and_fit(mdl, *obs_model);
|
|
|
|
// --- is there anything for the model to decide? ---
|
|
// There are only two: the space-group label, where the model asserts the other enantiomorph, and
|
|
// the indexing, where a relabelling of the data fits it better. A model that asserts neither -
|
|
// already in the group the data were merged in, and preferring their indexing - has made no claim
|
|
// that needs arbitrating, and the null below would be several seconds spent gating a decision
|
|
// nobody is making. That is the isomorphous case a screening campaign is made of, and it is the
|
|
// one that has to be fast. Nothing else changes: the R-factors, the maps and the placement are
|
|
// computed and reported exactly as in every other case.
|
|
const bool decision_pending = result.model_enantiomorph_candidate
|
|
|| !(indexing.op == gemmi::Op::identity());
|
|
if (decision_pending) {
|
|
// --- the null: this model, this data, in random orientations ---
|
|
// An R-factor on its own says nothing about whether a model belongs to a crystal. What a model
|
|
// that explains nothing reaches against these data depends on its atom count and its B-factors as
|
|
// much as on the data, so a threshold on R calibrated on one dataset misjudges the next; and the
|
|
// classical random-structure value does not apply either, because the scale here is fitted to
|
|
// minimise the very sum the R is made of. The only null that fits is made of the same model: it is
|
|
// reoriented about its own centroid and run through the identical fit and rigid-body placement,
|
|
// and the real fit is asked how far above the resulting distribution it sits.
|
|
//
|
|
// R-work, not R-free: nothing is refined against the working set here - the scale has four
|
|
// parameters and the placement six - so R-work carries no optimism, and it is decided on an order
|
|
// of magnitude more reflections than R-free.
|
|
//
|
|
// The replicates are fitted to the data as merged even where a reindexing won above: a reindexing
|
|
// is a relabelling of the same intensities, and a model in a random orientation has no more to do
|
|
// with one labelling than with the other.
|
|
logger.Info("Model validation: scoring the fit against a null of {} random placements of the same "
|
|
"model", NULL_REPLICATES);
|
|
// The orientations are drawn here, up front and in order, and not inside the loop: the
|
|
// replicates run at the same time below, and a draw taken on whichever thread reached the
|
|
// generator first would make the answer depend on the scheduling. Replicate i gets rotation i
|
|
// on every run, so the sigma reported is a property of the data and not of the machine.
|
|
std::vector<gemmi::Mat33> null_rotation;
|
|
{
|
|
std::mt19937 rng(NULL_SEED);
|
|
for (int i = 0; i < NULL_REPLICATES; i++)
|
|
null_rotation.push_back(random_rotation(rng));
|
|
}
|
|
// Each replicate works on its own copy of the model and of its structure factors, and the
|
|
// replicates share nothing else, so they run concurrently and the real model - which the
|
|
// maps and the atom-density readout below are still taken from - is never moved.
|
|
std::vector<double> null_r_work(NULL_REPLICATES), null_margin(NULL_REPLICATES);
|
|
ParallelFor(NULL_REPLICATES, nthreads, [&](int i) {
|
|
ModelState rep;
|
|
rep.st = mdl.st;
|
|
SetModelPositions(rep.st.models[0], as_read);
|
|
reorient_about_centroid(rep.st.models[0], null_rotation[i]);
|
|
compute_model_factors(rep);
|
|
if (!reindex_ops.empty())
|
|
null_margin[i] = probe_indexing(rep).margin;
|
|
null_r_work[i] = place_and_fit(rep, obs).fit.r_work;
|
|
});
|
|
const auto [null_mean, null_sd] = mean_sd(null_r_work);
|
|
result.fit_tested = true;
|
|
result.null_replicates = NULL_REPLICATES;
|
|
result.null_r_work_mean = null_mean;
|
|
result.null_r_work_sd = null_sd;
|
|
|
|
auto record_fit = [&](const Placement &p) {
|
|
result.r_work_sigma = null_sd > 0 ? (null_mean - p.fit.r_work) / null_sd : 0.0;
|
|
result.model_fits = result.r_work_sigma >= MODEL_FIT_SIGMA;
|
|
};
|
|
record_fit(real);
|
|
|
|
// --- does the model get to decide anything? ---
|
|
// The R-factors, the maps and the rigid-body placement are statements about the MODEL and are
|
|
// always computed and always reported, fit or no fit: a model that does not describe these data
|
|
// still has an R against them, and that is the negative result. What the gate below controls is
|
|
// the two decisions that rewrite the DATA - the space-group label and the indexing - which a
|
|
// wrong model must not be able to make.
|
|
if (!reindex_ops.empty()) {
|
|
result.indexing_probed = true;
|
|
result.indexing_margin = indexing.margin;
|
|
const auto [margin_mean, margin_sd] = mean_sd(null_margin);
|
|
result.indexing_margin_null_mean = margin_mean;
|
|
result.indexing_margin_null_sd = margin_sd;
|
|
result.indexing_margin_sigma =
|
|
margin_sd > 0 ? (indexing.margin - margin_mean) / margin_sd : 0.0;
|
|
// The margin is judged against the margin a random placement produces, not against a value.
|
|
// A random model also picks a winner, and on these data it picks one by a comparable lead
|
|
// (measured), so the raw margin says nothing on its own.
|
|
const bool decided = result.model_fits && result.indexing_margin_sigma >= MODEL_FIT_SIGMA;
|
|
if (decided)
|
|
result.indexing_op = indexing.op;
|
|
logger.Info("Model validation: probed {} indexing solution(s) against the model; the winner "
|
|
"leads the runner-up by {:.4f} in R-free, against {:.4f} +- {:.4f} for a random "
|
|
"placement of the same model ({:+.2f} sigma) => {}",
|
|
reindex_ops.size() + 1, indexing.margin, margin_mean, margin_sd,
|
|
result.indexing_margin_sigma,
|
|
!decided ? "not decided; the data keep the indexing they were merged in"
|
|
: (result.indexing_op == gemmi::Op::identity()
|
|
? "kept the current indexing" : "reindexed to the model's"));
|
|
if (!decided && !(indexing.op == gemmi::Op::identity())) {
|
|
// The R-factors and the maps have to describe the reflections the files carry, and those
|
|
// are now the ones the merge produced, so the fit is remade on them.
|
|
SetModelPositions(st.models[0], as_read);
|
|
compute_model_factors(mdl);
|
|
real = place_and_fit(mdl, obs);
|
|
record_fit(real);
|
|
}
|
|
}
|
|
logger.Info("Model validation: R-work {:.4f} against a null of {:.4f} +- {:.4f} = {:+.2f} sigma "
|
|
"=> the model {} these data",
|
|
real.fit.r_work, null_mean, null_sd, result.r_work_sigma,
|
|
result.model_fits ? "FITS" : "DOES NOT FIT");
|
|
} else {
|
|
logger.Info("Model validation: the model is already in the space group the data were merged in "
|
|
"and prefers their indexing, so it claims nothing about the written reflections and "
|
|
"no null was run; R-work {:.4f}, R-free {:.4f}", real.fit.r_work, real.fit.r_free);
|
|
}
|
|
|
|
Fit &best = real.fit;
|
|
result.rigid_body_applied = real.rb_applied;
|
|
result.rigid_body_angle_deg = real.rb_angle_deg;
|
|
result.rigid_body_shift_A = real.rb_shift_A;
|
|
result.r_free_before_rigid_body = real.r_free_before_rb;
|
|
|
|
gemmi::AsuData<std::complex<float>> &fmodel = best.fmodel;
|
|
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;
|
|
|
|
// --- sigma_A weighting: the maps are 2mFo-DFc and mFo-DFc, not 2Fo-Fc and Fo-Fc ---
|
|
// m and D come from a maximum-likelihood sigma_A per resolution shell, so a shell the model
|
|
// describes badly is damped rather than carried into the map at full weight, and the difference
|
|
// map is correspondingly less biased towards the model that made its phases.
|
|
//
|
|
// m and D are estimated on THIS dataset, so two datasets of one crystal form get slightly
|
|
// different weights, and to that extent their maps are no longer scaled identically - the same
|
|
// property the scaling above deliberately protects. It is kept anyway: the difference between two
|
|
// datasets' sigma_A curves is the difference in how well the model explains each of them, which is
|
|
// real and is what a screening campaign is looking for, and a PanDDA-style analysis consumes
|
|
// 2mFo-DFc maps and normalizes each dataset's map against the ensemble before comparing them. The
|
|
// per-reflection FOM is written to the MTZ so the weighting can be undone.
|
|
// Following Read (1986) Acta Cryst. A42, 140-149
|
|
struct MapTerm { gemmi::Miller hkl; double fo, fc, phi; bool free, centric; };
|
|
std::vector<MapTerm> terms;
|
|
std::vector<SigmaAReflection> sa_input;
|
|
for (const auto &hv : fmodel.v) {
|
|
const auto it = obs_by_hkl.find(hkl_key(hv.hkl));
|
|
if (it == obs_by_hkl.end()) continue;
|
|
const double Fo = it->second.first;
|
|
const double Fc = std::abs(hv.value);
|
|
const bool centric = gops.is_reflection_centric(hv.hkl);
|
|
terms.push_back({hv.hkl, Fo, Fc, std::arg(hv.value), it->second.second, centric});
|
|
sa_input.push_back({Fo, Fc, ucell.calculate_1_d2(hv.hkl), gops.epsilon_factor(hv.hkl),
|
|
centric, it->second.second});
|
|
}
|
|
const SigmaAResult sigma_a = EstimateSigmaA(sa_input, ucell);
|
|
result.mean_fom = sigma_a.mean_fom;
|
|
result.sigma_a_shells = sigma_a.shells;
|
|
|
|
gemmi::AsuData<std::complex<float>> map2fofc, mapfofc;
|
|
map2fofc.unit_cell_ = ucell; map2fofc.spacegroup_ = sg;
|
|
mapfofc.unit_cell_ = ucell; mapfofc.spacegroup_ = sg;
|
|
std::vector<float> fwt(terms.size()), delfwt(terms.size());
|
|
for (size_t i = 0; i < terms.size(); i++) {
|
|
const double m = sigma_a.weight[i].m, D = sigma_a.weight[i].d;
|
|
// A centric reflection's phase is either exactly right or 180 degrees wrong, never in
|
|
// between, so its bias-free coefficient is mFo and not 2mFo - DFc.
|
|
fwt[i] = static_cast<float>(terms[i].centric ? m * terms[i].fo
|
|
: 2 * m * terms[i].fo - D * terms[i].fc);
|
|
delfwt[i] = static_cast<float>(m * terms[i].fo - D * terms[i].fc);
|
|
const std::complex<float> ph = std::polar(1.0f, static_cast<float>(terms[i].phi));
|
|
map2fofc.v.push_back({terms[i].hkl, fwt[i] * ph});
|
|
mapfofc.v.push_back({terms[i].hkl, delfwt[i] * ph});
|
|
}
|
|
|
|
// --- write the maps and score the 2mFo-DFc 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 != nullptr ? data_space_group : 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) - 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);
|
|
}
|
|
}
|
|
|
|
// --- the enantiomorph, now that there is something to decide it on ---
|
|
// Two things have to hold before the model's hand is written on these data. The model has to
|
|
// describe them at all - the label is an assertion about the crystal, and a model that fits no
|
|
// better than its own random placements is in no position to make one. And the anomalous map, the
|
|
// only measurement here that is sensitive to the hand at all, must not contradict it: R-free
|
|
// cannot (measured, inverting the model through the origin moves R-work by less than 1e-4, because
|
|
// |F(h)| of the inverted structure is |F(-h)| = |F(h)| on Friedel-averaged data), so where the
|
|
// anomalous differences do say something they say it alone, and they get a veto.
|
|
result.adopted_model_enantiomorph = result.model_enantiomorph_candidate && result.model_fits
|
|
&& !result.anomalous_hands_disagree;
|
|
if (result.model_enantiomorph_candidate && !result.adopted_model_enantiomorph)
|
|
logger.Warning("Model validation: the model asserts the enantiomorph {} against the data's {}, "
|
|
"and that assertion is NOT taken up: {}. The reflections are written in the "
|
|
"group they were merged in",
|
|
sg->hm, data_space_group != nullptr ? data_space_group->hm : "?",
|
|
result.anomalous_hands_disagree
|
|
? "the anomalous density at the model's atoms is inverted"
|
|
: "the model does not fit these data");
|
|
|
|
// --- MTZ of map coefficients so the maps can be re-opened / rebuilt in Coot etc. ---
|
|
try {
|
|
gemmi::Mtz mtz(true);
|
|
// The group the REFLECTIONS end up in, which is the model's only where its hand was adopted -
|
|
// AdoptModelFrame decides this the same way a few lines below. An enantiomorphic pair indexes
|
|
// identically, so the coefficients are the same numbers either way and only the label moves;
|
|
// but the two groups have different screw translations, so a reader that expands symmetry out
|
|
// of this file works in the wrong one if the label disagrees with the .mtz beside it.
|
|
mtz.spacegroup = result.adopted_model_enantiomorph
|
|
? sg
|
|
: (data_space_group != nullptr ? data_space_group : sg);
|
|
mtz.set_cell_for_all(ucell);
|
|
mtz.add_dataset("model_validation");
|
|
mtz.datasets.back().wavelength = wavelength_A;
|
|
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);
|
|
// The figure of merit the coefficients carry, so the weighting can be read off - and undone -
|
|
// from the file rather than having to be taken on trust.
|
|
mtz.add_column("FOM", 'W', -1, -1, false);
|
|
mtz.add_column("FREE", 'I', -1, -1, false);
|
|
std::vector<float> data;
|
|
for (size_t i = 0; i < terms.size(); i++) {
|
|
const auto phi_deg = static_cast<float>(terms[i].phi * 180.0 / PI);
|
|
data.insert(data.end(), {static_cast<float>(terms[i].hkl[0]), static_cast<float>(terms[i].hkl[1]),
|
|
static_cast<float>(terms[i].hkl[2]),
|
|
static_cast<float>(terms[i].fo), static_cast<float>(terms[i].fc),
|
|
phi_deg,
|
|
fwt[i], phi_deg,
|
|
delfwt[i], phi_deg,
|
|
static_cast<float>(sigma_a.weight[i].m),
|
|
terms[i].free ? 0.0f : 1.0f});
|
|
}
|
|
mtz.nreflections = static_cast<int>(terms.size());
|
|
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;
|
|
// The placed coordinates, for the caller to write out beside the maps once the frame is settled.
|
|
result.placed_model = std::make_shared<gemmi::Structure>(st);
|
|
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 2mFo-DFc density at atom centres = {:.2f} sigma", result.mean_atom_density_sigma);
|
|
logger.Info("Model validation: sigma_A weighting over {} resolution shell(s) (estimated on the {} "
|
|
"free reflections): sigma_A {:.2f} at low resolution, {:.2f} at high, mean FOM {:.3f}",
|
|
sigma_a.shells, sigma_a.free_reflections, sigma_a.sigma_a_lowest_shell,
|
|
sigma_a.sigma_a_highest_shell, sigma_a.mean_fom);
|
|
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);
|
|
if (result.fit_tested && !result.model_fits)
|
|
logger.Info("Model validation: the model was rejected, so nothing downstream moved - the "
|
|
"reflection files are byte for byte the ones a run with no model would have "
|
|
"written. The R-factors and the maps above still describe this model against "
|
|
"these data, and are the negative result rather than a failure");
|
|
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
|
|
|
|
const gemmi::SpaceGroup *AdoptModelFrame(const ModelValidationResult &validation,
|
|
std::vector<MergedReflection> &merged,
|
|
const gemmi::SpaceGroup &data_space_group,
|
|
bool merge_friedel,
|
|
Logger &logger) {
|
|
const gemmi::SpaceGroup *space_group = &data_space_group;
|
|
// A model that could not be read, or that did not fit, decides nothing: the reflections are
|
|
// written exactly as a run with no model at all would have written them. Both flags are already
|
|
// folded into the fields below, and the guard states the rule where the frame is actually adopted.
|
|
if (!validation.ok || !validation.model_fits)
|
|
return space_group;
|
|
|
|
// 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 = gemmi::find_spacegroup_by_number(validation.model_space_group_number);
|
|
logger.Info("Model validation: the written reflections take the model's enantiomorph, {} ({}), "
|
|
"as a label - no reflection moved",
|
|
space_group ? space_group->short_name() : "?", validation.model_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, merge_friedel);
|
|
logger.Info("Model validation: the written reflections take the model's indexing, reindexed by {}",
|
|
hkl_triplet(validation.indexing_op));
|
|
}
|
|
return space_group;
|
|
}
|
|
|
|
std::vector<MergedReflection> ModelReferenceIntensities(const std::string &model_path,
|
|
const std::optional<UnitCell> &cell,
|
|
const gemmi::SpaceGroup *space_group,
|
|
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 != nullptr)
|
|
st.spacegroup_hm = space_group->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;
|
|
}
|