diff --git a/image_analysis/indexing/FFTIndexer.cpp b/image_analysis/indexing/FFTIndexer.cpp index 7d2aa035..86cabc4e 100644 --- a/image_analysis/indexing/FFTIndexer.cpp +++ b/image_analysis/indexing/FFTIndexer.cpp @@ -170,9 +170,10 @@ std::vector FFTIndexer::FilterFFTResults(size_t max_vectors) const { std::vector fft_result_filtered; int count = 0; - for (auto it = fft_result_map.rbegin(); - it != fft_result_map.rend() && count < max_vectors; - ++it, ++count) { + // `it` outlives the loop: the extra scan at the end of the function carries on from the + // first direction the budget did not reach. + auto it = fft_result_map.rbegin(); + for (; it != fft_result_map.rend() && count < max_vectors; ++it, ++count) { fft_result_filtered.emplace_back(it->second); } @@ -237,6 +238,36 @@ std::vector FFTIndexer::FilterFFTResults(size_t max_vectors) const { return A.Length() < B.Length(); }); + // max_vectors counts RAW search directions, but one lattice row is sampled by many neighbouring + // directions of the 16k half-sphere, so nearly all the strongest entries belong to the same two or + // three rows: 30 raw peaks prune down to only four or five distinct directions in practice, and + // those are the strongest - hence the shortest, densest - rows. When a crystal's densest rows share + // one plane, every direction that survives is coplanar, every triple ReduceResults forms from them + // is degenerate, and the indexer returns no cell at all. Keep walking the same magnitude order for + // a few more directions that are 5 deg clear of everything kept - the weak long axis that closes + // such a cell sits well past where the budget stops. They are appended AFTER the sort, so the + // entries above hold their positions and ReduceResults still forms every triple it formed before; + // the shortlist only gains candidates at its end. Four, because the standard reduction combines + // results[0..8] and a five-vector shortlist leaves most of that window unusable. + constexpr int EXTRA_DISTINCT_DIRECTIONS = 4; + + for (int extra = 0; it != fft_result_map.rend() && extra < EXTRA_DISTINCT_DIRECTIONS; ++it) { + Coord dir = direction_vectors.at(it->second.direction); + + bool distinct = true; + for (const auto &v: ret) { + if (std::fabs(dir * v.Normalize()) > COS_5_DEG) { + distinct = false; + break; + } + } + + if (distinct) { + ret.push_back(dir * it->second.length); + extra++; + } + } + return ret; } diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 20cf5603..52b81559 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -20,6 +20,25 @@ namespace { constexpr float ROT_SUBCELL_VOLUME_RATIO = 1.5f; constexpr float ROT_SUBCELL_FRAC_SLACK = 0.02f; + // How much better a lower-symmetry SETTING of an already-chosen lattice has to index before it is + // taken (see the selection below). A subgroup setting holds fewer cell parameters fixed, so it can + // never index less - only a decisively better fit is evidence that the higher symmetry is wrong. + constexpr float ROT_SUBGROUP_FRAC_RATIO = 1.5f; + + // Order of the lattice point group, so "lower symmetry" is a well-defined comparison + // (gemmi's enum orders Trigonal after Tetragonal, which have 6 and 8 rotations). + int LatticePointGroupOrder(gemmi::CrystalSystem s) { + switch (s) { + case gemmi::CrystalSystem::Monoclinic: return 2; + case gemmi::CrystalSystem::Orthorhombic: return 4; + case gemmi::CrystalSystem::Trigonal: return 6; + case gemmi::CrystalSystem::Tetragonal: return 8; + case gemmi::CrystalSystem::Hexagonal: return 12; + case gemmi::CrystalSystem::Cubic: return 24; + default: return 1; // Triclinic + } + } + // Re-express a primitive hexagonal/trigonal lattice in the conventional hexagonal setting // (a = b, gamma = 120). The Niggli-reduced primitive cell carries the two equal-length axes // at gamma = 60; replacing b with b - a opens that angle to 120 without changing the lattice. @@ -236,7 +255,7 @@ void RotationIndexer::RunIndexing() { // Assemble and select serially, in candidate order - identical to refining them one by one. float best_frac = -1.0f; - float best_vol = 0.0f; + float best_prim_vol = 0.0f; bool have_best = false; size_t best_ci = 0; XtalOptimizerData best_data; @@ -285,14 +304,27 @@ void RotationIndexer::RunIndexing() { // hysteresis margin is overridden by its own primitive. A real superstructure's true // (larger) cell indexes MORE than its sub-cell and is kept by the clearly-more branch; // twin lattices share the cell volume, so this never disturbs twin selection. - const float cand_vol = std::abs(data.latt.CalcVolume()); - const bool clearly_more = frac > best_frac + 0.05f && frac > 0.15f; + // Volumes are compared PRIMITIVE. A centred conventional cell is an exact integer multiple + // of its primitive one, so two settings of the same lattice differ by that factor and + // conventional volumes read a mere change of setting as a sub-cell. And two such settings + // are not comparable on the indexed fraction either: the lower-symmetry one holds fewer + // cell parameters fixed, so it can only index more. An F-cubic lattice contains an + // I-tetragonal cell of the same volume, and refining that cell frees the c/a ratio the + // cubic one holds at sqrt(2), buying back the spots a fraction of a percent of strain had + // put out of tolerance. A slightly higher fraction is therefore no evidence against the + // higher symmetry - only a decisively better fit is. + const float cand_prim_vol = std::abs(data.latt.ToPrimitive(sr.centering).CalcVolume()); + const bool lower_symmetry_setting = have_best + && LatticePointGroupOrder(sr.system) < LatticePointGroupOrder(best_sr.system) + && std::abs(cand_prim_vol - best_prim_vol) < 0.05f * best_prim_vol; + const bool clearly_more = frac > best_frac + 0.05f && frac > 0.15f + && (!lower_symmetry_setting || frac > ROT_SUBGROUP_FRAC_RATIO * best_frac); const bool smaller_subcell = have_best && frac > 0.15f && frac >= best_frac - ROT_SUBCELL_FRAC_SLACK - && cand_vol < best_vol / ROT_SUBCELL_VOLUME_RATIO; + && cand_prim_vol < best_prim_vol / ROT_SUBCELL_VOLUME_RATIO; if (!have_best || clearly_more || smaller_subcell) { best_frac = frac; - best_vol = cand_vol; + best_prim_vol = cand_prim_vol; have_best = true; best_data = std::move(data); best_sr = sr; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 51380139..752fd298 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -817,8 +817,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (!geometry_prepass && prepass_result_.has_value()) { const auto *reuse_sg = prepass_merge_sg_.has_value() ? gemmi::find_spacegroup_by_number(static_cast(*prepass_merge_sg_)) : nullptr; - const double v1 = prepass_result_->lattice.CalcVolume(); - const double v2 = best.result->lattice.CalcVolume(); + // Primitive volumes, like the scheme comparison above: a centred conventional cell is an + // exact integer multiple of its primitive one, so two settings of the SAME lattice differ + // by that factor and comparing the conventional cells reads a mere change of setting as a + // supercell - which then forces pass 1's setting and, with it, its lower symmetry. + const double v1 = std::abs(prepass_result_->lattice + .ToPrimitive(prepass_result_->search_result.centering).CalcVolume()); + const double v2 = std::abs(best.result->lattice + .ToPrimitive(best.result->search_result.centering).CalcVolume()); const bool supercell = v1 > 1.0 && v2 > 1.5 * v1; const bool wrong_centering = reuse_sg != nullptr && best.result->search_result.system != gemmi::CrystalSystem::Triclinic