leonarski_fandClaude Opus 5.5
1fd28672bb
The null is nine fixed-seed random orientations of the model, each placed by the same rigid body as
the real fit and fitted to the data AS MERGED. Where the indexing probe reindexed the data to the
model's, the model's solution against the merged data is its twin-related orientation - and more
generally the crystal looks the same from every orientation the lattice's rotations (the space
group's, composed with the twin laws) take the model to. A draw within the rigid body's reach of one
of those is walked onto it and is not a sample of the null.
On 5epe (F 2 3, twin law y,x,-z, 24 equivalent orientations) replicate 8 of 9 was drawn 7.9 deg from
one of them; the rigid body moved it 7.9 deg, held-out R-free 0.53 -> 0.20, R-work 0.197 against
0.56-0.58 for the other eight. The null's sd became 0.125 and the crystal's own model read +0.01
sigma, DOES NOT FIT, so the indexing it had decided at +25.8 sigma was refused and the files kept the
merged indexing (R-free 0.52). This happens on every binary whose indexer picks the other hand.
The fix is at the draw, not the score: a draw whose angle to the nearest equivalent orientation is
under the rigid body's reach is drawn again from the same generator, so everywhere no draw falls
there the nine rotations - and every output - are exactly the ones before. The reach is the rotation
about the centroid that moves the model's rms-radius atom by the ladder's first zone (6 A, or d_min
where coarser): 18.7 deg on 5epe, whose other draws sat at 21.7 deg and beyond and did not converge.
Twin laws come from the lattice (ReindexAmbiguityOperators on the data cell and group) regardless of
whether the probe ran. Redraws are capped at 1000 so that a model small enough for the reach to cover
the whole of SO(3) takes its draws as they come instead of looping.
Scoring the real model against the best of both hands instead was not taken: the real fit is
already scored in the hand that wins, and the contamination is in the null, not in the real fit.
Measured, battery commands with --model, rc173 27e154e5d vs this:
5epe 1 redraw; DOES NOT FIT +0.01 -> FITS +55.86 sigma (null 0.5730 +- 0.0072); indexing
decided at +40.4 sigma and reindexed y,x,-z; R-free 0.5246 -> 0.1749 (deposited 0.1671)
9fhc 1 redraw (reach 8.8 deg; the replaced draw had not converged); files byte-identical,
MODEL_FIT_SIGMA +201.27 -> +199.39, verdict unchanged
9hnc, 8sa8, 8xtg, 9bn8: no redraw; p.hkl, p.mtz, p.cif (bar the audit date), p_maps.mtz,
p_model.cif and the model-validation report byte-identical
[ModelValidation] tests pass.
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
1485 lines
84 KiB
C++
1485 lines
84 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 "ModelScaling.h"
|
|
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
#include <complex>
|
|
#include <array>
|
|
#include <atomic>
|
|
#include <future>
|
|
#include <numeric>
|
|
#include <optional>
|
|
#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> // get_size_for_hkl, get_f_phi_on_grid (the maps)
|
|
#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/CorrelationCoefficient.h"
|
|
#include "../common/JFJochMath.h" // PI (M_PI is not standard, and MSVC does not define it)
|
|
#include "../common/Logger.h"
|
|
#include "../image_analysis/scale_merge/ReindexAmbiguity.h" // ReindexReflections
|
|
#include "../image_analysis/scale_merge/CrystalSetting.h" // CellMappingOperators
|
|
#include "ModelFFT.h"
|
|
#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 MapFromFPhi(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 many draws the null may throw away for lying within the rigid body's reach of the model's own
|
|
// orientation. Only a model small enough for that reach to cover most orientations gets near it, and
|
|
// there the remaining draws are taken as they come rather than the run searching forever.
|
|
constexpr int NULL_MAX_REDRAWS = 1000;
|
|
|
|
// 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. Measured at nine replicates, the crystal's own model reads +17.7 sigma and an
|
|
// unrelated protein +1.7, so the cut sits in a gap an order of magnitude wider than the sd it is
|
|
// measured in.
|
|
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)};
|
|
}
|
|
|
|
// A rotation of the lattice, given on fractional coordinates as a symmetry operator is, in Cartesian
|
|
// coordinates - the frame the null's orientations are drawn in.
|
|
gemmi::Mat33 cartesian_rotation(const gemmi::Op &op, const gemmi::UnitCell &cell) {
|
|
gemmi::Mat33 f;
|
|
for (int i = 0; i < 3; i++)
|
|
for (int j = 0; j < 3; j++)
|
|
f.a[i][j] = static_cast<double>(op.rot[i][j]) / gemmi::Op::DEN;
|
|
return cell.orth.mat.multiply(f).multiply(cell.frac.mat);
|
|
}
|
|
|
|
// The angle, in degrees, of the rotation from `rot` to the nearest of `orientations`.
|
|
double angle_to_nearest_deg(const gemmi::Mat33 &rot, const std::vector<gemmi::Mat33> &orientations) {
|
|
double best = 180.0;
|
|
for (const gemmi::Mat33 &o : orientations) {
|
|
const gemmi::Mat33 d = rot.multiply(o.transpose());
|
|
const double c = std::clamp((d.a[0][0] + d.a[1][1] + d.a[2][2] - 1.0) / 2.0, -1.0, 1.0);
|
|
best = std::min(best, std::acos(c) * 180.0 / PI);
|
|
}
|
|
return best;
|
|
}
|
|
|
|
// 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);
|
|
}
|
|
|
|
// Move a model into another cell by keeping its fractional coordinates. Correct only when the two
|
|
// cells describe the same axes in the same order, which is what the probe below is there to arrange.
|
|
void refractionalize_into(gemmi::Structure &st, const gemmi::UnitCell &target) {
|
|
const gemmi::UnitCell from = 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(from.fractionalize(a.pos));
|
|
st.cell = target;
|
|
}
|
|
|
|
// --- putting the model in the data's description of the lattice ---------------------------------
|
|
//
|
|
// A model arrives in the cell its depositor chose and rugnux indexes in the cell its own reduction
|
|
// chose, and the two are often different descriptions of the SAME lattice: I-centred where the other
|
|
// is C-centred, unique axis c where the other took b, a cyclic permutation of an orthorhombic cell.
|
|
// The space-group NUMBER is identical in every one of those, so no comparison of numbers can see it,
|
|
// and re-fractionalizing straight across such a pair scrambles the model. The rigid body further
|
|
// down cannot undo it either - six parameters about a centroid are not a change of basis - so the
|
|
// run would otherwise report a placement R-free near 0.6 for data that are perfectly good.
|
|
|
|
// Two candidates that differ by a rotation the model's own group already has describe the same
|
|
// structure, so only one of each class is worth scoring. On a holohedral cell that collapses two
|
|
// dozen candidates to one, which is what keeps the probe free on the ordinary isomorphous run.
|
|
bool same_frame(const gemmi::Op &a, const gemmi::Op &b, const gemmi::GroupOps &gops) {
|
|
const gemmi::Op::Rot d = a.inverse().combine(b).rot;
|
|
for (const gemmi::Op &s : gops.sym_ops)
|
|
if (s.rot == d)
|
|
return true;
|
|
return false;
|
|
}
|
|
|
|
// Do two groups have the same rotations - the same point group in the same setting? I 2 3 and
|
|
// I 21 3 do, and so does an enantiomorphic pair; P 4 and P 4 2 2 do not, nor do P 1 2 1 and P 1 1 2.
|
|
bool same_rotations(const gemmi::GroupOps &a, const gemmi::GroupOps &b) {
|
|
if (a.sym_ops.size() != b.sym_ops.size())
|
|
return false;
|
|
for (const gemmi::Op &op : a.sym_ops)
|
|
if (b.find_by_rotation(op.rot) == nullptr)
|
|
return false;
|
|
return true;
|
|
}
|
|
|
|
// A change of basis is a matrix AND an origin shift, and the shift is not optional: an odd
|
|
// permutation of a screw-axis group lands on a group that is the same group on a moved origin, whose
|
|
// operator list GEMMI cannot name because it compares those lists exactly. Swapping b and c in
|
|
// P 21 21 21 is the everyday example - it needs (1/4, 1/4, 1/4) before it reads as P 21 21 21 again.
|
|
// So where the bare matrix names nothing, the shift that makes it name something is searched for, on
|
|
// the twelfths every crystallographic origin shift lies on. The order tries the common shifts first,
|
|
// so the search almost always ends on one of its first few candidates.
|
|
gemmi::Op with_origin_shift(const gemmi::SpaceGroup *sg, const gemmi::Op &op,
|
|
const gemmi::SpaceGroup **named) {
|
|
static const int TWELFTHS[] = {0, 6, 3, 9, 4, 8, 2, 10, 1, 5, 7, 11};
|
|
for (int i : TWELFTHS)
|
|
for (int j : TWELFTHS)
|
|
for (int k : TWELFTHS) {
|
|
gemmi::Op shifted = op;
|
|
shifted.tran = {i * gemmi::Op::DEN / 12, j * gemmi::Op::DEN / 12,
|
|
k * gemmi::Op::DEN / 12};
|
|
gemmi::GroupOps gops = sg->operations();
|
|
gops.change_basis_forward(shifted);
|
|
if (const gemmi::SpaceGroup *found = gemmi::find_spacegroup_by_ops(gops)) {
|
|
*named = found;
|
|
return shifted;
|
|
}
|
|
}
|
|
*named = nullptr;
|
|
return op;
|
|
}
|
|
|
|
// Put a model through a change of basis: coordinates, cell and space group together. Returns false -
|
|
// leaving the model untouched - when no origin shift makes the transformed group one GEMMI can name,
|
|
// which is how a basis that would leave a standard setting is refused rather than adopted.
|
|
bool change_model_basis(gemmi::Structure &st, const gemmi::SpaceGroup *&sg, gemmi::Op &op) {
|
|
const gemmi::SpaceGroup *moved = nullptr;
|
|
op = with_origin_shift(sg, op, &moved);
|
|
if (moved == nullptr)
|
|
return false;
|
|
gemmi::UnitCell old_cell = st.cell;
|
|
gemmi::Op rot_only = op;
|
|
rot_only.tran = {0, 0, 0}; // the cell follows the axes; only the atoms feel the origin shift
|
|
const gemmi::UnitCell new_cell = old_cell.changed_basis_forward(rot_only, false);
|
|
for (gemmi::Model &m : st.models)
|
|
for (gemmi::Chain &ch : m.chains)
|
|
for (gemmi::Residue &r : ch.residues)
|
|
for (gemmi::Atom &a : r.atoms) {
|
|
const gemmi::Fractional f = old_cell.fractionalize(a.pos);
|
|
const std::array<double, 3> t = op.apply_to_xyz({{f.x, f.y, f.z}});
|
|
a.pos = new_cell.orthogonalize(gemmi::Fractional(t[0], t[1], t[2]));
|
|
}
|
|
st.cell = new_cell;
|
|
st.spacegroup_hm = moved->xhm();
|
|
sg = moved;
|
|
return true;
|
|
}
|
|
|
|
// The coarse shell the frame is decided on. A frame that is wrong is wrong at low resolution, so the
|
|
// probe never goes near the resolution the real fit uses - that is what makes trying every candidate
|
|
// affordable. The very lowest resolution is left out with it: there a bulk solvent this scorer does
|
|
// not model would dominate, equally for every candidate, and only add noise to the comparison.
|
|
constexpr double FRAME_PROBE_D_MIN = 3.5;
|
|
constexpr double FRAME_PROBE_D_MAX = 8.0;
|
|
constexpr size_t FRAME_PROBE_MIN_REFLECTIONS = 200;
|
|
|
|
// R of the model against the observed amplitudes after an overall scale and an isotropic B, over
|
|
// that coarse shell. Only the ranking is ever used, never the value.
|
|
double frame_probe_r(const gemmi::Structure &st, const gemmi::SpaceGroup *sg,
|
|
const std::vector<MergedReflection> &obs, double d_min) {
|
|
const double probe_d_min = std::max(d_min, FRAME_PROBE_D_MIN);
|
|
gemmi::DensityCalculator<Table, float> dc;
|
|
dc.d_min = probe_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 =
|
|
MapToFPhi(dc.grid)
|
|
.prepare_asu_data(probe_d_min, dc.blur, false, false, false);
|
|
|
|
gemmi::GroupOps gops = sg->operations();
|
|
gemmi::ReciprocalAsu asu(sg);
|
|
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
|
|
fobs.unit_cell_ = st.cell;
|
|
fobs.spacegroup_ = sg;
|
|
// The low-resolution cap is dropped, rather than the probe abandoned, when a small cell does not
|
|
// put enough reflections in the shell.
|
|
for (double d_max : {FRAME_PROBE_D_MAX, 1e9}) {
|
|
fobs.v.clear();
|
|
for (const MergedReflection &r : obs) {
|
|
if (std::isnan(r.F) || r.d < probe_d_min || r.d > d_max)
|
|
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}});
|
|
}
|
|
if (fobs.v.size() >= FRAME_PROBE_MIN_REFLECTIONS)
|
|
break;
|
|
}
|
|
if (fobs.v.empty())
|
|
return 1.0;
|
|
fobs.ensure_asu();
|
|
fobs.ensure_sorted();
|
|
gemmi::Scaling<float> scaling(st.cell, sg);
|
|
scaling.use_solvent = false;
|
|
scaling.prepare_points(fcalc, fobs, nullptr);
|
|
scaling.fit_isotropic_b_approximately();
|
|
return scaling.calculate_r_factor();
|
|
}
|
|
|
|
|
|
} // 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,
|
|
const std::vector<float> &report_shell_d_min) {
|
|
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;
|
|
}
|
|
|
|
const gemmi::UnitCell data_cell = cell; // UnitCell -> gemmi::UnitCell
|
|
|
|
// --- indexing (merohedral) ambiguity: the candidates ---
|
|
// 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. They are proper
|
|
// rotations of the lattice, never the inversion, so relabelling by one keeps the hand.
|
|
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>{};
|
|
|
|
// --- same lattice, whose description of it? ---
|
|
// Scored, never asserted. The identity is always among the candidates, so a probe that finds
|
|
// nothing better than leaving the model where it is leaves it there; and the ordinary isomorphous
|
|
// run - where the only operators mapping the cell are the identity and its own symmetry
|
|
// equivalents - collapses to a single candidate and never reaches the scoring at all.
|
|
//
|
|
// An alternative indexing maps the cell onto itself too, but it is not a description of the
|
|
// lattice for the model to be moved into. Where the model and the data share a point group, it
|
|
// is left out here and the DATA are relabelled into the model's indexing by the probe further
|
|
// down instead - so that every dataset of one crystal form comes out in the model's convention,
|
|
// and their maps and reflection files can be compared directly. Moving the model would leave
|
|
// each dataset in whichever indexing its own run happened to pick.
|
|
if (data_cell.is_crystal()) {
|
|
const gemmi::GroupOps model_gops = sg->operations();
|
|
const bool same_point_group =
|
|
data_space_group != nullptr && same_rotations(data_space_group->operations(), model_gops);
|
|
std::vector<gemmi::Op> frames{gemmi::Op::identity()};
|
|
for (const gemmi::Op &op : CellMappingOperators(st.cell, data_cell)) {
|
|
bool seen = false;
|
|
for (const gemmi::Op &kept : frames)
|
|
seen = seen || same_frame(kept, op, model_gops);
|
|
if (same_point_group)
|
|
for (const gemmi::Op &law : reindex_ops)
|
|
seen = seen || same_frame(law, op, model_gops);
|
|
if (!seen)
|
|
frames.push_back(op);
|
|
}
|
|
if (frames.size() > 1) {
|
|
logger.Info("Model validation: the model's cell {:.2f} {:.2f} {:.2f} {:.1f} {:.1f} {:.1f} ({}) "
|
|
"is not how the data describe this lattice; scoring {} change(s) of basis",
|
|
st.cell.a, st.cell.b, st.cell.c, st.cell.alpha, st.cell.beta, st.cell.gamma,
|
|
sg->xhm(), frames.size() - 1);
|
|
size_t best = 0;
|
|
double best_r = 0;
|
|
for (size_t i = 0; i < frames.size(); i++) {
|
|
gemmi::Structure trial = st;
|
|
const gemmi::SpaceGroup *trial_sg = sg;
|
|
if (i > 0 && !change_model_basis(trial, trial_sg, frames[i]))
|
|
continue; // no origin shift names the transformed group; not a basis we can take
|
|
refractionalize_into(trial, data_cell);
|
|
trial.setup_cell_images();
|
|
// Each frame at the data's best indexing: the alternative indexings are only probed
|
|
// further down, so a frame scored against the data in the wrong one reads as random
|
|
// as every other and the ranking is noise. Measured on a pseudo-orthorhombic P2
|
|
// crystal: every frame read R 0.77-0.81 in the merged indexing, and the frame chosen
|
|
// on that noise swapped two axes the model shares with the data.
|
|
double r = frame_probe_r(trial, trial_sg, obs, d_min);
|
|
if (same_point_group)
|
|
for (const gemmi::Op &law : reindex_ops)
|
|
r = std::min(r, frame_probe_r(trial, trial_sg, ReindexReflections(obs, law), d_min));
|
|
logger.Info("Model validation: {:<12} -> {:<12} R {:.4f} on the coarse shell",
|
|
frames[i].triplet(), trial_sg->xhm(), r);
|
|
if (i == 0 || r < best_r) {
|
|
best_r = r;
|
|
best = i;
|
|
}
|
|
}
|
|
if (best > 0) {
|
|
result.setting_as_read = sg->xhm();
|
|
if (change_model_basis(st, sg, frames[best])) {
|
|
result.change_of_basis_op = frames[best];
|
|
logger.Info("Model validation: model put through {} into {} - the data's own "
|
|
"description of the same lattice", frames[best].triplet(), sg->xhm());
|
|
} else {
|
|
result.setting_as_read.clear();
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Re-fractionalize the model into the data cell (rigid cell adjustment; no refinement).
|
|
if (data_cell.is_crystal())
|
|
refractionalize_into(st, 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 = MapToFPhi(dc.grid)
|
|
.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 = MapToFPhi(mask_grid).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.
|
|
// What one observed reflection contributes: the amplitude the R-factors and the maps are built
|
|
// from, the intensity CC(model, data) correlates, and the resolution that bins it.
|
|
struct Obs { double F; double I; float d; bool free; };
|
|
struct Fit {
|
|
gemmi::AsuData<std::complex<float>> fmodel;
|
|
gemmi::AsuData<gemmi::ValueSigma<float>> fobs_work; // what it was fitted to
|
|
std::unordered_map<long, Obs> obs_by_hkl;
|
|
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.
|
|
//
|
|
// Every reflection goes into the map that the R-factors and the maps are read from, but only
|
|
// the WORKING ones into the amplitudes the scale is fitted to. R-free is only worth quoting if
|
|
// no parameter the model was scaled by ever saw a free reflection, and the scale below is
|
|
// eleven of them - a scale, an anisotropic B and two bulk-solvent constants - fitted to
|
|
// minimise exactly the sum R is made of. Fitting them on all the data pulls Fmodel towards the
|
|
// free set as well and reports an R-free that is a little too good, by an amount nobody
|
|
// downstream can subtract off again. The same working set is what the rigid-body placement
|
|
// below is refined against.
|
|
gemmi::AsuData<gemmi::ValueSigma<float>> &fobs_work = out.fobs_work;
|
|
fobs_work.unit_cell_ = ucell;
|
|
fobs_work.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;
|
|
if (!r.rfree_flag)
|
|
fobs_work.v.push_back({h, {r.F, 1.0f}});
|
|
out.obs_by_hkl[hkl_key(h)] = {r.F, r.I, r.d, r.rfree_flag};
|
|
}
|
|
fobs_work.ensure_asu();
|
|
fobs_work.ensure_sorted();
|
|
|
|
// --- scale Fmodel(+solvent) to the WORKING Fobs: k_overall, anisotropic B, k_sol, b_sol ---
|
|
// Fitted on the working set, then applied to every reflection: scale_data walks the whole of
|
|
// fmodel, so the free reflections are put on the same scale as the rest without having had a
|
|
// say in what that scale is, and R-free below is computed against them.
|
|
gemmi::Scaling<float> scaling(ucell, sg);
|
|
scaling.use_solvent = true;
|
|
scaling.prepare_points(out.fmodel, fobs_work, &ms.fmask);
|
|
FitModelScale(scaling, {}, nthreads);
|
|
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.F;
|
|
double Fc = std::abs(hv.value);
|
|
if (it->second.free) { 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.
|
|
//
|
|
// `fitted`: fit_model(ms, obs_in) where the caller already has it - the indexing probe fits the
|
|
// identity labelling first, on the same structure factors and the same reflections - so the same
|
|
// fit is not made twice in a row.
|
|
//
|
|
// A rejected placement puts the atoms back but not ms.fcalc/ms.fmask, which are left at the
|
|
// rejected placement: nothing reads them afterwards - every later fit_model on this state is
|
|
// preceded by its own compute_model_factors - and recomputing them was a full Fcalc and mask
|
|
// for nothing.
|
|
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,
|
|
std::optional<Fit> fitted) -> Placement {
|
|
Placement out;
|
|
out.fit = fitted ? std::move(*fitted) : 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.
|
|
const std::vector<gemmi::Position> before = ModelPositions(ms.st.models[0]);
|
|
const RigidBodyRefineResult rb =
|
|
RefineRigidBody(ms.st.models[0], ucell, *sg, out.fit.fobs_work, d_min, logger, nthreads);
|
|
if (!rb.converged) {
|
|
// RefineRigidBody leaves the model wherever the solver left it, usable answer or not, so
|
|
// the restore cannot be conditional on the same flag the re-fit is.
|
|
SetModelPositions(ms.st.models[0], before);
|
|
return out;
|
|
}
|
|
{
|
|
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);
|
|
}
|
|
}
|
|
return out;
|
|
};
|
|
|
|
// --- indexing (merohedral) ambiguity: the probe over the candidates reindex_ops holds ---
|
|
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
|
|
std::optional<Fit> identity_fit; // fit_model(ms, obs), for place_and_fit to reuse
|
|
};
|
|
auto probe_indexing = [&](ModelState &ms) {
|
|
IndexingProbe out;
|
|
out.identity_fit = fit_model(ms, obs);
|
|
std::vector<double> r_free{out.identity_fit->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;
|
|
}
|
|
}
|
|
// The probe's identity fit is this fit exactly when the identity won: same structure factors
|
|
// (nothing has moved the model since), same reflections.
|
|
std::optional<Fit> identity_fit;
|
|
if (obs_model == &obs)
|
|
identity_fit = std::move(indexing.identity_fit);
|
|
indexing.identity_fit.reset();
|
|
Placement real = place_and_fit(mdl, *obs_model, std::move(identity_fit));
|
|
|
|
// --- is there anything for the model to decide? ---
|
|
// There are three: the space-group label, where the model asserts the other enantiomorph; the
|
|
// indexing, where a relabelling of the data fits it better; and the setting, where the model is
|
|
// written on other axes, which the files then take (Rugnux.cpp). A model that asserts none -
|
|
// already in the group and on the axes 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())
|
|
|| !(result.change_of_basis_op == gemmi::Op::identity());
|
|
// A null that cannot be built is a question that could not be put, and this design already has a
|
|
// state for that. Without this the run dies here - model validation runs BEFORE the reflection
|
|
// files are written, so one failed replicate would take the .mtz, .cif, .hkl and _unmerged.mtz
|
|
// with it, after the merge has already been paid for. The replicates are rotated models fed to a
|
|
// scaling path that throws on data it cannot pair up, so this is the adversarial input for it.
|
|
try {
|
|
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. Not because nothing is refined against it - the placement's six parameters
|
|
// are, and the scale's four - but because every null replicate is placed and scaled the same
|
|
// way, so whatever that optimism is worth is bought on both sides and cancels in
|
|
// (mean - real)/sd. R-work is then decided on an order of magnitude more reflections.
|
|
//
|
|
// 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.
|
|
//
|
|
// A replicate has to start from an orientation the model does NOT have, and there is more than
|
|
// one orientation to keep away from. The crystal looks the same from every orientation its
|
|
// lattice's rotations take the model to: the space group's own, and those composed with a twin
|
|
// law, which is where the model sits against data merged in the other indexing - and the
|
|
// replicates are fitted to the data as merged. A draw within the rigid body's reach of any of
|
|
// them is walked onto it and scores like the real model: measured on a cubic crystal with a
|
|
// twin law, one replicate of nine drawn 7.9 deg from the twin-related orientation was placed
|
|
// to R-work 0.20 against 0.57 for the other eight, and the spread that one gave the null was
|
|
// enough for the crystal's own model to read as not fitting. Such a draw is not a sample of
|
|
// the null, so it is drawn again. Everywhere else the draws are exactly the ones taken before.
|
|
std::vector<gemmi::Op> lattice_ops{gemmi::Op::identity()};
|
|
for (const gemmi::Op &law : ReindexAmbiguityOperators(cell, ambiguity_sg))
|
|
lattice_ops.push_back(law);
|
|
std::vector<gemmi::Mat33> equivalent;
|
|
for (const gemmi::Op &law : lattice_ops)
|
|
for (const gemmi::Op &s : gops.sym_ops)
|
|
equivalent.push_back(cartesian_rotation(s.combine(law), ucell));
|
|
const double reach_deg = RigidBodyReachDeg(st.models[0], d_min);
|
|
std::vector<gemmi::Mat33> null_rotation;
|
|
{
|
|
std::mt19937 rng(NULL_SEED);
|
|
int redraws = 0;
|
|
while (static_cast<int>(null_rotation.size()) < NULL_REPLICATES) {
|
|
const gemmi::Mat33 rot = random_rotation(rng);
|
|
if (angle_to_nearest_deg(rot, equivalent) < reach_deg && redraws < NULL_MAX_REDRAWS) {
|
|
redraws++;
|
|
continue;
|
|
}
|
|
null_rotation.push_back(rot);
|
|
}
|
|
if (redraws > 0)
|
|
logger.Info("Model validation: {} random placement(s) fell within the rigid body's reach "
|
|
"({:.1f} deg) of an orientation equivalent to the model's, and were drawn again",
|
|
redraws, reach_deg);
|
|
}
|
|
// 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);
|
|
auto run_replicate = [&](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);
|
|
std::optional<Fit> identity_fit;
|
|
if (!reindex_ops.empty()) {
|
|
IndexingProbe probe = probe_indexing(rep);
|
|
null_margin[i] = probe.margin;
|
|
identity_fit = std::move(probe.identity_fit); // replicates are placed against obs
|
|
}
|
|
null_r_work[i] = place_and_fit(rep, obs, std::move(identity_fit)).fit.r_work;
|
|
};
|
|
// On threads of their own, not on ParallelFor's pool. A parallel pass reached from a pool
|
|
// worker runs inline (ParallelFor.h), so there the scale fits and the rigid-body Jacobian of
|
|
// every replicate but the one on the calling thread ran serially, and the null lasted as long
|
|
// as its slowest serial replicate. From threads outside the pool those passes spread over the
|
|
// pool as the real fit's do. Each pass splits its work the same way wherever it runs, so the
|
|
// numbers are the same.
|
|
std::atomic<int> next_replicate{0};
|
|
std::vector<std::future<void>> replicate_threads;
|
|
for (size_t t = 0; t < std::min<size_t>(std::max<size_t>(nthreads, 1), NULL_REPLICATES); t++)
|
|
replicate_threads.push_back(std::async(std::launch::async, [&] {
|
|
for (int i = next_replicate++; i < NULL_REPLICATES; i = next_replicate++)
|
|
run_replicate(i);
|
|
}));
|
|
for (std::future<void> &f : replicate_threads)
|
|
f.get();
|
|
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) {
|
|
// The spread is floored, not trusted down to zero: sd == 0 and sd == 1e-7 are one ulp apart
|
|
// and would give sigmas of nothing and of millions. A null with no spread - a model too
|
|
// small or too symmetric for a rotation about its own centroid to move |Fcalc| - fits no
|
|
// better than its own random placements, so it still reads near zero here. The floor is
|
|
// not a gate: a large, well-measured crystal has a null as narrow as 0.0006-0.0010
|
|
// (measured), and a real fit 0.3 below it is hundreds of sigma, not none.
|
|
constexpr double NULL_SD_FLOOR = 1e-3;
|
|
result.r_work_sigma = (null_mean - p.fit.r_work) / std::max(null_sd, NULL_SD_FLOOR);
|
|
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;
|
|
constexpr double MARGIN_SD_FLOOR = 1e-4;
|
|
result.indexing_margin_sigma = (indexing.margin - margin_mean) / std::max(margin_sd, MARGIN_SD_FLOOR);
|
|
// 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;
|
|
result.indexing_decided = decided;
|
|
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, std::nullopt);
|
|
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 asserts no other enantiomorph for these data and "
|
|
"prefers the indexing they were merged in, 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);
|
|
}
|
|
} catch (const std::exception &e) {
|
|
// The null could not be built, so no decision is licensed and none is taken - which is the
|
|
// NOT_TESTED state, not a failure of the run. If the indexing probe had already won on a
|
|
// relabelling, the fit that survives describes those reflections and the files will not carry
|
|
// them, so it is remade on the data as merged.
|
|
logger.Warning("Model validation: the null could not be built ({}), so the model decides "
|
|
"nothing; the reflections keep the group and the indexing they were merged in",
|
|
e.what());
|
|
result.fit_tested = false;
|
|
result.model_fits = false;
|
|
result.indexing_probed = false;
|
|
result.indexing_decided = false;
|
|
result.indexing_op = gemmi::Op::identity();
|
|
result.r_work_sigma = 0.0;
|
|
SetModelPositions(st.models[0], as_read);
|
|
compute_model_factors(mdl);
|
|
real = place_and_fit(mdl, obs, std::nullopt);
|
|
}
|
|
|
|
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, Obs> &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;
|
|
|
|
// --- CC(model, data) by resolution shell ---
|
|
// The correlation of the merged intensities with |F_model|^2, binned on the merge table's own
|
|
// shells so the two tables line up row for row. Nearly free: the scaled Fmodel is already here,
|
|
// fitted with the eleven parameters above and nothing more, and this is the statistic it is best
|
|
// suited to. See the note on cc_model_shells in ModelValidation.h for what it can and cannot
|
|
// decide - it is a one-sided test, and it only ever argues for MORE resolution.
|
|
if (!report_shell_d_min.empty()) {
|
|
// Shells run coarse to fine and each is labelled by the resolution it reaches, so a reflection
|
|
// belongs to the first shell whose bound it has not passed.
|
|
auto shell_of = [&](float d) {
|
|
for (size_t i = 0; i < report_shell_d_min.size(); i++)
|
|
if (d > report_shell_d_min[i])
|
|
return i;
|
|
return report_shell_d_min.size();
|
|
};
|
|
std::vector<CorrelationCoefficient> shell_cc(report_shell_d_min.size());
|
|
std::vector<int> shell_n(report_shell_d_min.size(), 0);
|
|
// Sums for the radial mismatch: the shell's mean observed intensity against the mean
|
|
// intensity of the already-scaled model. See ModelDataCCShell::i_over_model.
|
|
std::vector<double> shell_sum_I(report_shell_d_min.size(), 0.0);
|
|
std::vector<double> shell_sum_Im(report_shell_d_min.size(), 0.0);
|
|
CorrelationCoefficient overall_cc;
|
|
int overall_n = 0;
|
|
for (const auto &hv : fmodel.v) {
|
|
const auto it = obs_by_hkl.find(hkl_key(hv.hkl));
|
|
if (it == obs_by_hkl.end() || !std::isfinite(it->second.I))
|
|
continue;
|
|
const size_t bin = shell_of(it->second.d);
|
|
if (bin >= shell_cc.size()) // finer than the finest shell the merge reported
|
|
continue;
|
|
const double Ic = std::norm(hv.value); // |F_model|^2
|
|
shell_cc[bin].Add(it->second.I, Ic);
|
|
shell_sum_I[bin] += it->second.I;
|
|
shell_sum_Im[bin] += Ic;
|
|
++shell_n[bin];
|
|
overall_cc.Add(it->second.I, Ic);
|
|
++overall_n;
|
|
}
|
|
// Fisher's transform against a null of zero correlation. Below four reflections there is no
|
|
// score to give, and a correlation of exactly +-1 has no finite one.
|
|
auto fisher_sigma = [](double cc, int n) {
|
|
return (n > 3 && std::fabs(cc) < 1.0) ? std::atanh(cc) * std::sqrt(n - 3.0) : NAN;
|
|
};
|
|
result.cc_model_shells.reserve(report_shell_d_min.size());
|
|
for (size_t i = 0; i < report_shell_d_min.size(); i++) {
|
|
const double cc = shell_cc[i].GetCC();
|
|
result.cc_model_shells.push_back({report_shell_d_min[i], cc, shell_n[i],
|
|
fisher_sigma(cc, shell_n[i]),
|
|
shell_sum_Im[i] > 0.0 ? shell_sum_I[i] / shell_sum_Im[i] : NAN});
|
|
}
|
|
result.cc_model_overall = overall_cc.GetCC();
|
|
result.cc_model_n = overall_n;
|
|
logger.Info("Model validation: CC(model,data) overall {:.3f} on {} reflections; "
|
|
"outermost shell {:.2f} A: {:.3f} on {} ({:+.1f} sigma)",
|
|
result.cc_model_overall, result.cc_model_n,
|
|
result.cc_model_shells.back().d_min, result.cc_model_shells.back().cc,
|
|
result.cc_model_shells.back().n, result.cc_model_shells.back().sigma);
|
|
|
|
// --- the same R with the radial profile taken out, and how much of it there was ---
|
|
// R above is what the maps are made of and is read WITHIN a run; this is the one to read
|
|
// BETWEEN two. The scale the maps use can only bend as k*exp(-B s^2) (see the note in
|
|
// fit_model, and the reason it stays that way), so where two reductions of one crystal differ
|
|
// in the radial profile of their amplitudes, the part of the difference that shape cannot
|
|
// follow lands in R - measured on this corpus, enough to move R by more than a real change in
|
|
// the data does. Giving the scale one free factor per shell removes the radial profile and
|
|
// only the radial profile: what is left is the agreement inside each shell, which is the
|
|
// question "do these data fit this model better or worse than those did".
|
|
//
|
|
// On ALL the reflections, not the free set. Nothing here is refined - the coordinates, the B
|
|
// factors and the occupancies are the model's own - so work and free estimate the same
|
|
// quantity (measured: the two differ by a median 0.001 over the corpus), and the 5% split
|
|
// therefore buys no cross-validation while costing a factor of sqrt(20) in precision. Read on
|
|
// all of them, the number is worth about +-0.002 per set instead of +-0.010, which is the
|
|
// difference between seeing a real change and seeing the sampling of a free set.
|
|
//
|
|
// Nothing is written from any of this. fmodel, the maps, the MTZ and the R-factors above are
|
|
// untouched; it is a second reading of the fit that has already happened.
|
|
double num = 0, den = 0, plain_num = 0, plain_den = 0;
|
|
int n_scored = 0;
|
|
std::vector<double> sum_oc(report_shell_d_min.size(), 0), sum_cc(report_shell_d_min.size(), 0);
|
|
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 size_t bin = shell_of(it->second.d);
|
|
if (bin >= sum_oc.size()) continue;
|
|
sum_oc[bin] += it->second.F * std::abs(hv.value);
|
|
sum_cc[bin] += static_cast<double>(std::abs(hv.value)) * std::abs(hv.value);
|
|
}
|
|
// The shell's scale by least squares on its own reflections. One parameter per shell against
|
|
// the thousands of reflections in it, so scoring the same reflections it was fitted on lowers
|
|
// R by far less than the difference it is there to expose.
|
|
std::vector<double> ln_k;
|
|
ln_k.reserve(report_shell_d_min.size());
|
|
for (size_t i = 0; i < sum_oc.size(); i++)
|
|
if (sum_cc[i] > 0 && sum_oc[i] > 0)
|
|
ln_k.push_back(std::log(sum_oc[i] / sum_cc[i]));
|
|
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 size_t bin = shell_of(it->second.d);
|
|
if (bin >= sum_oc.size()) continue;
|
|
const double Fo = it->second.F, Fc = std::abs(hv.value);
|
|
const double k = (sum_cc[bin] > 0 && sum_oc[bin] > 0) ? sum_oc[bin] / sum_cc[bin] : 1.0;
|
|
num += std::fabs(Fo - k * Fc);
|
|
den += Fo;
|
|
plain_num += std::fabs(Fo - Fc);
|
|
plain_den += Fo;
|
|
++n_scored;
|
|
}
|
|
if (den > 0 && ln_k.size() >= 2) {
|
|
result.r_model = plain_num / plain_den;
|
|
result.r_model_shell_scaled = num / den;
|
|
result.n_model = n_scored;
|
|
result.shell_scale_bins = static_cast<int>(ln_k.size());
|
|
const double mean_ln_k = std::accumulate(ln_k.begin(), ln_k.end(), 0.0) / ln_k.size();
|
|
double var = 0;
|
|
for (double v : ln_k) var += (v - mean_ln_k) * (v - mean_ln_k);
|
|
result.radial_misfit = std::sqrt(var / ln_k.size());
|
|
logger.Info("Model validation: R(all reflections) {:.4f}, with a per-shell scale {:.4f} "
|
|
"over {} shells; radial misfit {:.3f}",
|
|
result.r_model, result.r_model_shell_scaled, result.shell_scale_bins,
|
|
result.radial_misfit);
|
|
}
|
|
// The radial mismatch, reported on its own because nothing else in the run can see it.
|
|
const double r_out = result.cc_model_shells.back().i_over_model;
|
|
if (std::isfinite(r_out))
|
|
logger.Info("Model validation: <I>/<I_model> after the model's own scale, anisotropic B and "
|
|
"bulk solvent runs to {:.2f} in the outermost shell ({:.2f} A) - the part of "
|
|
"this run's resolution dependence the model does not agree with, and the only "
|
|
"check on it that is not blind to a factor shared by a whole ASU group",
|
|
r_out, result.cc_model_shells.back().d_min);
|
|
}
|
|
|
|
// --- 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.F;
|
|
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.free, centric});
|
|
sa_input.push_back({Fo, Fc, ucell.calculate_1_d2(hv.hkl), gops.epsilon_factor(hv.hkl),
|
|
centric, it->second.free});
|
|
}
|
|
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
|
|
// No data group means the caller merged in P1, and P1 is what the
|
|
// reflection files beside this one carry - not the model's group.
|
|
: (data_space_group != nullptr ? data_space_group
|
|
: gemmi::find_spacegroup_by_number(1));
|
|
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++) {
|
|
// [0, 360), the convention the MTZ 'P' columns carried before sigma_A weighting: std::arg
|
|
// returns (-180, 180] and a P column is not supposed to.
|
|
const double phi_wrapped = terms[i].phi * 180.0 / PI;
|
|
const auto phi_deg = static_cast<float>(phi_wrapped < 0 ? phi_wrapped + 360.0 : phi_wrapped);
|
|
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, fitted on the working set only]",
|
|
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;
|
|
}
|
|
|
|
void KeepModelVerdict(ModelValidationResult &remade, const ModelValidationResult &decided) {
|
|
remade.fit_tested = decided.fit_tested;
|
|
remade.model_fits = decided.model_fits;
|
|
remade.null_r_work_mean = decided.null_r_work_mean;
|
|
remade.null_r_work_sd = decided.null_r_work_sd;
|
|
remade.r_work_sigma = decided.r_work_sigma;
|
|
remade.null_replicates = decided.null_replicates;
|
|
remade.indexing_op = decided.indexing_op;
|
|
remade.indexing_probed = decided.indexing_probed;
|
|
remade.indexing_decided = decided.indexing_decided;
|
|
remade.indexing_margin = decided.indexing_margin;
|
|
remade.indexing_margin_null_mean = decided.indexing_margin_null_mean;
|
|
remade.indexing_margin_null_sd = decided.indexing_margin_null_sd;
|
|
remade.indexing_margin_sigma = decided.indexing_margin_sigma;
|
|
remade.model_enantiomorph_candidate = decided.model_enantiomorph_candidate;
|
|
remade.adopted_model_enantiomorph = decided.adopted_model_enantiomorph;
|
|
remade.model_space_group_number = decided.model_space_group_number;
|
|
}
|
|
|
|
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 =
|
|
MapToFPhi(dc.grid).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 =
|
|
MapToFPhi(mask_grid).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;
|
|
}
|