mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 81 Wiki Activity

rugnux: do not make the first validation's maps where the forecast one is kept

Browse Source 
Where the validation in the model's setting, started beside the first on
the forecast (b5980c9ba), is going to be kept, it writes the maps, the map
coefficients and the placed model on the written axes, and the first
validation's p_2fofc/p_fofc/p_anom.ccp4 and p_maps.mtz were overwritten
unread. ModelValidationSchedule::maps_superseded is asked once the
verdict, the setting and the indexing are decided, before any map is
made; the first validation then makes none (the map coefficients and the
MTZ are not built, the three map transforms not run, nothing written).
The answer is the same rule that keeps the forecast validation
afterwards (forecast_kept, now one lambda for both). The kept validation
also brings back the reflections it relabelled, free set included, so
relabel_output does not reindex and redraw them a second time. Where the
forecast is discarded or was never started, everything runs as before,
and modelpar's delete-before-rewrite stays for those cases. If the kept
validation then fails, a warning says no maps were written.

Log lines that disappear (first validation only, kept-forecast runs):
"mean 2mFo-DFc density at atom centres = ...", "anomalous difference map
from ... Bijvoet pairs; strongest density ...", "wrote p_2fofc.ccp4, ...";
in their place: "no maps made here - the validation in the model's
setting writes them". (The "anomalous density ... inverted" warning would
also go; it cannot be raised where a forecast is kept.)

Verified on 8sa8, 8xtg, 9ea5 (forecast kept), 8tyy, 7qis, myob_x10sa (no
forecast): p.mtz, p.hkl, p.cif, the three maps, p_maps.mtz and
p_model.cif/pdb md5-identical to the merge commit, repeated and at -N 8;
the discard path (8sa8 with a rotated model) identical files and log.
Time saved is small with the maps already on the GPU: 8sa8 7.6 -> 7.4 s,
8xtg 5.1 -> 5.0 s, 9ea5 within noise.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
This commit is contained in:
leonarski_f
2026-10-09 14:34:46 +02:00
co-authored by Claude Opus 5.5
parent 278c488c37
commit 9a1f729993
3 changed files with 198 additions and 161 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/structure_refinement/ModelValidation.cpp
+162 -150
View File
@@ -1409,113 +1409,119 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
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);
// --- 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");
});
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();
// --- 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;
}
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());
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;
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) {
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);
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;
}
});
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();
}
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
@@ -1536,55 +1542,57 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
? "the anomalous density at the model's atoms is inverted"
: "the model does not fit these data");
// --- MTZ of map coefficients so the maps can be re-opened / rebuilt in Coot etc. ---
try {
gemmi::Mtz mtz(true);
// The group the REFLECTIONS end up in, which is the model's only where its hand was adopted -
// AdoptModelFrame decides this the same way a few lines below. An enantiomorphic pair indexes
// identically, so the coefficients are the same numbers either way and only the label moves;
// but the two groups have different screw translations, so a reader that expands symmetry out
// of this file works in the wrong one if the label disagrees with the .mtz beside it.
mtz.spacegroup = result.adopted_model_enantiomorph
? sg
// No data group means the caller merged in P1, and P1 is what the
// reflection files beside this one carry - not the model's group.
: (data_space_group != nullptr ? data_space_group
: gemmi::find_spacegroup_by_number(1));
mtz.set_cell_for_all(ucell);
mtz.add_dataset("model_validation");
mtz.datasets.back().wavelength = wavelength_A;
mtz.add_column("FP", 'F', -1, -1, false);
mtz.add_column("FC", 'F', -1, -1, false);
mtz.add_column("PHIC", 'P', -1, -1, false);
mtz.add_column("FWT", 'F', -1, -1, false);
mtz.add_column("PHWT", 'P', -1, -1, false);
mtz.add_column("DELFWT", 'F', -1, -1, false);
mtz.add_column("PHDELWT", 'P', -1, -1, false);
// The figure of merit the coefficients carry, so the weighting can be read off - and undone -
// from the file rather than having to be taken on trust.
mtz.add_column("FOM", 'W', -1, -1, false);
mtz.add_column("FREE", 'I', -1, -1, false);
std::vector<float> data;
for (size_t i = 0; i < terms.size(); i++) {
// [0, 360), the convention the MTZ 'P' columns carried before sigma_A weighting: std::arg
// returns (-180, 180] and a P column is not supposed to.
const double phi_wrapped = terms[i].phi * 180.0 / PI;
const auto phi_deg = static_cast<float>(phi_wrapped < 0 ? phi_wrapped + 360.0 : phi_wrapped);
data.insert(data.end(), {static_cast<float>(terms[i].hkl[0]), static_cast<float>(terms[i].hkl[1]),
static_cast<float>(terms[i].hkl[2]),
static_cast<float>(terms[i].fo), static_cast<float>(terms[i].fc),
phi_deg,
fwt[i], phi_deg,
delfwt[i], phi_deg,
static_cast<float>(sigma_a.weight[i].m),
terms[i].free ? 0.0f : 1.0f});
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());
}
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;
@@ -1596,7 +1604,8 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
result.r_work, result.n_work, result.r_free, result.n_free);
logger.Info("Model validation: bulk solvent k_sol={:.3f} b_sol={:.1f}, k_overall={:.3f}",
result.k_sol, result.b_sol, result.k_overall);
logger.Info("Model validation: mean 2mFo-DFc density at atom centres = {:.2f} sigma", result.mean_atom_density_sigma);
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,
@@ -1615,10 +1624,13 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
"agree - either the model is the wrong enantiomorph for this crystal, or the "
"data were indexed in the wrong hand, and reindexing would hide which",
result.anomalous_deepest_site, result.anomalous_deepest_sigma);
logger.Info("Model validation: wrote {}_2fofc.ccp4, {}_fofc.ccp4{}, {}_maps.mtz",
output_prefix, output_prefix,
result.anomalous_sites.empty() ? "" : ", " + output_prefix + "_anom.ccp4",
output_prefix);
if (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 "