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

rugnux --af-model: tNCS half-vector moves checked, the merge's own group not ruled out beside tNCS, P1's first copy judged by its orientation (WIP, validation pending)

Browse Source 
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3
This commit is contained in:
leonarski_f
2026-10-12 16:49:02 +02:00
co-authored by Claude Opus 5.5
parent 6c8ead466a
commit c6912e2a06
2 changed files with 168 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/ACKNOWLEDGEMENT.md
+6 -1
View File
@@ -481,7 +481,12 @@ boundaries of molecular replacement with maximum likelihood" (2001), Acta Cryst.
Z-score over the same orientation at other positions, as the likelihood-enhanced translation function's
is: A. J. McCoy, R. W. Grosse-Kunstleve, L. C. Storoni and R. J. Read, "Likelihood-enhanced fast
translation functions" (2005), Acta Cryst. D61, 458-464
[doi:10.1107/S0907444905001617](https://doi.org/10.1107/S0907444905001617). The logarithm of the Bessel
[doi:10.1107/S0907444905001617](https://doi.org/10.1107/S0907444905001617). A model whose expected
log-likelihood gain - the sum of sigma_A^4 / 2 over the reflections - is too small for a placement to stand
out is not searched for: A. J. McCoy, R. D. Oeffner, A. G. Wrobel, J. R. M. Ojala, K. Tryggvason,
B. Lohkamp and R. J. Read, "Ab initio solution of macromolecular crystal structures without direct
methods" (2017), Proc. Natl. Acad. Sci. USA 114, 3637-3641
[doi:10.1073/pnas.1701640114](https://doi.org/10.1073/pnas.1701640114). The logarithm of the Bessel
function I0 in it is evaluated with the polynomial approximations of M. Abramowitz and I. A. Stegun,
"Handbook of Mathematical Functions" (1964), National Bureau of Standards, 9.8.1-9.8.2. A predicted
model is cut into rigid domains by its predicted aligned error, as Oeffner et al. (above) do before
image_analysis/structure_refinement/MolecularReplacement.cpp
+162 -10
View File
@@ -102,6 +102,14 @@ constexpr double NULL_SIGMA_MIN = 4.51;
// A model whose expected log-likelihood gain (ExpectedLogLikelihoodGain) is below this is not searched for. The
// lowest of a model the search did place, over 193 survey sets: 26; every PAE domain below 24 placed nowhere.
constexpr double ELLG_MIN = 10.0;
// tNCS by half a lattice vector v: copies v apart make the reflections with h.v half-integral weak, and a
// placement moved by (v + o) / 2, o an origin shift of the group, fits the strong ones as well as the right one
// does; only the weak ones tell them apart. A placement is compared with every such move of itself on those
// reflections (R_tncs_weak), and one of them that fits them better by TNCS_SHIFT_MARGIN is taken instead; where the
// two best are closer than that, the placement is not decided on. A vector is taken as half-integral within
// TNCS_HALF_TOLERANCE.
constexpr double TNCS_SHIFT_MARGIN = 0.01;
constexpr double TNCS_HALF_TOLERANCE = 0.05;
double seconds_since(std::chrono::steady_clock::time_point t0) {
return std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count();
@@ -161,6 +169,9 @@ struct PlacementScore {
// Against the data as twinned by a law (score_placement's `law`): at the fraction that gives the lowest
// R-work.
double r_work_twin = NAN, r_free_twin = NAN, twin_fraction = NAN;
// Over the reflections a tNCS vector v weakens - h.v half-integral - work and free together (score_placement's
// `tncs`).
double r_tncs_weak = NAN;
};
gemmi::AsuData<gemmi::ValueSigma<float>> working_amplitudes(const std::vector<MergedReflection> &merged,
@@ -205,37 +216,54 @@ void model_factors(const gemmi::Structure &st, double d_min, gemmi::AsuData<std:
// R-work and R-free of the scaled Fcalc over the reflections to d_min.
PlacementScore r_factors(const gemmi::AsuData<std::complex<float>> &fmodel, const gemmi::Scaling<float> &scaling,
const gemmi::SpaceGroup &sg, const std::vector<MergedReflection> &merged, double d_min) {
const gemmi::SpaceGroup &sg, const std::vector<MergedReflection> &merged, double d_min,
const gemmi::Fractional *tncs = nullptr) {
const gemmi::GroupOps gops = sg.operations();
const gemmi::ReciprocalAsu asu(&sg);
std::unordered_map<long, std::pair<double, bool>> obs;
struct Obs {
double f;
bool free;
bool weak; // h.v half-integral
};
std::unordered_map<long, Obs> obs;
for (const MergedReflection &r : merged) {
if (std::isnan(r.F) || r.d < d_min)
continue;
gemmi::Miller h{{r.h, r.k, r.l}};
if (!asu.is_in(h))
h = asu.to_asu(h, gops).first;
obs[hkl_key(h)] = {r.F, r.rfree_flag};
bool weak = false;
if (tncs != nullptr) {
const double hv = h[0] * tncs->x + h[1] * tncs->y + h[2] * tncs->z;
weak = std::fabs(hv - std::round(hv)) > 0.25;
}
obs[hkl_key(h)] = {r.F, r.rfree_flag, weak};
}
double num_w = 0, den_w = 0, num_f = 0, den_f = 0;
double num_w = 0, den_w = 0, num_f = 0, den_f = 0, num_t = 0, den_t = 0;
for (const auto &hv : fmodel.v) {
const auto it = obs.find(hkl_key(hv.hkl));
if (it == obs.end())
continue;
const double fo = it->second.first, fc = std::abs(hv.value);
if (it->second.second) {
const double fo = it->second.f, fc = std::abs(hv.value);
if (it->second.free) {
num_f += std::fabs(fo - fc);
den_f += fo;
} else {
num_w += std::fabs(fo - fc);
den_w += fo;
}
if (it->second.weak) {
num_t += std::fabs(fo - fc);
den_t += fo;
}
}
PlacementScore out;
out.k_sol = scaling.k_sol;
out.b_sol = scaling.b_sol;
out.r_work = den_w > 0 ? num_w / den_w : 1;
out.r_free = den_f > 0 ? num_f / den_f : 1;
if (tncs != nullptr)
out.r_tncs_weak = den_t > 0 ? num_t / den_t : 1;
return out;
}
@@ -423,7 +451,7 @@ void twin_r_factors(const gemmi::AsuData<std::complex<float>> &fmodel, const gem
PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGroup &sg,
const std::vector<MergedReflection> &merged, double d_min, size_t nthreads,
ScoreEngines *engines = nullptr, const PlacementScore *solvent = nullptr,
const gemmi::Op *law = nullptr) {
const gemmi::Op *law = nullptr, const gemmi::Fractional *tncs = nullptr) {
gemmi::AsuData<std::complex<float>> fmask, fmodel;
#ifdef JFJOCH_USE_CUDA
ScoreEngines::Engine *gpu = engines != nullptr ? engines->get(st, sg, d_min) : nullptr;
@@ -454,7 +482,7 @@ PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGro
FitModelScale(scaling, {.b_hi = SolventBMax(st.models[0])}, nthreads);
}
scaling.scale_data(fmodel, &fmask);
PlacementScore out = r_factors(fmodel, scaling, sg, merged, d_min);
PlacementScore out = r_factors(fmodel, scaling, sg, merged, d_min, tncs);
if (law != nullptr)
twin_r_factors(fmodel, sg, merged, d_min, *law, out);
return out;
@@ -651,6 +679,33 @@ bool enantiomorphs(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b) {
return other != nullptr && other->number == b.number;
}
// The tNCS vector where it is half a lattice vector - every component within TNCS_HALF_TOLERANCE of 0 or 1/2 -
// on those halves; empty otherwise.
std::optional<gemmi::Fractional> half_tncs_vector(const TranslationalNcs &tncs) {
if (!tncs.present)
return std::nullopt;
double v[3] = {tncs.vector.x, tncs.vector.y, tncs.vector.z};
for (double &x : v) {
const double h = std::round(2 * x) / 2;
if (std::fabs(x - h) > TNCS_HALF_TOLERANCE)
return std::nullopt;
x = h - std::floor(h);
}
return gemmi::Fractional(v[0], v[1], v[2]);
}
// Whether moving the origin by `s` (a multiple of 1/24 in each component) leaves `sg`'s operators as they are:
// a placement moved by s is then the same crystal.
bool is_origin_shift(const gemmi::SpaceGroup &sg, const gemmi::Fractional &s) {
gemmi::Op shift = gemmi::Op::identity();
const double x[3] = {s.x, s.y, s.z};
for (int a = 0; a < 3; a++)
shift.tran[a] = static_cast<int>(std::lround(x[a] * gemmi::Op::DEN)) % gemmi::Op::DEN;
gemmi::GroupOps ops = sg.operations();
ops.change_basis_forward(shift);
return ops.is_same_as(sg.operations());
}
} // namespace
// --- the model ---
@@ -2531,8 +2586,11 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec
cand.verdict = "CLASH";
// The likelihood's Z-score stands in for the translation function's only where the copies placed
// are as many as the cell most probably holds: beside content the search has not found, a copy's
// likelihood over random positions of its own is no measure of the placement.
else if (!(cand.tf_z >= TF_Z_PLACED || (cand.llg_z >= TF_Z_PLACED && cand.copies >= cand.copies_likely)))
// likelihood over random positions of its own is no measure of the placement. In P1 the first copy sits
// at the origin, which is arbitrary there and which no function of its position can single out: one copy
// is judged by its orientation alone, against the null.
else if (!(cand.tf_z >= TF_Z_PLACED || (cand.llg_z >= TF_Z_PLACED && cand.copies >= cand.copies_likely)
|| (gops.order() == 1 && cand.copies == 1)))
cand.verdict = "NOT_PLACED";
else
cand.verdict = "PLACED";
@@ -2630,6 +2688,79 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec
// by side, each taken where the choice comes to it.
std::vector<double> nulls(result.candidates.size(), NAN);
bool first = true;
// With tNCS by half a lattice vector, candidate i against its moves by half that vector (TNCS_SHIFT_MARGIN):
// the better taken into it, and whether the two best are too close to tell apart.
// The tNCS vector is the answer's merge's, in its cell.
const std::optional<gemmi::Fractional> half = merge == 0 ? half_tncs_vector(result.tncs) : std::nullopt;
std::vector<char> shift_checked(result.candidates.size(), 0);
const auto check_shifts = [&](int i) -> bool {
MRCandidate &c = result.candidates[i];
const gemmi::Fractional v = *half;
const std::vector<MergedReflection> &merged = *merges[merge].merged;
const double d_min = result.d_min_score;
struct Option {
gemmi::Fractional shift;
std::shared_ptr<gemmi::Structure> placed;
PlacementScore score;
double clash = 0;
};
std::vector<Option> options{{gemmi::Fractional(0, 0, 0), c.placed,
score_placement(*c.placed, c.space_group, merged, d_min, nthreads, &scorers,
nullptr, nullptr, &v), c.clash_fraction}};
RigidBodyEngines engines;
for (int o = 0; o < 8; o++) {
const gemmi::Fractional shift((v.x + 0.5 * (o & 1)) / 2, (v.y + 0.5 * ((o >> 1) & 1)) / 2,
(v.z + 0.5 * ((o >> 2) & 1)) / 2);
// The same crystal as one already tried, or as no move at all.
if (std::any_of(options.begin(), options.end(), [&](const Option &p) {
return is_origin_shift(c.space_group, gemmi::Fractional(shift.x - p.shift.x, shift.y - p.shift.y,
shift.z - p.shift.z));
}))
continue;
auto moved = std::make_shared<gemmi::Structure>(*c.placed);
const gemmi::Position by = moved->cell.orthogonalize(shift);
std::vector<gemmi::Position> pos = ModelPositions(moved->models[0]);
for (gemmi::Position &p : pos)
p = gemmi::Position(p + by);
SetModelPositions(moved->models[0], pos);
Option opt{shift, moved, {}, ClashFraction(*moved, c.copies)};
if (opt.clash <= 2 * CLASH_MAX) {
if (engines.get() == nullptr)
engines = make_engines(*moved, c.space_group, merged, d_min, 1, logger);
refine_rigid_body(*moved, c.space_group, merged, d_min, engines.get(), nthreads);
opt.clash = ClashFraction(*moved, c.copies);
opt.score = score_placement(*moved, c.space_group, merged, d_min, nthreads, &scorers, nullptr,
nullptr, &v);
}
logger.Info("Model prior: {} moved by {:.2f} {:.2f} {:.2f}: clash {:.1f}%, R-free {:.3f}, R on the "
"reflections tNCS weakens {:.3f} against {:.3f}", c.space_group.xhm(), shift.x, shift.y,
shift.z, 100 * opt.clash, opt.score.r_free, opt.score.r_tncs_weak,
options[0].score.r_tncs_weak);
options.push_back(std::move(opt));
}
std::vector<const Option *> packing;
for (const Option &p : options)
if (p.clash <= CLASH_MAX && std::isfinite(p.score.r_tncs_weak))
packing.push_back(&p);
std::stable_sort(packing.begin(), packing.end(), [](const Option *a, const Option *b) {
return a->score.r_tncs_weak < b->score.r_tncs_weak;
});
if (packing.empty() || packing[0] == &options[0])
return packing.size() > 1 && packing[1]->score.r_tncs_weak - packing[0]->score.r_tncs_weak < TNCS_SHIFT_MARGIN;
const Option &b = *packing[0];
if (options[0].score.r_tncs_weak - b.score.r_tncs_weak < TNCS_SHIFT_MARGIN
|| (packing.size() > 1 && packing[1]->score.r_tncs_weak - b.score.r_tncs_weak < TNCS_SHIFT_MARGIN))
return true;
logger.Info("Model prior: {} taken moved by {:.2f} {:.2f} {:.2f}", c.space_group.xhm(), b.shift.x,
b.shift.y, b.shift.z);
c.placed = b.placed;
c.clash_fraction = b.clash;
c.r_work = b.score.r_work;
c.r_free = b.score.r_free;
c.k_sol = b.score.k_sol;
c.b_sol = b.score.b_sol;
return false;
};
int best = -1;
while (true) {
best = -1;
@@ -2642,6 +2773,15 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec
if (best < 0)
return -1;
MRCandidate &c = result.candidates[best];
if (half && !shift_checked[best]) {
shift_checked[best] = 1;
if (check_shifts(best)) {
result.decision = fmt::format("UNDECIDED: {}'s placement and its move by half the tNCS vector fit "
"alike", c.space_group.xhm());
return -2;
}
continue; // its R-free may have changed
}
if (std::isnan(nulls[best])) {
const auto t_null = std::chrono::steady_clock::now();
std::vector<int> todo{best};
@@ -2665,6 +2805,18 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec
}
for (size_t i = 0; i < result.candidates.size(); i++) {
const MRCandidate &c = result.candidates[i];
// With tNCS, the group the merge was made in, searched to fewer copies than the chosen one - whether it
// packs or not: its first copy alone can sit where half the tNCS vector puts it - is not ruled out by
// any placement of the chosen one, which may sit there itself. (Its enantiomorph is the hand check's.)
if (result.tncs.present && static_cast<int>(i) != best && c.merge == merge
&& c.space_group.xhm() == merges[merge].space_group.xhm()
&& result.candidates[best].space_group.xhm() != c.space_group.xhm()
&& !enantiomorphs(c.space_group, result.candidates[best].space_group)
&& c.copies < result.candidates[best].copies) {
result.decision = fmt::format("UNDECIDED: {} placed fewer copies than {} beside tNCS",
c.space_group.xhm(), result.candidates[best].space_group.xhm());
return -2;
}
if (static_cast<int>(i) == best || c.merge != merge || c.clash_fraction > CLASH_MAX)
continue;
// With tNCS one copy alone fits no better than at random, so a group whose search stopped at