From 54cf04d1594f442b41206d7ae847610538655dee Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 17:53:00 +0100 Subject: [PATCH] Map coefficients: FWT and DELFWT written as non-negative amplitudes, the phase turned by 180 degrees where the coefficient is negative The refinement's and model validation's _maps.mtz carried 2mFo-DFc and mFo-DFc as signed numbers on the model phase, so DELFWT (and FWT where DFc > 2mFo) went negative. Coot reads that, but an MTZ amplitude column is non-negative by convention and other readers need not cope. The structure factors - and so every map made from them, in rugnux and in the broker - are the same numbers. The phase columns are also wrapped to [0, 360) after the rounding to float, which could give exactly 360. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- .../structure_refinement/ModelValidation.cpp | 19 ++++++++++++++----- .../structure_refinement/Refinement.cpp | 18 ++++++++++++------ tests/ModelValidationTest.cpp | 10 ++++++++++ tests/RefinementTest.cpp | 11 +++++++++++ 4 files changed, 47 insertions(+), 11 deletions(-) diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index c2b2ca749..b53d4e210 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -1835,15 +1835,24 @@ ModelValidationResult Validate(const std::vector &merged, anom_by_hkl[hkl_key(hv.hkl)] = hv.value; 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(phi_wrapped < 0 ? phi_wrapped + 360.0 : phi_wrapped); + // returns (-180, 180] and a P column is not supposed to - nor 360 itself, which the rounding to float + // can give. + const double deg = terms[i].phi * 180.0 / PI; + const auto wrap = [](double d) { + const auto f = static_cast(d - 360.0 * std::floor(d / 360.0)); + return f < 360.0f ? f : 0.0f; + }; + const float phi_deg = wrap(deg); + // A negative coefficient is written as its magnitude on the phase turned by 180 degrees: an MTZ + // amplitude column is non-negative by convention, and not every reader accepts otherwise. + const float phwt = wrap(fwt[i] < 0 ? deg + 180.0 : deg); + const float phdelwt = wrap(delfwt[i] < 0 ? deg + 180.0 : deg); const float row[] = {static_cast(terms[i].hkl[0]), static_cast(terms[i].hkl[1]), static_cast(terms[i].hkl[2]), static_cast(terms[i].fo), static_cast(terms[i].fc), phi_deg, - fwt[i], phi_deg, - delfwt[i], phi_deg, + std::abs(fwt[i]), phwt, + std::abs(delfwt[i]), phdelwt, static_cast(sigma_a.weight[i].m), terms[i].free ? 0.0f : 1.0f}; for (size_t c = 0; c < mtz.columns.size(); c++) diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index 11ea25145..c273f43ab 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -1150,6 +1150,12 @@ double BalancedADPWeight(Context &ctx, const RestraintSet &b_restraints, const s return nx > 0 ? std::sqrt(ng / nx) : 1.0; } +// A phase in degrees as an MTZ 'P' column holds it: in [0, 360), also after the rounding to float. +float MtzPhase(double deg) { + const auto f = static_cast(deg - 360.0 * std::floor(deg / 360.0)); + return f < 360.0f ? f : 0.0f; +} + // The model's map coefficients and its difference map, weighted by the refinement's own likelihood: m its figure // of merit and D its D (2mFo - DFc, mFo - DFc). With a twin law, Fo is detwinned by the model in proportion to its // own share of the twinned intensity, Fo |F_model(h)| / |F_twin(h)|, and Fc is the model's own |F_model(h)|. @@ -1200,14 +1206,14 @@ void MakeMaps(Context &ctx, const std::vector &merged, double const double fwt = ctx.refl.centric[m] ? mw * fo : 2 * mw * fo - d * fc; const double delfwt = mw * fo - d * fc; const double phi = std::arg(fmodel[m]); - double phi_deg = phi * 180.0 / PI; - if (phi_deg < 0) - phi_deg += 360.0; + const double phi_deg = phi * 180.0 / PI; + // A negative coefficient is written as its magnitude on the phase turned by 180 degrees: an MTZ amplitude + // column is non-negative by convention, and not every reader accepts otherwise. const float row[] = {static_cast(ctx.refl.hkl[m][0]), static_cast(ctx.refl.hkl[m][1]), static_cast(ctx.refl.hkl[m][2]), static_cast(fo), static_cast(fc), - static_cast(phi_deg), static_cast(fwt), static_cast(phi_deg), - static_cast(delfwt), static_cast(phi_deg), static_cast(mw), - ctx.refl.free[m] ? 0.0f : 1.0f}; + MtzPhase(phi_deg), static_cast(std::abs(fwt)), MtzPhase(fwt < 0 ? phi_deg + 180.0 : phi_deg), + static_cast(std::abs(delfwt)), MtzPhase(delfwt < 0 ? phi_deg + 180.0 : phi_deg), + static_cast(mw), ctx.refl.free[m] ? 0.0f : 1.0f}; for (size_t c = 0; c < mtz.columns.size(); c++) mtz.columns[c].values.push_back(row[c]); diff.v.push_back({ctx.refl.hkl[m], std::polar(static_cast(delfwt), static_cast(phi))}); diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index 7b622a91d..aff6193e6 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -613,6 +613,16 @@ TEST_CASE("ModelValidation_MapCoefficientsCarryTheAnomalousMap", "[ModelValidati REQUIRE(plain.ok); REQUIRE(plain.maps); CHECK(plain.maps->coefficients.columns.size() == 12); + // FWT and DELFWT as amplitude and phase: non-negative, their phases in [0, 360) + for (size_t c : {6, 8}) { + const auto &f = plain.maps->coefficients.columns[c], &phi = plain.maps->coefficients.columns[c + 1]; + CHECK(f.label == (c == 6 ? "FWT" : "DELFWT")); + for (size_t i = 0; i < f.values.size(); i++) { + CHECK(f.values[i] >= 0.0f); + CHECK(phi.values[i] >= 0.0f); + CHECK(phi.values[i] < 360.0f); + } + } // With it: ANOM / PHANOM beside the others, missing on the reflections without a pair. for (size_t i = 0; i < obs.size(); i++) diff --git a/tests/RefinementTest.cpp b/tests/RefinementTest.cpp index 215da9d20..1dd5b3990 100644 --- a/tests/RefinementTest.cpp +++ b/tests/RefinementTest.cpp @@ -834,6 +834,17 @@ TEST_CASE("Refinement_restrained_recovers_model", "[Refinement]") { rms += refined[i]->pos.dist_sq(ref[i]->pos); rms = std::sqrt(rms / ref.size()); CHECK(rms < 0.15); + // the map coefficients as amplitude and phase: FWT and DELFWT non-negative, their phases in [0, 360) + REQUIRE(r.maps); + const auto &columns = r.maps->coefficients.columns; + REQUIRE(columns[6].label == "FWT"); + REQUIRE(columns[8].label == "DELFWT"); + for (size_t c : {6, 8}) + for (size_t i = 0; i < columns[c].values.size(); i++) { + CHECK(columns[c].values[i] >= 0.0f); + CHECK(columns[c + 1].values[i] >= 0.0f); + CHECK(columns[c + 1].values[i] < 360.0f); + } // one trial at a time instead of three side by side (sharing a card's tables, mask and scale fit): the same // model, bit for bit