leonarski_fandClaude Opus 5
f7cb701f43
ValidateAgainstModel was called after WriteReflections, so the two relabelings it makes - the change of hand where the data were merged in the model's enantiomorph, and the alternative indexing it picks by R-free where no reference MTZ had already fixed it - reached the R-factors, the maps and _maps.mtz and nothing else. The .mtz/.cif/.hkl beside them kept whichever indexing the merge happened to land in, so a file whose maps came from a model could not be refined against that model without being reindexed first. Validation now runs before the reflection files are written and returns the operators it applied. AdoptModelFrame puts the merged reflections through them with a proper ASU reduction - ReindexReflections only relabels, which is all its scoring callers need - and the Bijvoet halves follow the Friedel sign of that reduction. The integrated observations are relabelled too, so --export-unmerged describes the same indexing as the merged files rather than a second one. change_of_hand_op is the inversion, so adopting a model's hand exchanges I(+) and I(-). On a tetragonal crystal whose enantiomorphic pair share their whole absence pattern, the search picks one of the two arbitrarily, and the exported anomalous differences were the wrong way round whenever it picked the other member from the model's. The written space group follows the hand, and so does the report that describes those files. The _process.h5 is deliberately left alone: its per-image reflections went to disk as they were integrated, and recording a group that did not match them would mis-merge on a later --mode scale. Nothing is lost - an enantiomorphic pair has the same Laue class and the same absences, so it merges identically either way. Serial stills could not use a model for this at all: each crystal is indexed in its own hand, so by the time there is a merge to fit a model to, the hands have already been averaged together, and reindexing the merged result as a whole cannot separate them again. ModelReferenceIntensities computes |F_model|^2 from the structure - Fcalc plus a flat solvent term at the standard constants, which are not fitted because there are no observations yet - and hands it to the per-image resolver that a reference MTZ already drove. It needs the cell and the group up front, which serial indexing wants anyway. Measured on a serial dataset in a merohedral trigonal group, at a fixed resolution limit so the shells match: CC1/2 better in 8 of 10 shells (79.7 to 82.5 overall), R_meas better in 9 of 10 (111.1 to 104.7), ISa 1.11 to 1.20, and R-free against a structure the merge never saw 0.393 to 0.375. The ambiguity probe now logs the runner-up R-free as well as the winner. On the same data the margin is 0.016 when the hands are mixed, where the two global indexings are tied because both are already inside every reflection, and 0.193 once they are not - the difference between a decision and a coin toss, which the old single number hid. The warning about an unresolved ambiguity no longer fires when a model will resolve it, and names the flags that would; and the summary line no longer offers the adopted group as its own alternative, which it did once --model had moved the run onto the member the search had listed as the alternative. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Vi1gV6Z45aZL5wLwe85Ksn
477 lines
23 KiB
C++
477 lines
23 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "ModelValidation.h"
|
|
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
#include <complex>
|
|
#include <array>
|
|
#include <string>
|
|
#include <vector>
|
|
#include <unordered_map>
|
|
|
|
#include <gemmi/pdb.hpp> // read_pdb
|
|
#include <gemmi/gz.hpp> // MaybeGzipped
|
|
#include <gemmi/it92.hpp> // IT92 x-ray form factors
|
|
#include <gemmi/dencalc.hpp> // DensityCalculator
|
|
#include <gemmi/fourier.hpp> // transform_map_to_f_phi, get_f_phi_on_grid, transform_f_phi_grid_to_map
|
|
#include <gemmi/solmask.hpp> // SolventMasker
|
|
#include <gemmi/scaling.hpp> // Scaling (bulk solvent + anisotropic B)
|
|
#include <gemmi/ccp4.hpp> // Ccp4 map I/O
|
|
#include <gemmi/mtz.hpp> // Mtz (map-coefficient output)
|
|
|
|
#include "../common/Logger.h"
|
|
#include "../image_analysis/scale_merge/ReindexAmbiguity.h" // ReindexReflections
|
|
|
|
namespace {
|
|
|
|
using Table = gemmi::IT92<float>;
|
|
|
|
// Stable key for a Miller index reduced into the ASU (indices are small, well within +/-512).
|
|
long hkl_key(const gemmi::Miller &h) {
|
|
return (h[0] + 512L) * 1048576 + (h[1] + 512L) * 1024 + (h[2] + 512L);
|
|
}
|
|
|
|
// Write an FFT of ASU map coefficients as a CCP4 map; return its RMS (for reporting / sigma units).
|
|
double write_ccp4(gemmi::AsuData<std::complex<float>> &coef, const std::string &path) {
|
|
coef.ensure_sorted();
|
|
std::array<int, 3> size = gemmi::get_size_for_hkl(coef, {{0, 0, 0}}, 3.0);
|
|
gemmi::Grid<float> map =
|
|
gemmi::transform_f_phi_grid_to_map(gemmi::get_f_phi_on_grid<float>(coef, size, true));
|
|
gemmi::Ccp4<float> ccp4;
|
|
ccp4.grid = map;
|
|
ccp4.update_ccp4_header(2);
|
|
ccp4.write_ccp4_map(path);
|
|
return ccp4.hstats.rms;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &merged,
|
|
const UnitCell &cell,
|
|
const std::string &model_path,
|
|
const std::string &output_prefix,
|
|
Logger &logger,
|
|
std::optional<int> data_space_group_number,
|
|
bool probe_indexing_ambiguity) {
|
|
ModelValidationResult result;
|
|
|
|
// --- read the atomic model ---
|
|
gemmi::Structure st;
|
|
try {
|
|
st = gemmi::read_pdb(gemmi::MaybeGzipped(model_path));
|
|
} catch (const std::exception &e) {
|
|
logger.Error("Model validation: cannot read model {}: {}", model_path, e.what());
|
|
return result;
|
|
}
|
|
if (st.models.empty() || !st.cell.is_crystal()) {
|
|
logger.Error("Model validation: model {} has no atoms or no unit cell", model_path);
|
|
return result;
|
|
}
|
|
const gemmi::SpaceGroup *sg = st.find_spacegroup();
|
|
if (!sg) {
|
|
logger.Error("Model validation: model {} has no usable space group", model_path);
|
|
return result;
|
|
}
|
|
result.model_space_group_number = sg->number;
|
|
|
|
// If the data was indexed in the enantiomorph of the model's space group (e.g. data P4(1)2(1)2,
|
|
// model P4(3)2(1)2 - the merged intensities cannot tell them apart), reindex the observed
|
|
// reflections into the model's hand so the two settings agree. This does not change the
|
|
// R-factors (which use |F|), but keeps the observed data consistent with the model.
|
|
std::vector<MergedReflection> reindexed;
|
|
const std::vector<MergedReflection> *obs_ptr = &merged;
|
|
if (data_space_group_number && *data_space_group_number != sg->number) {
|
|
const gemmi::SpaceGroup *dsg = gemmi::find_spacegroup_by_number(*data_space_group_number);
|
|
if (dsg && dsg->is_enantiomorphic() && sg->is_enantiomorphic()) {
|
|
gemmi::GroupOps eops = dsg->operations();
|
|
eops.change_basis_forward(dsg->change_of_hand_op());
|
|
const gemmi::SpaceGroup *enant = gemmi::find_spacegroup_by_ops(eops);
|
|
if (enant && enant->number == sg->number) {
|
|
reindexed = ReindexReflections(merged, dsg->change_of_hand_op());
|
|
obs_ptr = &reindexed;
|
|
result.hand_op = dsg->change_of_hand_op();
|
|
logger.Info("Model validation: data space group {} is the enantiomorph of the model {}; "
|
|
"reindexed the observed reflections into the model's hand",
|
|
dsg->short_name(), sg->hm);
|
|
}
|
|
}
|
|
}
|
|
const std::vector<MergedReflection> &obs = *obs_ptr;
|
|
|
|
// Resolution limit from the data (the merged set is already resolution-trimmed).
|
|
double d_min = 0.0;
|
|
for (const MergedReflection &r : obs)
|
|
if (r.d > 0 && (d_min == 0.0 || r.d < d_min))
|
|
d_min = r.d;
|
|
if (d_min <= 0.0) {
|
|
logger.Error("Model validation: merged reflections carry no resolution");
|
|
return result;
|
|
}
|
|
|
|
// Re-fractionalize the model into the data cell (rigid cell adjustment; no refinement).
|
|
const gemmi::UnitCell data_cell = cell; // UnitCell -> gemmi::UnitCell
|
|
if (data_cell.is_crystal()) {
|
|
gemmi::UnitCell old = st.cell;
|
|
for (gemmi::Model &m : st.models)
|
|
for (gemmi::Chain &ch : m.chains)
|
|
for (gemmi::Residue &r : ch.residues)
|
|
for (gemmi::Atom &a : r.atoms)
|
|
a.pos = data_cell.orthogonalize(old.fractionalize(a.pos));
|
|
st.cell = data_cell;
|
|
}
|
|
st.setup_cell_images();
|
|
const gemmi::UnitCell &ucell = st.cell;
|
|
|
|
logger.Info("Model validation: {} atoms, cell a={:.2f} b={:.2f} c={:.2f}, sg {}, to {:.2f} A",
|
|
gemmi::count_atom_sites(st.models[0]), ucell.a, ucell.b, ucell.c, sg->hm, d_min);
|
|
|
|
// --- Fcalc (atomic) via electron density on a grid + FFT ---
|
|
gemmi::DensityCalculator<Table, float> dc;
|
|
dc.d_min = d_min;
|
|
dc.rate = 1.5;
|
|
dc.set_grid_cell_and_spacegroup(st);
|
|
dc.set_refmac_compatible_blur(st.models[0]);
|
|
dc.put_model_density_on_grid(st.models[0]);
|
|
gemmi::AsuData<std::complex<float>> fcalc =
|
|
gemmi::transform_map_to_f_phi(dc.grid, true).prepare_asu_data(dc.d_min, dc.blur, false, false, false);
|
|
|
|
// --- flat bulk-solvent mask -> Fmask ---
|
|
// Refmac radii give a slightly lower R than the Cctbx set on our test cases, at no cost.
|
|
gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac);
|
|
gemmi::Grid<float> mask_grid;
|
|
mask_grid.unit_cell = dc.grid.unit_cell;
|
|
mask_grid.spacegroup = dc.grid.spacegroup;
|
|
mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
|
|
masker.put_mask_on_grid(mask_grid, st.models[0]);
|
|
gemmi::AsuData<std::complex<float>> fmask =
|
|
gemmi::transform_map_to_f_phi(mask_grid, true).prepare_asu_data(dc.d_min, 0);
|
|
|
|
gemmi::GroupOps gops = sg->operations();
|
|
gemmi::ReciprocalAsu asu(sg);
|
|
|
|
// --- fit the (scaled, solvent-corrected) model to one observed set and score it ---
|
|
// Factored into a lambda so we can probe indexing (merohedral) ambiguities: run the same scale +
|
|
// R computation on each reindexing of the observed reflections and keep the lowest-R-free one.
|
|
struct Fit {
|
|
gemmi::AsuData<std::complex<float>> fmodel, map2fofc, mapfofc;
|
|
std::unordered_map<long, std::pair<double, bool>> obs_by_hkl; // hkl -> (Fobs, is_free)
|
|
double r_work = 1, r_free = 1, k_sol = 0, b_sol = 0, k_overall = 0;
|
|
int n_w = 0, n_f = 0;
|
|
};
|
|
auto fit_model = [&](const std::vector<MergedReflection> &obs_in) -> Fit {
|
|
Fit out;
|
|
out.fmodel = fcalc; // copy the atomic structure factors; scaling mutates them in place
|
|
|
|
// --- observed amplitudes into the model ASU, keyed by hkl (also remember free flag) ---
|
|
// Observed amplitudes are the French-Wilson |F| already computed at the end of the merge
|
|
// (MergedReflection.F), so the model R-free / maps use exactly the same amplitudes as the
|
|
// written reflection file.
|
|
gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
|
|
fobs.unit_cell_ = ucell;
|
|
fobs.spacegroup_ = sg;
|
|
for (const MergedReflection &r : obs_in) {
|
|
if (std::isnan(r.F)) continue;
|
|
gemmi::Miller h{{r.h, r.k, r.l}};
|
|
if (!asu.is_in(h)) h = asu.to_asu(h, gops).first;
|
|
fobs.v.push_back({h, {r.F, 1.0f}});
|
|
out.obs_by_hkl[hkl_key(h)] = {r.F, r.rfree_flag};
|
|
}
|
|
fobs.ensure_asu();
|
|
fobs.ensure_sorted();
|
|
|
|
// --- scale Fmodel(+solvent) to Fobs: k_overall, anisotropic B, k_sol, b_sol ---
|
|
gemmi::Scaling<float> scaling(ucell, sg);
|
|
scaling.use_solvent = true;
|
|
scaling.prepare_points(out.fmodel, fobs, &fmask);
|
|
scaling.fit_isotropic_b_approximately();
|
|
scaling.fit_parameters();
|
|
scaling.scale_data(out.fmodel, &fmask); // out.fmodel now holds the scaled, solvent-corrected Fmodel
|
|
out.k_sol = scaling.k_sol;
|
|
out.b_sol = scaling.b_sol;
|
|
out.k_overall = scaling.k_overall;
|
|
|
|
// The model is scaled to the data with an overall scale, an anisotropic B and a flat bulk
|
|
// solvent only - the standard, few-parameter model that refinement programs use. A dataset-
|
|
// specific free-form per-resolution-shell rescale would lower this dataset's R a little, but
|
|
// it reshapes each map's radial amplitude profile differently, so a batch of maps would no
|
|
// longer be directly comparable. For a fragment-screening / PanDDA campaign, comparable maps
|
|
// across datasets matter more than the last bit of per-dataset R, so it is deliberately omitted.
|
|
|
|
// --- R-work / R-free and map coefficients (2Fo-Fc and Fo-Fc, model phases) ---
|
|
out.map2fofc.unit_cell_ = ucell; out.map2fofc.spacegroup_ = sg;
|
|
out.mapfofc.unit_cell_ = ucell; out.mapfofc.spacegroup_ = sg;
|
|
double num_w = 0, den_w = 0, num_f = 0, den_f = 0;
|
|
for (const auto &hv : out.fmodel.v) {
|
|
auto it = out.obs_by_hkl.find(hkl_key(hv.hkl));
|
|
if (it == out.obs_by_hkl.end()) continue;
|
|
double Fo = it->second.first;
|
|
double Fc = std::abs(hv.value);
|
|
double phi = std::arg(hv.value);
|
|
if (it->second.second) { num_f += std::fabs(Fo - Fc); den_f += Fo; ++out.n_f; }
|
|
else { num_w += std::fabs(Fo - Fc); den_w += Fo; ++out.n_w; }
|
|
std::complex<float> ph = std::polar(1.0f, static_cast<float>(phi));
|
|
out.map2fofc.v.push_back({hv.hkl, static_cast<float>(2 * Fo - Fc) * ph});
|
|
out.mapfofc.v.push_back({hv.hkl, static_cast<float>(Fo - Fc) * ph});
|
|
}
|
|
out.r_work = den_w > 0 ? num_w / den_w : 1;
|
|
out.r_free = den_f > 0 ? num_f / den_f : 1;
|
|
return out;
|
|
};
|
|
|
|
// --- indexing (merohedral) ambiguity ---
|
|
// When a reference MTZ was supplied, the data were already reindexed to agree with the reference
|
|
// intensities (at the merge stage for rotation data, per image in stills scaling), and that
|
|
// choice is authoritative - we keep it. Only with a model and NO reference do we resolve the
|
|
// ambiguity here, as a fallback, by fitting each candidate reindexing and keeping the lowest
|
|
// R-free. A no-op either way for a holohedral crystal (no twin laws). The
|
|
// enantiomorph/screw ambiguity is never probed by R-free: |Fcalc| is the same for both hands, so
|
|
// it cannot distinguish them - that is taken from the model hand above.
|
|
Fit best = fit_model(obs);
|
|
if (probe_indexing_ambiguity) {
|
|
const auto reindex_ops = ReindexAmbiguityOperators(cell, sg->number);
|
|
std::vector<double> candidate_r_free{best.r_free}; // identity first, then the twin laws
|
|
for (const auto &op : reindex_ops) {
|
|
Fit cand = fit_model(ReindexReflections(obs, op));
|
|
candidate_r_free.push_back(cand.r_free);
|
|
if (cand.r_free < best.r_free) { best = std::move(cand); result.indexing_op = op; }
|
|
}
|
|
if (!reindex_ops.empty()) {
|
|
// The runner-up as well as the winner: the margin between them is what says whether the
|
|
// choice was made on evidence, and on weak data the two can come out within noise.
|
|
std::sort(candidate_r_free.begin(), candidate_r_free.end());
|
|
logger.Info("Model validation: no reference - probed {} indexing solution(s) against the model; "
|
|
"{} (R-free {:.4f}, runner-up {:.4f})", candidate_r_free.size(),
|
|
result.indexing_op == gemmi::Op::identity() ? "kept the current indexing"
|
|
: "reindexed to the lower-R-free solution",
|
|
candidate_r_free[0], candidate_r_free[1]);
|
|
}
|
|
}
|
|
|
|
gemmi::AsuData<std::complex<float>> &fmodel = best.fmodel;
|
|
gemmi::AsuData<std::complex<float>> &map2fofc = best.map2fofc;
|
|
gemmi::AsuData<std::complex<float>> &mapfofc = best.mapfofc;
|
|
std::unordered_map<long, std::pair<double, bool>> &obs_by_hkl = best.obs_by_hkl;
|
|
|
|
result.r_work = best.r_work;
|
|
result.r_free = best.r_free;
|
|
result.n_work = best.n_w;
|
|
result.n_free = best.n_f;
|
|
result.k_sol = best.k_sol;
|
|
result.b_sol = best.b_sol;
|
|
result.k_overall = best.k_overall;
|
|
|
|
// --- write the maps and score the 2Fo-Fc map at atom centres (a real map peaks there) ---
|
|
const std::string p2 = output_prefix + "_2fofc.ccp4";
|
|
const std::string pd = output_prefix + "_fofc.ccp4";
|
|
double rms2 = write_ccp4(map2fofc, p2);
|
|
write_ccp4(mapfofc, pd);
|
|
|
|
// Recompute the 2Fo-Fc real-space grid once more to sample it at atom positions.
|
|
{
|
|
map2fofc.ensure_sorted();
|
|
std::array<int, 3> size = gemmi::get_size_for_hkl(map2fofc, {{0, 0, 0}}, 3.0);
|
|
gemmi::Grid<float> grid =
|
|
gemmi::transform_f_phi_grid_to_map(gemmi::get_f_phi_on_grid<float>(map2fofc, size, true));
|
|
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 += grid.interpolate_value(a.pos); ++n; }
|
|
result.mean_atom_density_sigma = (n > 0 && rms2 > 0) ? (s / n) / rms2 : 0;
|
|
}
|
|
|
|
// --- MTZ of map coefficients so the maps can be re-opened / rebuilt in Coot etc. ---
|
|
try {
|
|
gemmi::Mtz mtz(true);
|
|
mtz.spacegroup = sg;
|
|
mtz.set_cell_for_all(ucell);
|
|
mtz.add_dataset("model_validation");
|
|
mtz.add_column("FP", 'F', -1, -1, false);
|
|
mtz.add_column("FC", 'F', -1, -1, false);
|
|
mtz.add_column("PHIC", 'P', -1, -1, false);
|
|
mtz.add_column("FWT", 'F', -1, -1, false);
|
|
mtz.add_column("PHWT", 'P', -1, -1, false);
|
|
mtz.add_column("DELFWT", 'F', -1, -1, false);
|
|
mtz.add_column("PHDELWT", 'P', -1, -1, false);
|
|
mtz.add_column("FREE", 'I', -1, -1, false);
|
|
std::vector<float> data;
|
|
int nref = 0;
|
|
for (const auto &hv : fmodel.v) {
|
|
auto it = obs_by_hkl.find(hkl_key(hv.hkl));
|
|
if (it == obs_by_hkl.end()) continue;
|
|
double Fo = it->second.first;
|
|
double Fc = std::abs(hv.value);
|
|
double phi_deg = gemmi::phase_in_angles(hv.value);
|
|
data.insert(data.end(), {static_cast<float>(hv.hkl[0]), static_cast<float>(hv.hkl[1]),
|
|
static_cast<float>(hv.hkl[2]),
|
|
static_cast<float>(Fo), static_cast<float>(Fc),
|
|
static_cast<float>(phi_deg),
|
|
static_cast<float>(2 * Fo - Fc), static_cast<float>(phi_deg),
|
|
static_cast<float>(Fo - Fc), static_cast<float>(phi_deg),
|
|
it->second.second ? 0.0f : 1.0f});
|
|
++nref;
|
|
}
|
|
mtz.nreflections = nref;
|
|
mtz.data = std::move(data);
|
|
mtz.write_to_file(output_prefix + "_maps.mtz");
|
|
} catch (const std::exception &e) {
|
|
logger.Warning("Model validation: could not write map MTZ: {}", e.what());
|
|
}
|
|
|
|
result.ok = true;
|
|
result.maps_prefix = output_prefix;
|
|
logger.Info("Model validation: R-work={:.4f} ({} refl) R-free={:.4f} ({} refl) "
|
|
"[overall + anisotropic B + bulk solvent]",
|
|
result.r_work, result.n_work, result.r_free, result.n_free);
|
|
logger.Info("Model validation: bulk solvent k_sol={:.3f} b_sol={:.1f}, k_overall={:.3f}",
|
|
result.k_sol, result.b_sol, result.k_overall);
|
|
logger.Info("Model validation: mean 2Fo-Fc density at atom centres = {:.2f} sigma", result.mean_atom_density_sigma);
|
|
logger.Info("Model validation: wrote {}_2fofc.ccp4, {}_fofc.ccp4, {}_maps.mtz",
|
|
output_prefix, output_prefix, output_prefix);
|
|
return result;
|
|
}
|
|
|
|
namespace {
|
|
|
|
// The reindexing operator as it reads on Miller indices ("k,h,-l" rather than "y,x,-z").
|
|
std::string hkl_triplet(const gemmi::Op &op) {
|
|
std::string t = op.triplet();
|
|
std::replace(t.begin(), t.end(), 'x', 'h');
|
|
std::replace(t.begin(), t.end(), 'y', 'k');
|
|
std::replace(t.begin(), t.end(), 'z', 'l');
|
|
return t;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
int AdoptModelFrame(const ModelValidationResult &validation,
|
|
std::vector<MergedReflection> &merged,
|
|
int data_space_group_number,
|
|
bool merge_friedel,
|
|
Logger &logger) {
|
|
int space_group_number = data_space_group_number;
|
|
if (!validation.ok)
|
|
return space_group_number;
|
|
|
|
// The change of hand also changes the space group the reflections are in: P4(1)2(1)2 data written
|
|
// in the model's hand are P4(3)2(1)2 data. The merged intensities are the same measurements either
|
|
// way - the enantiomorph is exactly what they cannot distinguish - so the model's hand is the only
|
|
// evidence there is, and adopting it is what makes the file refinable against that model.
|
|
if (!(validation.hand_op == gemmi::Op::identity()) && validation.model_space_group_number > 0) {
|
|
merged = ReindexMergedIntoAsu(merged, validation.hand_op,
|
|
validation.model_space_group_number, merge_friedel);
|
|
space_group_number = validation.model_space_group_number;
|
|
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number);
|
|
logger.Info("Model validation: the written reflections take the model's enantiomorph, {} ({}), "
|
|
"reindexed by {}", sg ? sg->short_name() : "?", space_group_number,
|
|
hkl_triplet(validation.hand_op));
|
|
}
|
|
|
|
// The alternative indexing, by contrast, is metric- and group-preserving: only the labels move.
|
|
if (!(validation.indexing_op == gemmi::Op::identity())) {
|
|
merged = ReindexMergedIntoAsu(merged, validation.indexing_op, space_group_number, merge_friedel);
|
|
logger.Info("Model validation: the written reflections take the model's indexing, reindexed by {}",
|
|
hkl_triplet(validation.indexing_op));
|
|
}
|
|
return space_group_number;
|
|
}
|
|
|
|
std::vector<MergedReflection> ModelReferenceIntensities(const std::string &model_path,
|
|
const std::optional<UnitCell> &cell,
|
|
std::optional<int> space_group_number,
|
|
double d_min,
|
|
Logger &logger) {
|
|
std::vector<MergedReflection> out;
|
|
if (!(d_min > 0.0)) {
|
|
logger.Warning("Model reference: no resolution limit to compute the model intensities to");
|
|
return out;
|
|
}
|
|
|
|
gemmi::Structure st;
|
|
try {
|
|
st = gemmi::read_pdb(gemmi::MaybeGzipped(model_path));
|
|
} catch (const std::exception &e) {
|
|
logger.Error("Model reference: cannot read model {}: {}", model_path, e.what());
|
|
return out;
|
|
}
|
|
if (st.models.empty() || !st.cell.is_crystal()) {
|
|
logger.Error("Model reference: model {} has no atoms or no unit cell", model_path);
|
|
return out;
|
|
}
|
|
|
|
// Put the model in the cell and group the run works in, where it knows them, so the reference is
|
|
// indexed the way the data are. The correlation that consumes this matches on hkl, so a small cell
|
|
// difference costs nothing; the space group is what has to agree.
|
|
if (cell.has_value()) {
|
|
const gemmi::UnitCell target = *cell;
|
|
if (target.is_crystal()) {
|
|
const gemmi::UnitCell old = st.cell;
|
|
for (gemmi::Model &m : st.models)
|
|
for (gemmi::Chain &ch : m.chains)
|
|
for (gemmi::Residue &r : ch.residues)
|
|
for (gemmi::Atom &a : r.atoms)
|
|
a.pos = target.orthogonalize(old.fractionalize(a.pos));
|
|
st.cell = target;
|
|
}
|
|
}
|
|
if (space_group_number.has_value())
|
|
if (const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(*space_group_number))
|
|
st.spacegroup_hm = sg->xhm();
|
|
const gemmi::SpaceGroup *sg = st.find_spacegroup();
|
|
if (!sg) {
|
|
logger.Error("Model reference: model {} has no usable space group", model_path);
|
|
return out;
|
|
}
|
|
st.setup_cell_images();
|
|
|
|
gemmi::DensityCalculator<Table, float> dc;
|
|
dc.d_min = d_min;
|
|
dc.rate = 1.5;
|
|
dc.set_grid_cell_and_spacegroup(st);
|
|
dc.set_refmac_compatible_blur(st.models[0]);
|
|
dc.put_model_density_on_grid(st.models[0]);
|
|
gemmi::AsuData<std::complex<float>> fcalc =
|
|
gemmi::transform_map_to_f_phi(dc.grid, true).prepare_asu_data(dc.d_min, dc.blur, false, false, false);
|
|
|
|
// Flat bulk solvent at the standard constants. Nothing here is fitted - there are no observations
|
|
// yet - but without it the few lowest-resolution reflections are the largest and the most wrong,
|
|
// and a correlation on raw intensities would be led by them.
|
|
constexpr double K_SOL = 0.35;
|
|
constexpr double B_SOL = 46.0;
|
|
gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac);
|
|
gemmi::Grid<float> mask_grid;
|
|
mask_grid.unit_cell = dc.grid.unit_cell;
|
|
mask_grid.spacegroup = dc.grid.spacegroup;
|
|
mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
|
|
masker.put_mask_on_grid(mask_grid, st.models[0]);
|
|
gemmi::AsuData<std::complex<float>> fmask =
|
|
gemmi::transform_map_to_f_phi(mask_grid, true).prepare_asu_data(dc.d_min, 0);
|
|
|
|
std::unordered_map<long, std::complex<float>> mask_by_hkl;
|
|
mask_by_hkl.reserve(fmask.v.size());
|
|
for (const auto &hv : fmask.v)
|
|
mask_by_hkl[hkl_key(hv.hkl)] = hv.value;
|
|
|
|
const gemmi::UnitCell &ucell = st.cell;
|
|
out.reserve(fcalc.v.size());
|
|
for (const auto &hv : fcalc.v) {
|
|
const double d = ucell.calculate_d(hv.hkl);
|
|
if (!(d > 0.0))
|
|
continue;
|
|
std::complex<float> f = hv.value;
|
|
const auto it = mask_by_hkl.find(hkl_key(hv.hkl));
|
|
if (it != mask_by_hkl.end())
|
|
f += static_cast<float>(K_SOL * std::exp(-B_SOL / (4.0 * d * d))) * it->second;
|
|
const double F = std::abs(f);
|
|
out.push_back(MergedReflection{.h = hv.hkl[0], .k = hv.hkl[1], .l = hv.hkl[2],
|
|
.I = static_cast<float>(F * F),
|
|
.d = static_cast<float>(d)});
|
|
}
|
|
|
|
logger.Info("Model reference: {} intensities computed from {} to {:.2f} A, space group {}",
|
|
out.size(), model_path, d_min, sg->short_name());
|
|
return out;
|
|
}
|