leonarski_fandClaude Opus 5.5
672e182d6a
A pooled CudaDevicePtr frees on the thread's allocation stream, not on the engine's stream, and the
pool may hand the memory to another engine - or, past its release threshold, unmap it - as soon as
that free is reached, which on an idle allocation stream is at once. An engine destroyed with work
still queued (FFTIndexerGPU after SearchCap's last DirectionsChanged upload, a spot finder between
DetectAt and Extract, a shadow accumulator after a pending fold, any engine on an exception path)
thus had kernels or copies writing memory that was someone else's or no longer mapped. Now:
- CudaStream synchronises before cudaStreamDestroy (destructor and move-assignment), which covers
engines whose own stream is declared after their buffers (FFTIndexerGPU, the gather buffer);
- every engine holding pooled buffers and a stream (shared or own, declared before the buffers)
synchronises it in its destructor; BraggIntegrationEngineGPU also before EnsureCapacity
reallocates, where a Run that threw leaves work queued.
A GPU failure that is handled no longer hides a lost context: ShadowFinder, BeamCenterFFT, the
rigid-body pool and model validation call cuda_throw_if_context_lost() before cuda_clear_error(),
as the device-decode fallbacks already did; RotationScaleMergeGPU's Alloc does so before waiting up
to ten minutes for GPU work beside it and then reporting a lost device as out of memory; and a
failed cudaMalloc says why. BeamCenterFFT logs the failure it used to drop silently, the
speculative geometry probe logs the exception it swallowed (its GPU fault was otherwise reported by
the merge beside it, under the merge's name), and RotationScaleMergeGPU's DeviceGuard no longer
throws from its destructor.
Only synchronisation and error paths change: p.hkl md5 and the MTZ data (gemmi) are identical to
the b530c2d full battery on myob_x10sa, cytc_x10sa, 8a1a, 9gdj, 11if, kdp_x10sa_20keV and 6z9g.
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
1902 lines
107 KiB
C++
1902 lines
107 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 <cstdio>
|
|
#include <array>
|
|
#include <atomic>
|
|
#include <functional>
|
|
#include <memory>
|
|
#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 "../../common/ParallelFor.h"
|
|
#include "../scale_merge/ReindexAmbiguity.h" // ReindexReflections
|
|
#include "../scale_merge/CrystalSetting.h" // CellMappingOperators
|
|
#include "ModelFFT.h"
|
|
#include "RigidBodyRefine.h"
|
|
#include "SigmaA.h"
|
|
#ifdef JFJOCH_USE_CUDA
|
|
#include "ModelStructureFactorsGPU.h"
|
|
#include "RigidBodyGPU.h"
|
|
#include "RigidBodyGPUEngine.h" // RigidBodyGPUEngine::CurrentDevice
|
|
#include "../../common/CUDAWrapper.h"
|
|
#include "../../common/JFJochException.h"
|
|
#else
|
|
class ModelStructureFactorsGPU;
|
|
#endif
|
|
|
|
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, on `gpu` where the validation has one.
|
|
gemmi::Grid<float> map_from_coefficients(gemmi::AsuData<std::complex<float>> &coef, ModelStructureFactorsGPU *gpu) {
|
|
#ifdef JFJOCH_USE_CUDA
|
|
if (gpu != nullptr) {
|
|
try {
|
|
return gpu->Map(coef);
|
|
} catch (const JFJochException &e) {
|
|
throw RigidBodyGPUFailure(e.what());
|
|
}
|
|
}
|
|
#else
|
|
(void) gpu;
|
|
#endif
|
|
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));
|
|
}
|
|
|
|
// A file this run may already have written - the validation in the model's setting writes over the
|
|
// first one's maps - is removed before it is written again, not truncated: XFS (and ext4) flush a file
|
|
// truncated and rewritten when it is closed, and three 150 MB maps forced out to a disk that way cost
|
|
// seconds where a fresh file costs nothing.
|
|
void remove_before_rewriting(const std::string &path) {
|
|
std::remove(path.c_str());
|
|
}
|
|
|
|
// 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);
|
|
remove_before_rewriting(path);
|
|
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;
|
|
|
|
// The lightest element counted as an anomalous scatterer: phosphorus. Below it (C, N, O, Na, Mg) f''
|
|
// is a fraction of sulfur's at any wavelength these data are taken at.
|
|
constexpr int ANOMALOUS_SCATTERER_MIN_Z = 15;
|
|
|
|
// 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;
|
|
|
|
// The resolution the null and the real fit's side of it are scored to (see the null below), and the
|
|
// fewest free reflections it may leave the indexing margin, which is read in R-free.
|
|
constexpr double NULL_D_MIN = 3.5;
|
|
constexpr size_t NULL_MIN_FREE = 1000;
|
|
|
|
// 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();
|
|
}
|
|
|
|
#ifdef JFJOCH_USE_CUDA
|
|
// The model's structure factors to the data's resolution, and the maps, on the GPU - or not, decided here,
|
|
// once, before any of them is computed, and never revisited: nothing is moved to the CPU part way through.
|
|
// Decided on what the GPU path needs against the card's TOTAL memory and not against what happens to be
|
|
// free at the moment, so the same input on the same machine always takes the same path, whatever else is
|
|
// running beside it. A validation plans at most half of its card: its structure-factor engines a quarter
|
|
// together - `reserved` is what one of them already holds there - and its rigid-body engines a quarter
|
|
// (RigidBodyGPUPool::Create). A second validation runs beside it only where twice that plan fits half of
|
|
// the cards together (see on_forecast below), so validations never plan more than half of the memory.
|
|
std::unique_ptr<ModelStructureFactorsGPU> StructureFactorsGPU(int device, const gemmi::Model &model,
|
|
const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
|
|
double d_min, size_t reserved, Logger &logger) {
|
|
std::unique_ptr<ModelStructureFactorsGPU> sf;
|
|
try {
|
|
sf = std::make_unique<ModelStructureFactorsGPU>(device, model, cell, sg, d_min);
|
|
} catch (const JFJochException &e) {
|
|
throw RigidBodyGPUFailure(e.what());
|
|
}
|
|
const std::array<int, 3> n = sf->GridSize();
|
|
if (!sf->Supported()) {
|
|
logger.Info("Model validation: structure factors to {:.2f} A on the CPU - the cell is too small for the "
|
|
"GPU's gridding", d_min);
|
|
return nullptr;
|
|
}
|
|
const size_t total = sf->DeviceTotalMemory();
|
|
if (reserved + sf->DeviceBytes() > total / 4) {
|
|
logger.Info("Model validation: structure factors to {:.2f} A on the CPU - on a {}x{}x{} grid they need "
|
|
"{:.2f} GB of GPU memory beside {:.2f} GB already reserved, over a quarter of the card's {:.2f} GB",
|
|
d_min, n[0], n[1], n[2], sf->DeviceBytes() / 1e9, reserved / 1e9, total / 1e9);
|
|
return nullptr;
|
|
}
|
|
try {
|
|
sf->Reserve();
|
|
} catch (const JFJochException &e) {
|
|
throw RigidBodyGPUFailure(e.what());
|
|
}
|
|
logger.Info("Model validation: structure factors to {:.2f} A on the GPU - a {}x{}x{} grid, {:.2f} GB of the "
|
|
"card's {:.2f} GB", d_min, n[0], n[1], n[2], sf->DeviceBytes() / 1e9, total / 1e9);
|
|
return sf;
|
|
}
|
|
#endif
|
|
|
|
} // namespace
|
|
|
|
namespace {
|
|
|
|
// ValidateAgainstModel() with the rigid body, the structure factors and the maps on the GPU where
|
|
// `rigid_body_gpu` allows it and a card is there to take them, on the CPU otherwise - for the whole
|
|
// validation either way.
|
|
ModelValidationResult Validate(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,
|
|
const ModelValidationSchedule &schedule,
|
|
const gemmi::Op *twin_law,
|
|
bool rigid_body_gpu) {
|
|
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);
|
|
// The candidates are scored at the same time, each on its own copy of the model, and read
|
|
// in their own order below, so the ranking and its tie-break are those of one after another.
|
|
std::vector<double> frame_r(frames.size(), 0.0);
|
|
std::vector<const gemmi::SpaceGroup *> frame_sg(frames.size(), nullptr);
|
|
ParallelFor(static_cast<int>(frames.size()), nthreads, [&](int i) {
|
|
gemmi::Structure trial = st;
|
|
const gemmi::SpaceGroup *trial_sg = sg;
|
|
if (i > 0 && !change_model_basis(trial, trial_sg, frames[i]))
|
|
return; // 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));
|
|
frame_r[i] = r;
|
|
frame_sg[i] = trial_sg;
|
|
});
|
|
size_t best = 0;
|
|
double best_r = 0;
|
|
for (size_t i = 0; i < frames.size(); i++) {
|
|
if (frame_sg[i] == nullptr)
|
|
continue;
|
|
const double r = frame_r[i];
|
|
logger.Info("Model validation: {:<12} -> {:<12} R {:.4f} on the coarse shell",
|
|
frames[i].triplet(), frame_sg[i]->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);
|
|
|
|
// The structure factors to d_min and the maps on the GPU, where `rigid_body_gpu` allows it and the
|
|
// card can take them; decided here, before any is computed. The null's, to its coarser limit, get an
|
|
// engine of their own further down, decided the same way before any replicate starts - for all of
|
|
// them and the real model's side of the null together, so that every one is computed the same way.
|
|
ModelStructureFactorsGPU *sf_gpu = nullptr, *sf_null = nullptr;
|
|
#ifdef JFJOCH_USE_CUDA
|
|
std::unique_ptr<ModelStructureFactorsGPU> sf_engine, sf_null_engine;
|
|
// On the current device, as the rigid body's engines take it: one card for the whole validation.
|
|
if (rigid_body_gpu && get_gpu_count() > 0)
|
|
sf_engine = StructureFactorsGPU(RigidBodyGPUEngine::CurrentDevice(), st.models[0], ucell, *sg, d_min, 0,
|
|
logger);
|
|
sf_gpu = sf_engine.get();
|
|
#endif
|
|
|
|
// --- 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, double to_d) {
|
|
#ifdef JFJOCH_USE_CUDA
|
|
ModelStructureFactorsGPU *gpu = sf_gpu != nullptr && to_d == sf_gpu->DMin() ? sf_gpu
|
|
: sf_null != nullptr && to_d == sf_null->DMin() ? sf_null
|
|
: nullptr;
|
|
if (gpu != nullptr) {
|
|
try {
|
|
gpu->Compute(ms.st.models[0], ms.fcalc, ms.fmask);
|
|
} catch (const JFJochException &e) {
|
|
throw RigidBodyGPUFailure(e.what());
|
|
}
|
|
return;
|
|
}
|
|
#endif
|
|
gemmi::DensityCalculator<Table, float> dc;
|
|
dc.d_min = to_d;
|
|
dc.rate = 1.5;
|
|
dc.set_grid_cell_and_spacegroup(ms.st);
|
|
dc.set_refmac_compatible_blur(ms.st.models[0]);
|
|
|
|
// The mask is made on a thread of its own, beside the density: both only read the model, each
|
|
// fills a grid of its own, and between them they are the whole cost of this step.
|
|
// Refmac radii give a slightly lower R than the Cctbx set on our test cases, at no cost.
|
|
std::future<gemmi::AsuData<std::complex<float>>> fmask = std::async(std::launch::async,
|
|
[&model = ms.st.models[0], unit_cell = dc.grid.unit_cell, spacegroup = dc.grid.spacegroup,
|
|
spacing = dc.requested_grid_spacing(), to_d] {
|
|
gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac);
|
|
gemmi::Grid<float> mask_grid;
|
|
mask_grid.unit_cell = unit_cell;
|
|
mask_grid.spacegroup = spacegroup;
|
|
mask_grid.set_size_from_spacing(spacing, gemmi::GridSizeRounding::Up);
|
|
masker.put_mask_on_grid(mask_grid, model);
|
|
return MapToFPhi(mask_grid).prepare_asu_data(to_d, 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);
|
|
ms.fmask = fmask.get();
|
|
};
|
|
// The model's structure factors are made beside the reservation of the GPU engines below, which
|
|
// does not need them.
|
|
std::future<void> model_factors = std::async(std::launch::async, [&] { compute_model_factors(mdl, d_min); });
|
|
|
|
// The rigid body's device, decided once for the whole validation - the real fit and every replicate
|
|
// of the null on the same one, so that a verdict never depends on how full the card was when a
|
|
// replicate started. The engines are reserved here, up front, and never grow.
|
|
RigidBodyGPUPool *rigid_body_pool = nullptr;
|
|
#ifdef JFJOCH_USE_CUDA
|
|
std::unique_ptr<RigidBodyGPUPool> rigid_body_engines;
|
|
if (rigid_body_gpu) {
|
|
const double finest_zone = RigidBodyLadder(d_min).back();
|
|
size_t zone_observations = 0;
|
|
for (const MergedReflection &r : obs)
|
|
if (!std::isnan(r.F) && r.d >= finest_zone)
|
|
++zone_observations;
|
|
// One engine for the real model and one for each replicate of the null placed beside it.
|
|
rigid_body_engines = RigidBodyGPUPool::Create(st.models[0], ucell, *sg, d_min, zone_observations,
|
|
std::min<size_t>(std::max<size_t>(nthreads, 1), NULL_REPLICATES) + 1,
|
|
logger);
|
|
rigid_body_pool = rigid_body_engines.get();
|
|
}
|
|
#else
|
|
(void) rigid_body_gpu;
|
|
#endif
|
|
model_factors.get();
|
|
|
|
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;
|
|
gemmi::Mat33 rb_rotation; // the committed rotation about the centroid, identity if none
|
|
};
|
|
auto place_and_fit = [&](ModelState &ms, const std::vector<MergedReflection> &obs_in,
|
|
std::optional<Fit> fitted, double to_d) -> 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, to_d, logger, nthreads, rigid_body_pool);
|
|
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, to_d);
|
|
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_rotation = rb.rotation;
|
|
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, const std::vector<MergedReflection> &obs_in) {
|
|
IndexingProbe out;
|
|
// Each relabelling is fitted on a thread of its own: the fits only read the structure factors.
|
|
std::vector<std::future<double>> relabelled;
|
|
for (const auto &op : reindex_ops)
|
|
relabelled.push_back(std::async(std::launch::async, [&, op] {
|
|
return fit_model(ms, ReindexReflections(obs_in, op)).r_free;
|
|
}));
|
|
out.identity_fit = fit_model(ms, obs_in);
|
|
std::vector<double> r_free{out.identity_fit->r_free}; // identity first, then the twin laws
|
|
double best_r_free = r_free.front();
|
|
for (size_t i = 0; i < reindex_ops.size(); i++) {
|
|
const gemmi::Op &op = reindex_ops[i];
|
|
const double cand = relabelled[i].get();
|
|
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, obs);
|
|
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();
|
|
|
|
// --- 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());
|
|
if (schedule.on_forecast) {
|
|
// Two validations at once take twice the engines, and are allowed where twice what this one plans -
|
|
// its rigid-body engines, its structure-factor engine and, where a null is coming, the null's,
|
|
// which is no larger - fits half of the cards' memory taken together (see StructureFactorsGPU).
|
|
bool may_run_beside = true;
|
|
#ifdef JFJOCH_USE_CUDA
|
|
size_t planned = 0, card = 0;
|
|
if (rigid_body_pool != nullptr) {
|
|
planned += rigid_body_pool->PlannedBytes();
|
|
card = rigid_body_pool->CardBytes();
|
|
}
|
|
if (sf_gpu != nullptr) {
|
|
planned += (decision_pending ? 2 : 1) * sf_gpu->DeviceBytes();
|
|
card = sf_gpu->DeviceTotalMemory();
|
|
}
|
|
if (planned > 0)
|
|
may_run_beside = 2 * planned <= static_cast<size_t>(get_gpu_count()) * card / 2;
|
|
#endif
|
|
schedule.on_forecast({result.change_of_basis_op, indexing.op, result.model_enantiomorph_candidate,
|
|
may_run_beside});
|
|
}
|
|
std::vector<gemmi::Mat33> equivalent; // the orientations the crystal cannot tell from the model's
|
|
double reach_deg = 0;
|
|
std::vector<gemmi::Mat33> null_rotation;
|
|
double null_d_min = d_min;
|
|
std::vector<MergedReflection> obs_null;
|
|
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);
|
|
for (const gemmi::Op &law : lattice_ops)
|
|
for (const gemmi::Op &s : gops.sym_ops)
|
|
equivalent.push_back(cartesian_rotation(s.combine(law), ucell));
|
|
reach_deg = RigidBodyReachDeg(st.models[0], d_min);
|
|
{
|
|
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);
|
|
}
|
|
// The null and the side of the real fit it is compared with are both scored to NULL_D_MIN,
|
|
// not to the data's resolution. The placement already stops there (RigidBodyRefine's ladder
|
|
// ends at 3.5 A), and whether a model is in the right orientation, or the data in its
|
|
// indexing, is decided at low resolution: the finer shells only add reflections to the scale
|
|
// fits, which were most of the null's cost. Both sides go through the same fits on the same
|
|
// reflections, so the comparison stays like for like. On a small cell 3.5 A leaves too few
|
|
// free reflections for a real indexing margin to be told from a random placement's, so the
|
|
// limit then moves out to where NULL_MIN_FREE of them are - to the data's own on the
|
|
// smallest, which are the cheap ones.
|
|
std::vector<float> free_d;
|
|
for (const MergedReflection &r : obs)
|
|
if (r.rfree_flag && !std::isnan(r.F))
|
|
free_d.push_back(r.d);
|
|
std::sort(free_d.begin(), free_d.end(), std::greater<>());
|
|
null_d_min = std::max(d_min, NULL_D_MIN);
|
|
if (free_d.size() < NULL_MIN_FREE)
|
|
null_d_min = d_min;
|
|
else
|
|
null_d_min = std::min<double>(null_d_min, free_d[NULL_MIN_FREE - 1]);
|
|
for (const MergedReflection &r : obs)
|
|
if (r.d >= null_d_min)
|
|
obs_null.push_back(r);
|
|
#ifdef JFJOCH_USE_CUDA
|
|
// On the card the d_min engine is on; the replicates' threads may sit on other cards, and each
|
|
// call works on the engine's own device. Where the null is scored to d_min it shares that engine.
|
|
if (sf_gpu != nullptr && null_d_min != d_min)
|
|
sf_null_engine = StructureFactorsGPU(sf_gpu->Device(), st.models[0], ucell, *sg, null_d_min,
|
|
sf_gpu->DeviceBytes(), logger);
|
|
sf_null = sf_null_engine.get();
|
|
#endif
|
|
}
|
|
const gemmi::Structure null_model = st; // the model as read, which the replicates are made from
|
|
// Each replicate works on its own copy of the model as read and of its structure factors, and the
|
|
// replicates share nothing else, so they run concurrently, and beside the real model's placement
|
|
// below, which moves the real model - the one the maps and the atom-density readout are still
|
|
// taken from - and nothing a replicate reads.
|
|
//
|
|
// The check before the fit is only a filter: the rigid body can walk further than its nominal
|
|
// reach, and a replicate it carries onto an orientation equivalent to the model's has been
|
|
// refined into the model's own solution (measured: 24 deg, R-free 0.56 -> 0.39, which took the
|
|
// null's spread from 0.016 to 0.066 and a fitting model below the gate). So where the
|
|
// replicate ENDS is checked as well, and one that ended within reach is drawn again and
|
|
// placed again - from a generator of its own, so the draws a replicate takes do not depend on
|
|
// the order the concurrent replicates finish in.
|
|
std::vector<double> null_r_work(NULL_REPLICATES), null_margin(NULL_REPLICATES);
|
|
std::atomic<int> placed_redraws{0};
|
|
auto run_replicate = [&](int i) {
|
|
std::mt19937 rng(NULL_SEED + 1 + i);
|
|
gemmi::Mat33 rot = null_rotation[i];
|
|
for (int redraw = 0;; redraw++) {
|
|
ModelState rep;
|
|
rep.st = null_model;
|
|
reorient_about_centroid(rep.st.models[0], rot);
|
|
compute_model_factors(rep, null_d_min);
|
|
std::optional<Fit> identity_fit;
|
|
if (!reindex_ops.empty()) {
|
|
IndexingProbe probe = probe_indexing(rep, obs_null);
|
|
null_margin[i] = probe.margin;
|
|
identity_fit = std::move(probe.identity_fit); // replicates are placed against obs
|
|
}
|
|
const Placement placed = place_and_fit(rep, obs_null, std::move(identity_fit), null_d_min);
|
|
null_r_work[i] = placed.fit.r_work;
|
|
const gemmi::Mat33 ended = placed.rb_rotation.multiply(rot);
|
|
if (redraw >= NULL_MAX_REDRAWS || angle_to_nearest_deg(ended, equivalent) >= reach_deg)
|
|
break;
|
|
++placed_redraws;
|
|
rot = random_rotation(rng);
|
|
for (int k = 0; k < NULL_MAX_REDRAWS && angle_to_nearest_deg(rot, equivalent) < reach_deg; k++)
|
|
rot = random_rotation(rng);
|
|
}
|
|
};
|
|
// 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;
|
|
if (decision_pending)
|
|
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, [&, t] {
|
|
// Each on a card in turn from the one after the real fit's, which its engines are
|
|
// spread the same way over (RigidBodyGPUPool::Create).
|
|
#ifdef JFJOCH_USE_CUDA
|
|
if (rigid_body_pool != nullptr)
|
|
pin_gpu((rigid_body_pool->Device() + 1 + static_cast<int>(t)) % get_gpu_count());
|
|
#else
|
|
(void) t;
|
|
#endif
|
|
for (int i = next_replicate++; i < NULL_REPLICATES; i = next_replicate++)
|
|
run_replicate(i);
|
|
}));
|
|
// The real model is placed while they run; what is compared with them is read off it below.
|
|
Placement real = place_and_fit(mdl, *obs_model, std::move(identity_fit), d_min);
|
|
|
|
// 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 real model's side of it: the indexing margin from the model as read, as a replicate's
|
|
// is, and the R-work of the model where its placement left it, against the reflections that
|
|
// placement was fitted to.
|
|
auto coarse_r_work = [&](const std::vector<MergedReflection> &placed_against) {
|
|
ModelState coarse;
|
|
coarse.st = mdl.st;
|
|
compute_model_factors(coarse, null_d_min);
|
|
std::vector<MergedReflection> in;
|
|
for (const MergedReflection &r : placed_against)
|
|
if (r.d >= null_d_min)
|
|
in.push_back(r);
|
|
return fit_model(coarse, in).r_work;
|
|
};
|
|
double real_margin = 0;
|
|
if (!reindex_ops.empty()) {
|
|
ModelState coarse;
|
|
coarse.st = mdl.st;
|
|
SetModelPositions(coarse.st.models[0], as_read);
|
|
compute_model_factors(coarse, null_d_min);
|
|
real_margin = probe_indexing(coarse, obs_null).margin;
|
|
}
|
|
double real_r_work = coarse_r_work(*obs_model);
|
|
|
|
for (std::future<void> &f : replicate_threads)
|
|
f.get();
|
|
if (placed_redraws > 0)
|
|
logger.Info("Model validation: {} random placement(s) were refined to within the rigid body's reach "
|
|
"({:.1f} deg) of an orientation equivalent to the model's, and were drawn and placed again",
|
|
placed_redraws.load(), reach_deg);
|
|
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 = [&](double r_work) {
|
|
// 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 - r_work) / std::max(null_sd, NULL_SD_FLOOR);
|
|
result.model_fits = result.r_work_sigma >= MODEL_FIT_SIGMA;
|
|
};
|
|
record_fit(real_r_work);
|
|
|
|
// --- 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 = real_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 = (real_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 to {:.2f} A, against {:.4f} +- {:.4f} for a "
|
|
"random placement of the same model ({:+.2f} sigma) => {}",
|
|
reindex_ops.size() + 1, real_margin, null_d_min, 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, d_min);
|
|
real = place_and_fit(mdl, obs, std::nullopt, d_min);
|
|
real_r_work = coarse_r_work(obs);
|
|
record_fit(real_r_work);
|
|
}
|
|
}
|
|
logger.Info("Model validation: R-work {:.4f} to {:.2f} A against a null of {:.4f} +- {:.4f} = "
|
|
"{:+.2f} sigma => the model {} these data",
|
|
real_r_work, null_d_min, 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);
|
|
}
|
|
#ifdef JFJOCH_USE_CUDA
|
|
} catch (const RigidBodyGPUFailure &) {
|
|
throw; // not the null's failure but the device's: the whole validation runs again on the CPU
|
|
#endif
|
|
} 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, d_min);
|
|
real = place_and_fit(mdl, obs, std::nullopt, d_min);
|
|
}
|
|
|
|
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;
|
|
|
|
// --- the twin target, where a twin law was given ---
|
|
// The intensity a twinned crystal records at h is (1-a) I(h) + a I(Th), so the model is compared
|
|
// against the data as twinned: |F_twin(h)|^2 = (1-a)|F_model(h)|^2 + a|F_model(Th)|^2 from the
|
|
// placed, scaled model, with one overall scale refitted at each a. The fraction is scanned in steps
|
|
// of 0.01 and chosen on the working set; R-free is read at it. Nothing else is refitted against the
|
|
// twinned target - the placement, the anisotropic B and the bulk solvent are those of the untwinned
|
|
// fit - so this is the same few-parameter comparison as R above, not a twin refinement.
|
|
// Following Yeates (1997) Methods Enzymol. 276, 344-358
|
|
if (twin_law != nullptr) {
|
|
std::unordered_map<long, double> model_i;
|
|
model_i.reserve(fmodel.v.size());
|
|
for (const auto &hv : fmodel.v)
|
|
model_i[hkl_key(hv.hkl)] = std::norm(hv.value);
|
|
struct TwinTerm { double fo, i1, i2; bool free; };
|
|
std::vector<TwinTerm> twin_terms;
|
|
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 gemmi::Miller mate = asu.to_asu(twin_law->apply_to_hkl(hv.hkl), gops).first;
|
|
const auto m = model_i.find(hkl_key(mate));
|
|
if (m == model_i.end()) continue;
|
|
twin_terms.push_back({it->second.F, std::norm(hv.value), m->second, it->second.free});
|
|
}
|
|
double best_r_work = 2.0;
|
|
for (int step = 0; step <= 50 && !twin_terms.empty(); ++step) {
|
|
const double a = 0.01 * step;
|
|
double fo_fc = 0, fc_fc = 0;
|
|
for (const auto &t : twin_terms)
|
|
if (!t.free) {
|
|
const double fc = std::sqrt((1.0 - a) * t.i1 + a * t.i2);
|
|
fo_fc += t.fo * fc;
|
|
fc_fc += fc * fc;
|
|
}
|
|
const double k = fc_fc > 0 ? fo_fc / fc_fc : 1.0;
|
|
double num_w = 0, den_w = 0, num_f = 0, den_f = 0;
|
|
for (const auto &t : twin_terms) {
|
|
const double diff = std::fabs(t.fo - k * std::sqrt((1.0 - a) * t.i1 + a * t.i2));
|
|
if (t.free) { num_f += diff; den_f += t.fo; }
|
|
else { num_w += diff; den_w += t.fo; }
|
|
}
|
|
if (den_w > 0 && den_f > 0 && num_w / den_w < best_r_work) {
|
|
best_r_work = num_w / den_w;
|
|
result.r_work_twin = best_r_work;
|
|
result.r_free_twin = num_f / den_f;
|
|
result.twin_fraction = a;
|
|
}
|
|
}
|
|
result.twin_target = std::isfinite(result.r_free_twin);
|
|
if (result.twin_target)
|
|
logger.Info("Model validation: against the data as twinned by {} - fraction {:.2f} (scanned on the "
|
|
"working set) - R-work {:.4f}, R-free {:.4f}; untwinned {:.4f} / {:.4f}",
|
|
twin_law->as_hkl().triplet('h'), result.twin_fraction, result.r_work_twin,
|
|
result.r_free_twin, result.r_work, result.r_free);
|
|
}
|
|
|
|
// --- 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});
|
|
}
|
|
|
|
// A validation started before the one it follows had decided writes nothing until that one has,
|
|
// and nothing at all where it decided otherwise.
|
|
if (schedule.write_gate.valid() && !schedule.write_gate.get()) {
|
|
result.failure_reason = "superseded before its files were written";
|
|
return result;
|
|
}
|
|
// Where the validation in the model's setting, started beside this one, is going to be kept, it
|
|
// writes the maps and the map coefficients on the axes the files are written on, and this one's would
|
|
// be overwritten unread: they are not made at all.
|
|
const bool maps_superseded = schedule.maps_superseded && schedule.maps_superseded(result);
|
|
|
|
if (!maps_superseded) {
|
|
// --- write the maps and score the 2mFo-DFc map at atom centres (a real map peaks there) ---
|
|
// The three maps - this one, the difference map and the anomalous map - share nothing but what
|
|
// they read, so the other two are made and written on threads of their own beside this one.
|
|
std::future<void> fofc_map = std::async(std::launch::async, [&] {
|
|
write_ccp4(map_from_coefficients(mapfofc, sf_gpu), output_prefix + "_fofc.ccp4");
|
|
});
|
|
|
|
// --- 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
|
|
std::future<void> anomalous_map = std::async(std::launch::async, [&] {
|
|
// 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, sf_gpu);
|
|
const double rms = write_ccp4(grid, output_prefix + "_anom.ccp4");
|
|
std::vector<ModelValidationResult::AnomalousSite> sites;
|
|
double scatterer_sum = 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) {
|
|
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});
|
|
if (a.element.atomic_number() >= ANOMALOUS_SCATTERER_MIN_Z) {
|
|
scatterer_sum += sites.back().sigma;
|
|
result.anomalous_scatterers++;
|
|
}
|
|
}
|
|
if (result.anomalous_scatterers > 0)
|
|
result.anomalous_scatterer_mean_sigma = scatterer_sum / result.anomalous_scatterers;
|
|
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);
|
|
}
|
|
});
|
|
|
|
const gemmi::Grid<float> grid2fofc = map_from_coefficients(map2fofc, sf_gpu);
|
|
const double rms2 = write_ccp4(grid2fofc, output_prefix + "_2fofc.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;
|
|
}
|
|
fofc_map.get();
|
|
anomalous_map.get();
|
|
}
|
|
|
|
// --- 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");
|
|
|
|
if (!maps_superseded) {
|
|
// --- 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);
|
|
remove_before_rewriting(output_prefix + "_maps.mtz");
|
|
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);
|
|
if (!maps_superseded)
|
|
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);
|
|
if (maps_superseded)
|
|
logger.Info("Model validation: no maps made here - the validation in the model's setting writes them");
|
|
else
|
|
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
|
|
|
|
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,
|
|
const ModelValidationSchedule &schedule,
|
|
const gemmi::Op *twin_law) {
|
|
#ifdef JFJOCH_USE_CUDA
|
|
// A CUDA failure ends the validation, not the run, and nothing is moved to the CPU part way through:
|
|
// the validation reports why it did not finish, and a validation that did not finish decides nothing,
|
|
// so the reflection files are those of a run without a model.
|
|
try {
|
|
return Validate(merged, cell, model_path, output_prefix, logger, data_space_group,
|
|
probe_indexing_ambiguity, nthreads, wavelength_A, report_shell_d_min, schedule, twin_law,
|
|
true);
|
|
} catch (const RigidBodyGPUFailure &e) {
|
|
ModelValidationResult failed;
|
|
failed.model_path = model_path;
|
|
failed.failure_reason = fmt::format("the GPU failed during the validation ({})", e.what());
|
|
logger.Error("Model validation: {}", failed.failure_reason);
|
|
cuda_throw_if_context_lost(); // a lost device ends the run here, not at its next user
|
|
cuda_clear_error();
|
|
return failed;
|
|
}
|
|
#else
|
|
return Validate(merged, cell, model_path, output_prefix, logger, data_space_group, probe_indexing_ambiguity,
|
|
nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, false);
|
|
#endif
|
|
}
|
|
|
|
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;
|
|
}
|