Anomalous difference map: ANOM/PHANOM in _maps.mtz, /rugnux/map.ccp4?type=anom, a magenta mesh in the viewer
Model validation already made the anomalous difference map's coefficients ((F(+) - F(-)) on the model phases less 90 degrees, where the merge kept the Bijvoet pairs) to read the anomalous sites, then threw them away. They now go into the map coefficients as ANOM/PHANOM (phenix.maps' labels), missing on the reflections without a pair, so _maps.mtz opens as an anomalous map in Coot. Nothing in scaling, merging or the validation's decisions changes; a merge without the Bijvoet split writes the twelve columns as before. The broker serves the map as /rugnux/map.ccp4?type=anom (404 without Bijvoet pairs), and rugnux_live_status.anomalous_map says it is there. The Mol* viewer loads it only then, as a magenta mesh at +4 sigma with its own level slider and switch. Tests: ANOM/PHANOM appear with Bijvoet pairs (as many as the reported pairs) and not without; the served anomalous map of a point peaks at it; the auth test covers type=anom. End to end (myoglobin 2M replay): anomalous_map true, the served map equals gemmi's FFT of ANOM/PHANOM in the written _maps.mtz to 3.4e-6 rms. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit c4d6ad74c706181cbb183cfcaa2481653cd7cad6)
This commit is contained in:
1 parent
ee016d3d13
commit
6abc44710d
21 files changed
+224
-52
No files matched your search
@@ -596,6 +596,51 @@ TEST_CASE("ModelValidation_FindsTheDatasDescriptionOfTheLattice", "[ModelValidat
|
||||
std::filesystem::remove(other_setting);
|
||||
}
|
||||
|
||||
TEST_CASE("ModelValidation_MapCoefficientsCarryTheAnomalousMap", "[ModelValidation]") {
|
||||
Logger logger("ModelValidation_MapCoefficientsCarryTheAnomalousMap");
|
||||
const auto data_model = WriteTemp("anom_coef_test.pdb", ClusterPdb().c_str());
|
||||
auto obs = ModelReferenceIntensities(ReadModel(data_model), {}, {}, 2.5, logger);
|
||||
REQUIRE(obs.size() > 500);
|
||||
for (size_t i = 0; i < obs.size(); i++) {
|
||||
obs[i].F = std::sqrt(std::max(0.0f, obs[i].I));
|
||||
obs[i].rfree_flag = (i % 20) == 0;
|
||||
}
|
||||
const UnitCell data_cell{.a = 30, .b = 34, .c = 38, .alpha = 90, .beta = 90, .gamma = 90};
|
||||
const gemmi::SpaceGroup *data_sg = gemmi::find_spacegroup_by_name("P 21 21 21");
|
||||
|
||||
// Without the Bijvoet split: no anomalous columns.
|
||||
const auto plain = ValidateAgainstModel(obs, data_cell, ReadModel(data_model), logger, data_sg, false, 1, 1.0);
|
||||
REQUIRE(plain.ok);
|
||||
REQUIRE(plain.maps);
|
||||
CHECK(plain.maps->coefficients.columns.size() == 12);
|
||||
|
||||
// With it: ANOM / PHANOM beside the others, missing on the reflections without a pair.
|
||||
for (size_t i = 0; i < obs.size(); i++)
|
||||
if (i % 3 != 0) {
|
||||
obs[i].F_plus = obs[i].F * 1.02f;
|
||||
obs[i].F_minus = obs[i].F * 0.98f;
|
||||
}
|
||||
const auto anom = ValidateAgainstModel(obs, data_cell, ReadModel(data_model), logger, data_sg, false, 1, 1.0);
|
||||
REQUIRE(anom.ok);
|
||||
REQUIRE(anom.maps);
|
||||
const auto &columns = anom.maps->coefficients.columns;
|
||||
REQUIRE(columns.size() == 14);
|
||||
CHECK(columns[12].label == "ANOM");
|
||||
CHECK(columns[13].label == "PHANOM");
|
||||
size_t present = 0;
|
||||
for (size_t i = 0; i < columns[12].values.size(); i++)
|
||||
if (std::isfinite(columns[12].values[i])) {
|
||||
++present;
|
||||
CHECK(columns[13].values[i] >= 0.0f);
|
||||
CHECK(columns[13].values[i] < 360.0f);
|
||||
}
|
||||
CHECK(present > 0);
|
||||
CHECK(present < columns[12].values.size());
|
||||
CHECK(static_cast<int>(present) == anom.anomalous_pairs);
|
||||
|
||||
std::filesystem::remove(data_model);
|
||||
}
|
||||
|
||||
// An alternative indexing is settled by relabelling the DATA into the model's indexing, not by moving
|
||||
// the model into the data's: that is what puts every dataset of one crystal form in one convention.
|
||||
// Only where a reference has already fixed the data's indexing is the model moved instead.
|
||||
|
||||
Reference in new issue
Block a user