From 1f27fba1f3aae0393eeb63121188f66e715b51cf Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 10 Oct 2026 11:38:26 +0200 Subject: [PATCH] Model validation expands strict NCS; the written cell follows the model's indexing - ModelValidation reads the model and, where it carries non-given NCS operators (_struct_ncs_oper 'generate' / MTRIX without the given flag), expands them (gemmi expand_ncs, vendored assembly.cpp from GEMMI 0.7.5, MPL-2.0 as the rest of gemmi_gph). Same in the battery's model_check.py and depdata_check.py. An icosahedral-virus entry deposited as 3 of the 45 chains of its ASU was scored against a fifteenth of its atoms: R-free 0.545 -> 0.297. - The indexing the model chooses (AdoptModelFrame) now moves the cell with the indices (reseat_integrated in the main path; CellInBasis in the forecast, twin-file and stills/scale paths). A near-symmetry law (b within 2% of c, beta a degree off 90, a~b in P1) relabelled the reflections only, so the files were written with two axes swapped or beta flipped relative to the indices. Where the law moved the cell, the validation is remade on the moved cell. Measured on 18 open/in-house sets with models: p.hkl and MTZ reflection data bit-identical on all; cell and model R unchanged where the law is exact; the cell now matches the deposited one on the 6 pseudo-merohedral cases (orthorhombic a/c swap, three monoclinic beta flips, a P1 alpha/beta swap), R-free within +-0.001 there. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- common/UnitCell.h | 1 + gemmi_gph/CMakeLists.txt | 4 +- gemmi_gph/assembly.cpp | 466 ++++++++++++++++++ .../structure_refinement/ModelValidation.cpp | 12 + rugnux/RugnuxScaleMerge.cpp | 27 +- rugnux/rugnux_cli.cpp | 6 +- tests/ModelValidationTest.cpp | 40 ++ tests/ReindexAmbiguityTest.cpp | 27 + tools/battery/depdata_check.py | 2 +- tools/battery/model_check.py | 12 +- 10 files changed, 584 insertions(+), 13 deletions(-) create mode 100644 gemmi_gph/assembly.cpp diff --git a/common/UnitCell.h b/common/UnitCell.h index 2e3f37020..787210050 100644 --- a/common/UnitCell.h +++ b/common/UnitCell.h @@ -17,6 +17,7 @@ struct UnitCell { float gamma; operator gemmi::UnitCell() const; + bool operator==(const UnitCell &) const = default; [[nodiscard]] bool is_finite() const; [[nodiscard]] bool is_close(const UnitCell &ref, float dist_tolerance, float angle_tolerance_deg) const; diff --git a/gemmi_gph/CMakeLists.txt b/gemmi_gph/CMakeLists.txt index c8de1d26e..c9c129ecf 100644 --- a/gemmi_gph/CMakeLists.txt +++ b/gemmi_gph/CMakeLists.txt @@ -9,11 +9,13 @@ # # to_mmcif/to_pdb are the writing side, which `rugnux --model` needs to put the model it has placed # against the data on disk. to_mmcif.cpp calls use_hetatm() out of to_pdb.cpp, so the two come together. +# +# assembly.cpp expands a deposited model's strict-NCS operators into the full asymmetric unit. ADD_LIBRARY(gemmi STATIC symmetry.cpp gz.cpp mtz.cpp sprintf.cpp xds_ascii.cpp pdb.cpp resinfo.cpp polyheur.cpp calculate.cpp eig3.cpp ccp4.cpp read_cif.cpp mmcif.cpp json.cpp mmread_gz.cpp - to_mmcif.cpp to_pdb.cpp + to_mmcif.cpp to_pdb.cpp assembly.cpp gemmi/cellred.hpp gemmi/symmetry.hpp gemmi/fail.hpp diff --git a/gemmi_gph/assembly.cpp b/gemmi_gph/assembly.cpp new file mode 100644 index 000000000..343686838 --- /dev/null +++ b/gemmi_gph/assembly.cpp @@ -0,0 +1,466 @@ +// Copyright 2020 Global Phasing Ltd. + +#include "gemmi/assembly.hpp" + +#include // unique_ptr +#include "gemmi/modify.hpp" // transform_pos_and_adp +#include "gemmi/neighbor.hpp" // NeighborSearch + +namespace gemmi { + +namespace { + +bool any_subchain_matches(const Chain& chain, const Assembly::Gen& gen) { + const std::string* prev_subchain = nullptr; + for (const Residue& res : chain.residues) + if (prev_subchain == nullptr || res.subchain != *prev_subchain) { + if (in_vector(res.subchain, gen.subchains)) + return true; + prev_subchain = &res.subchain; + } + return false; +} + +// chain name mappings old->new for each Assembly::Operator used +struct ChainMap { + bool uses_segments = false; + std::string id; + std::map names; +}; +struct AssemblyMapping { + HowToNameCopiedChain how; + // records subchain name correspondence (new->old), used for updating Entity::subchain + std::map sub; + std::vector chain_maps; +}; + +void update_address(AtomAddress& a, const ChainMap& chain_map) { + auto it = chain_map.names.find(a.chain_name); + if (it != chain_map.names.end()) + a.chain_name = it->second; +} + +void remove_cras(Model& model, std::vector& vec) { + // sort in reverse order, so items can be erased without invalidating pointers + std::sort(vec.begin(), vec.end(), [](const CRA& a, const CRA& b) { + return std::tie(a.chain, a.residue, a.atom) > std::tie(b.chain, b.residue, b.atom); + }); + const Atom* prev_a = nullptr; + for (CRA& cra : vec) { + if (cra.atom == prev_a) + continue; + prev_a = cra.atom; + auto atom_idx = cra.atom - cra.residue->atoms.data(); + cra.residue->atoms.erase(cra.residue->atoms.begin() + atom_idx); + if (cra.residue->atoms.empty()) { + auto res_idx = cra.residue - cra.chain->residues.data(); + cra.chain->residues.erase(cra.chain->residues.begin() + res_idx); + if (cra.chain->residues.empty()) { + auto chain_idx = cra.chain - model.chains.data(); + model.chains.erase(model.chains.begin() + chain_idx); + } + } + } +} + +Model make_assembly_(const Assembly& assembly, const Model& model, + HowToNameCopiedChain how, const Logger& logger, + AssemblyMapping* mapping) { + Model new_model(model.num); + ChainNameGenerator namegen(how); + std::map subs = model.subchain_to_chain(); + int counter = 0; + for (const Assembly::Gen& gen : assembly.generators) { + bool all_chains = (!gen.chains.empty() && gen.chains[0] == "(all)"); + for (const Assembly::Operator& oper : gen.operators) { + if (logger.callback) { + std::string note = cat("Applying ", oper.name, " to"); + if (all_chains) + note += " all chains"; + else if (!gen.chains.empty()) + cat_to(note, " chains: ", join_str(gen.chains, ',')); + else if (!gen.subchains.empty()) + cat_to(note, " subchains: ", join_str(gen.subchains, ',')); + logger.note(note); + } + if (!all_chains) { + for (const std::string& chain_name : gen.chains) + if (!model.find_chain(chain_name)) + logger.err("no chain ", chain_name); + for (const std::string& subchain_name : gen.subchains) + if (subs.find(subchain_name) == subs.end()) + logger.err("no subchain ", subchain_name); + } + // chains are not merged here, multiple chains may have the same name + ChainMap chain_map; + if (counter != 0) { + chain_map.uses_segments = (how == HowToNameCopiedChain::Dup); + chain_map.id = std::to_string(counter); + } + for (const Chain& chain : model.chains) { + // PDB files specify bioassemblies in terms of chains, + // mmCIF files in terms of subchains. + bool whole_chain = (all_chains || in_vector(chain.name, gen.chains)); + if (whole_chain || + (!gen.subchains.empty() && any_subchain_matches(chain, gen))) { + // try add a new empty chain, but first figure out the name for it + auto result = chain_map.names.emplace(chain.name, ""); + if (result.second) // insertion happened - generate a new chain name + result.first->second = namegen.make_new_name(chain.name, counter+1); + new_model.chains.emplace_back(result.first->second); + Chain& new_chain = new_model.chains.back(); + // add residues to the chain + for (const Residue& res : chain.residues) + if (whole_chain || in_vector(res.subchain, gen.subchains)) { + new_chain.residues.push_back(res); + Residue& new_res = new_chain.residues.back(); + transform_pos_and_adp(new_res, oper.transform); + if (!new_res.subchain.empty()) { + // change subchain name for the residue + if (how == HowToNameCopiedChain::Short) + new_res.subchain = new_chain.name + ":" + new_res.subchain; + else if (how == HowToNameCopiedChain::AddNumber) + new_res.subchain += new_chain.name.substr(chain.name.size()); + if (mapping) + mapping->sub.emplace(new_res.subchain, res.subchain); + } + if (chain_map.uses_segments) + new_res.segment = chain_map.id; + } + } + } + if (mapping) + mapping->chain_maps.push_back(std::move(chain_map)); + ++counter; + } + } + return new_model; +} + + +void expand_ncs_model_(Model& model, const std::vector& ncs, + HowToNameCopiedChain how, AssemblyMapping* mapping) { + // For HowToNameCopiedChain::Dup we used to set segment="0" for original + // residues. Now it's left blank. Not sure which is better. + size_t orig_size = model.chains.size(); + ChainNameGenerator namegen(model, how); + if (mapping) { // add identities to AssemblyMapping + mapping->chain_maps.emplace_back(); + ChainMap& chain_map = mapping->chain_maps.back(); + // chain_map.id is left empty + for (const Chain& chain : model.chains) { + chain_map.names.emplace(chain.name, chain.name); + for (const ConstResidueSpan& span : chain.subchains()) { + const std::string& sub_id = span.subchain_id(); + mapping->sub.emplace(sub_id, sub_id); + } + } + } + int counter = 0; + for (const NcsOp& op : ncs) + if (!op.given) { + ChainMap chain_map; + chain_map.uses_segments = (how == HowToNameCopiedChain::Dup); + chain_map.id = op.id; + ++counter; + for (size_t i = 0; i != orig_size; ++i) { + model.chains.push_back(model.chains[i]); + Chain& new_chain = model.chains.back(); + const std::string& old_name = model.chains[i].name; + auto result = chain_map.names.emplace(old_name, ""); + if (how != HowToNameCopiedChain::Dup) { + if (result.second) // if insertion happened - generate a new chain name + result.first->second = namegen.make_new_name(old_name, counter); + new_chain.name = result.first->second; + } + for (Residue& new_res : new_chain.residues) { + transform_pos_and_adp(new_res, op.tr); + if (!new_res.subchain.empty()) { + std::string old_subchain = new_res.subchain; + new_res.subchain = new_chain.name + ":" + new_res.subchain; + if (mapping) + mapping->sub.emplace(new_res.subchain, old_subchain); + } + if (chain_map.uses_segments) + new_res.segment = chain_map.id; + } + } + if (mapping) + mapping->chain_maps.push_back(std::move(chain_map)); + } +} + + +void finalize_expansion(Structure& st, const AssemblyMapping& mapping, + double merge_dist, bool expanding_ncs) { + if (merge_dist > 0) + for (Model& model : st.models) { + merge_atoms_in_expanded_model(model, gemmi::UnitCell(), merge_dist); + assign_serial_numbers(model); + } + + // update Entity::subchains + if (!mapping.sub.empty()) + for (Entity& ent : st.entities) { + std::vector new_subchains; + for (const std::string& s : ent.subchains) + for (const auto& new_old : mapping.sub) + if (new_old.second == s) + new_subchains.push_back(new_old.first); + ent.subchains = std::move(new_subchains); + } + + // connections + std::vector new_connections; + for (const Connection& conn : st.connections) + if (expanding_ncs || conn.asu == Asu::Same) { + bool first = true; + for (const ChainMap& chain_map : mapping.chain_maps) { + auto ch1 = chain_map.names.find(conn.partner1.chain_name); + auto ch2 = chain_map.names.find(conn.partner2.chain_name); + if (ch1 != chain_map.names.end() && + ch2 != chain_map.names.end()) { + Connection new_conn = conn; + new_conn.partner1.chain_name = ch1->second; + new_conn.partner2.chain_name = ch2->second; + if (chain_map.uses_segments) + new_conn.partner1.res_id.segment = + new_conn.partner2.res_id.segment = chain_map.id; + if (st.models[0].find_atom(new_conn.partner1) && + st.models[0].find_atom(new_conn.partner2)) { + if (!first) { + cat_to(new_conn.name, '.', chain_map.id); + first = false; + } + new_connections.push_back(new_conn); + } + } + } + } else { + // connections other than 1_555 are lost when making assembly + // if it's needed - get in touch + } + st.connections = std::move(new_connections); + + if (mapping.how == HowToNameCopiedChain::Dup) + return; + + // cispeps + std::vector new_cispeps; + new_cispeps.reserve(st.cispeps.size() * mapping.chain_maps.size()); + for (const CisPep& cispep : st.cispeps) + for (const ChainMap& chain_map : mapping.chain_maps) { + new_cispeps.push_back(cispep); + update_address(new_cispeps.back().partner_c, chain_map); + update_address(new_cispeps.back().partner_n, chain_map); + } + st.cispeps = std::move(new_cispeps); + + // secondary structure - helices + std::vector new_helices; + new_helices.reserve(st.helices.size() * mapping.chain_maps.size()); + for (const Helix& helix : st.helices) + for (const ChainMap& chain_map : mapping.chain_maps) { + new_helices.push_back(helix); + update_address(new_helices.back().start, chain_map); + update_address(new_helices.back().end, chain_map); + } + st.helices = std::move(new_helices); + + // secondary structure - sheets + std::vector new_sheets; + new_sheets.reserve(st.sheets.size() * mapping.chain_maps.size()); + for (const Sheet& sheet : st.sheets) + for (const ChainMap& chain_map : mapping.chain_maps) { + new_sheets.push_back(sheet); + for (Sheet::Strand& strand : new_sheets.back().strands) { + update_address(strand.start, chain_map); + update_address(strand.end, chain_map); + update_address(strand.hbond_atom2, chain_map); + update_address(strand.hbond_atom1, chain_map); + } + } + st.sheets = new_sheets; +} + +} // anonymous namespace + +Model make_assembly(const Assembly& assembly, const Model& model, + HowToNameCopiedChain how, const Logger& logging) { + return make_assembly_(assembly, model, how, logging, nullptr); +} + +void transform_to_assembly(Structure& st, const std::string& assembly_name, + HowToNameCopiedChain how, const Logger& logging, + bool keep_spacegroup, double merge_dist) { + const Assembly* assembly = st.find_assembly(assembly_name); + std::unique_ptr p1_assembly; + if (!assembly) { + if (assembly_name == "unit_cell") { + p1_assembly.reset(new Assembly(pseudo_assembly_for_unit_cell(st.cell))); + assembly = p1_assembly.get(); + } else if (st.assemblies.empty()) { + fail("no bioassemblies are listed for this structure"); + } else { + fail("wrong assembly name, use one of: " + + join_str(st.assemblies, ' ', [](const Assembly& a) { return a.name; })); + } + } + + AssemblyMapping mapping; + mapping.how = how; + AssemblyMapping* mapping_ptr = &mapping; + for (Model& model : st.models) { + model = make_assembly_(*assembly, model, how, logging, mapping_ptr); + mapping_ptr = nullptr; // AssemblyMapping is based only on the first model + } + finalize_expansion(st, mapping, merge_dist, false); + + // Should Assembly instructions be kept or removed? Currently - removing. + st.assemblies.clear(); + + if (!keep_spacegroup) { + st.spacegroup_hm = "P 1"; // maybe it should be empty? + if (assembly_name != "unit_cell") + st.cell = UnitCell(); + } else { + st.cell.images.clear(); + } +} + +Model expand_ncs_model(const Model& model, const std::vector& ncs, + HowToNameCopiedChain how) { + Model model_copy = model; + expand_ncs_model_(model_copy, ncs, how, nullptr); + return model_copy; +} + +void expand_ncs(Structure& st, HowToNameCopiedChain how, double merge_dist) { + AssemblyMapping mapping; + mapping.how = how; + AssemblyMapping* mapping_ptr = &mapping; + for (Model& model : st.models) { + expand_ncs_model_(model, st.ncs, how, mapping_ptr); + mapping_ptr = nullptr; // AssemblyMapping is based only on the first model + } + finalize_expansion(st, mapping, merge_dist, true); + + for (NcsOp& op : st.ncs) + op.given = true; + + st.setup_cell_images(); +} + + +void merge_atoms_in_expanded_model(Model& model, const UnitCell& cell, double max_dist, + bool compare_serial) { + using Mark = NeighborSearch::Mark; + NeighborSearch ns(model, cell, 4.0); + ns.populate(true); + std::vector to_be_deleted; + for (int n_ch = 0; n_ch != (int) model.chains.size(); ++n_ch) { + Chain& chain = model.chains[n_ch]; + for (int n_res = 0; n_res != (int) chain.residues.size(); ++n_res) { + Residue& res = chain.residues[n_res]; + for (int n_atom = 0; n_atom != (int) res.atoms.size(); ++n_atom) { + Atom& atom = res.atoms[n_atom]; + std::vector> equiv; + ns.for_each_cell(atom.pos, [&](std::vector& marks, const Fractional& fr) { + for (Mark& m : marks) { + // We look for the same atoms, but copied to a different chain. + // First quick check that filters out most of non-matching pairs. + if (m.altloc != atom.altloc || m.element != atom.element || + m.chain_idx == n_ch || m.atom_idx != n_atom) + continue; + // Now check if everything else matches. + CRA cra = m.to_cra(model); + if (cra.atom && + (!compare_serial || cra.atom->serial == atom.serial) && + cra.atom->name == atom.name && + cra.atom->b_iso == atom.b_iso && + cra.residue->matches_noseg(res) && + m.pos.dist_sq(ns.grid.unit_cell.orthogonalize(fr)) < sq(max_dist)) + equiv.emplace_back(cra, m.image_idx); + } + }); + if (!equiv.empty()) { + Position pos_sum = atom.pos; + for (auto& t : equiv) { + CRA& cra = t.first; + pos_sum += ns.grid.unit_cell.find_nearest_pbc_position( + atom.pos, cra.atom->pos, t.second); + // The atoms in equiv are to be discarded later. + // Deleting now would invalidate indices in NeighborSearch. + to_be_deleted.push_back(cra); + // Modify the atoms to avoid processing them again. + cra.atom->serial = -1; + cra.atom->name.clear(); + } + size_t n = 1 + equiv.size(); + atom.pos = pos_sum / double(n); + atom.occ = std::min(1.f, n * atom.occ); + } + } + } + } + remove_cras(model, to_be_deleted); +} + + +void shorten_chain_names(Structure& st) { + ChainNameGenerator namegen(HowToNameCopiedChain::Short); + Model& model0 = st.models[0]; + size_t max_len = model0.chains.size() < 63 ? 1 : 2; + for (const Chain& chain : model0.chains) + if (chain.name.length() <= max_len) + namegen.used_names.push_back(chain.name); + for (Chain& chain : model0.chains) + if (chain.name.length() > max_len) + rename_chain(st, chain.name, + namegen.make_short_name(chain.name.substr(0, max_len))); +} + + +static std::vector split_chain_by_segments(Chain& orig, ChainNameGenerator& namegen) { + std::vector chains; + std::vector orig_res; + orig_res.swap(orig.residues); + int n = 0; + for (auto start = orig_res.begin(); start != orig_res.end(); ) { + const std::string& seg = start->segment; + auto ch = std::find_if(chains.begin(), chains.end(), [&](Chain& c) { + return !c.residues.empty() && c.residues[0].segment == seg; }); + if (ch == chains.end()) { + chains.push_back(orig); + ch = chains.end() - 1; + // Naming. Here, Dup means chain name + segment. + switch (namegen.how) { + case HowToNameCopiedChain::Short: + ch->name = namegen.make_short_name(ch->name + seg); + break; + case HowToNameCopiedChain::AddNumber: + ch->name = namegen.make_name_with_numeric_postfix(ch->name, ++n); + break; + case HowToNameCopiedChain::Dup: + ch->name += seg; + break; + } + } + auto end = std::find_if(start, orig_res.end(), + [&](Residue& r) { return r.segment != seg; }); + ch->residues.insert(ch->residues.end(), std::make_move_iterator(start), + std::make_move_iterator(end)); + start = end; + } + return chains; +} + +void split_chains_by_segments(Model& model, HowToNameCopiedChain how) { + ChainNameGenerator namegen(how); + std::vector new_chains; + for (Chain& chain : model.chains) + vector_move_extend(new_chains, split_chain_by_segments(chain, namegen)); + model.chains = std::move(new_chains); +} + +} // namespace gemmi diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index ab8797727..3a3eaaa1b 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -22,6 +22,7 @@ #include #include // read_structure_gz +#include // expand_ncs #include // MaybeGzipped #include // IT92 x-ray form factors #include // DensityCalculator @@ -212,6 +213,15 @@ void reorient_about_centroid(gemmi::Model &model, const gemmi::Mat33 &rot) { SetModelPositions(model, pos); } +// A deposition can hold only part of the asymmetric unit, with the strict-NCS operators that generate +// the rest (mmCIF _struct_ncs_oper code 'generate', PDB MTRIX without the given flag): an icosahedral +// virus deposits a few chains of the dozens in the crystal's ASU. Read as it is, the model is that +// fraction of the crystal, and every R and map made from it describes the fraction. +void expand_strict_ncs(gemmi::Structure &st) { + if (std::any_of(st.ncs.begin(), st.ncs.end(), [](const gemmi::NcsOp &op) { return !op.given; })) + gemmi::expand_ncs(st, gemmi::HowToNameCopiedChain::Short); +} + // Move a model into another cell by keeping its fractional coordinates. Correct only when the two // cells describe the same axes in the same order, which is what the probe below is there to arrange. void refractionalize_into(gemmi::Structure &st, const gemmi::UnitCell &target) { @@ -436,6 +446,7 @@ ModelValidationResult Validate(const std::vector &merged, // falls back to the content when it does not recognise one. A model arrives named however // whoever produced it named it, so the file itself is the better authority. st = gemmi::read_structure_gz(model_path, gemmi::CoorFormat::Detect); + expand_strict_ncs(st); } catch (const std::exception &e) { result.failure_reason = fmt::format("cannot read model {}: {}", model_path, e.what()); logger.Error("Model validation: {}", result.failure_reason); @@ -1817,6 +1828,7 @@ std::vector ModelReferenceIntensities(const std::string &model // falls back to the content when it does not recognise one. A model arrives named however // whoever produced it named it, so the file itself is the better authority. st = gemmi::read_structure_gz(model_path, gemmi::CoorFormat::Detect); + expand_strict_ncs(st); } catch (const std::exception &e) { logger.Error("Model reference: cannot read model {}: {}", model_path, e.what()); return out; diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 4b3b058c0..a2befb19b 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -3727,8 +3727,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { gate = ahead.write.get_future().share()]() mutable { if (const int32_t cards = get_gpu_count(); cards > 0) pin_gpu(1 % cards); - if (!(ahead.indexing_op == gemmi::Op::identity())) + if (!(ahead.indexing_op == gemmi::Op::identity())) { merged = ReindexMergedIntoAsu(merged, ahead.indexing_op, sg, friedel); + cell = CellInBasis(cell, HklOperator(ahead.indexing_op)); + } gemmi::Op to_model = ahead.change_of_basis_op; to_model.tran = {0, 0, 0}; const gemmi::Op cob = to_model.inverse(); @@ -3792,6 +3794,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { result.merge_statistics_text.erase(pos, ambiguity_warning_text.size()); logger.Info("Indexing ambiguity resolved against the model; warning withdrawn"); } + bool indexing_moved_cell = false; if (data_sg.has_value()) { const gemmi::SpaceGroup *adopted = AdoptModelFrame(validation, sm.merged, *data_sg, @@ -3808,13 +3811,14 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // The unmerged export below comes from the integrated observations, not from the merge, // so it follows the same relabeling - otherwise the two files would describe different // indexings of one dataset. The enantiomorph is a label, so only the indexing moves it. + // The cell moves with the indices: an alternative indexing that is only a near-symmetry of + // the lattice (b within 2% of c, beta within a degree of 90) permutes the axes or flips an + // angle, and the written cell has to be the one the written indices are on. const gemmi::Op &op = validation.indexing_op; if (!(op == gemmi::Op::identity())) { - for (auto &io : indexer->GetIntegrationOutcome()) - for (auto &r : io.reflections) { - const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}}); - r.h = h[0]; r.k = h[1]; r.l = h[2]; - } + const UnitCell before = *result.consensus_cell; + reseat_integrated(HklOperator(op)); + indexing_moved_cell = !(*result.consensus_cell == before); merge_to_written = merge_to_written.combine(op); determined_to_written = determined_to_written.combine(op); result.setting_operator = IndexTriplet(determined_to_written); @@ -3826,8 +3830,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // validation is then made again on the relabelled data, so the maps and the placed model // come out on the same axes as the reflections; what it had decided is kept, since the // second one, facing a model already in the data's setting, has nothing left to decide. + // The same where the model's indexing moved the cell (a near-symmetry, above): the first + // validation placed the model in the cell as it was before. const bool remake = config_.reference_data.empty() && validation.model_fits - && !(validation.change_of_basis_op == gemmi::Op::identity()); + && (!(validation.change_of_basis_op == gemmi::Op::identity()) || indexing_moved_cell); const bool forecast_held = forecast_kept(validation); if (ahead.run.valid() && !forecast_held) { ahead.write.set_value(false); @@ -3843,7 +3849,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { relabel_output(to_model.inverse(), "MODEL", &ahead.merged); ahead.log.ReplayInto(logger); } else { - relabel_output(to_model.inverse(), "MODEL"); + if (!(to_model == gemmi::Op::identity())) + relabel_output(to_model.inverse(), "MODEL"); const auto sg_now = experiment_.GetGemmiSpaceGroup(); remade = ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path, config_.output_prefix, logger, @@ -3899,8 +3906,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { twin_file->space_group, experiment_.GetScalingSettings().GetMergeFriedel(), logger); - if (!(op == gemmi::Op::identity())) + if (!(op == gemmi::Op::identity())) { twin_file->law = op.inverse().combine(twin_file->law).combine(op); + twin_file->cell = CellInBasis(twin_file->cell, HklOperator(op)); + } result.twin_subgroup_merge->twin_law = twin_file->law.as_hkl().triplet(); result.twin_subgroup_merge->space_group = twin_file->space_group.xhm(); } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 4fb623ca4..9a8ea0e21 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -52,6 +52,7 @@ #include "../image_analysis/SensorAbsorption.h" #include "../image_analysis/UpdateReflectionResolution.h" #include "../image_analysis/WriteReflections.h" +#include "../image_analysis/scale_merge/CrystalSetting.h" // CellInBasis, HklOperator #include "../image_analysis/scale_merge/Merge.h" #include "../image_analysis/scale_merge/RfreeFlags.h" #include "../image_analysis/scale_merge/ReindexAmbiguity.h" // MatchReference @@ -2221,13 +2222,16 @@ static int RunRugnux(int argc, char **argv) { // The unmerged export below writes the integrated observations, so it follows the same // relabeling as the merged set - only the indexing moves a reflection; the enantiomorph // is a label - or the two files would describe different indexings of one dataset. + // The cell moves with the indices too, for an indexing that is only a near-symmetry. const gemmi::Op &op = validation.indexing_op; - if (!(op == gemmi::Op::identity())) + if (!(op == gemmi::Op::identity())) { for (auto &io : reflections) for (auto &r : io.reflections) { const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}}); r.h = h[0]; r.k = h[1]; r.l = h[2]; } + experiment.SetUnitCell(CellInBasis(*experiment.GetUnitCell(), HklOperator(op))); + } } // The coordinates that go with the maps just written, in the frame the reflection files // beside them carry. Written whenever the validation ran, not only where the rigid-body diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index a7e0eb588..47720f21a 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -134,6 +134,46 @@ TEST_CASE("ModelValidation_ReadsPdbAndMmcif", "[ModelValidation]") { std::filesystem::remove(misnamed); } +// A deposition can hold part of the asymmetric unit and the strict-NCS operators that generate the rest +// (MTRIX without the given flag). The model it describes is the expanded one: here the two atoms and +// their copy half a cell along a, which the explicit four-atom model spells out. +TEST_CASE("ModelValidation_ExpandsStrictNcs", "[ModelValidation]") { + Logger logger("ModelValidation_ExpandsStrictNcs"); + + const char *with_ncs = + "CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n" + "MTRIX1 1 1.000000 0.000000 0.000000 0.00000 1\n" + "MTRIX2 1 0.000000 1.000000 0.000000 0.00000 1\n" + "MTRIX3 1 0.000000 0.000000 1.000000 0.00000 1\n" + "MTRIX1 2 1.000000 0.000000 0.000000 20.00000\n" + "MTRIX2 2 0.000000 1.000000 0.000000 0.00000\n" + "MTRIX3 2 0.000000 0.000000 1.000000 0.00000\n" + "ATOM 1 CA GLY A 1 10.000 12.000 14.000 1.00 20.00 C\n" + "ATOM 2 CB GLY A 1 12.000 14.000 16.000 1.00 20.00 C\n" + "END\n"; + const char *expanded = + "CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n" + "ATOM 1 CA GLY A 1 10.000 12.000 14.000 1.00 20.00 C\n" + "ATOM 2 CB GLY A 1 12.000 14.000 16.000 1.00 20.00 C\n" + "ATOM 3 CA GLY B 1 30.000 12.000 14.000 1.00 20.00 C\n" + "ATOM 4 CB GLY B 1 32.000 14.000 16.000 1.00 20.00 C\n" + "END\n"; + const auto ncs_path = WriteTemp("model_validation_test_ncs.pdb", with_ncs); + const auto full_path = WriteTemp("model_validation_test_ncs_full.pdb", expanded); + + const auto from_ncs = ModelReferenceIntensities(ncs_path, {}, {}, 4.0, logger); + const auto from_full = ModelReferenceIntensities(full_path, {}, {}, 4.0, logger); + REQUIRE(from_ncs.size() == from_full.size()); + REQUIRE_FALSE(from_ncs.empty()); + for (size_t i = 0; i < from_ncs.size(); i++) { + CHECK(from_ncs[i].h == from_full[i].h); + CHECK(from_ncs[i].I == Catch::Approx(from_full[i].I).margin(1e-3 * from_full[i].I + 1e-6)); + } + + std::filesystem::remove(ncs_path); + std::filesystem::remove(full_path); +} + // A model that cannot be used must say why. Returning an empty result and logging was enough to lose // the fact entirely: the run finished successfully with no R-free and no maps, which is exactly what // a run without --model looks like. diff --git a/tests/ReindexAmbiguityTest.cpp b/tests/ReindexAmbiguityTest.cpp index 9ab69ef46..e092c464d 100644 --- a/tests/ReindexAmbiguityTest.cpp +++ b/tests/ReindexAmbiguityTest.cpp @@ -3,10 +3,12 @@ #include +#include #include #include #include +#include "../image_analysis/scale_merge/CrystalSetting.h" #include "../image_analysis/scale_merge/HKLKey.h" #include "../image_analysis/scale_merge/ReindexAmbiguity.h" @@ -212,3 +214,28 @@ TEST_CASE("Reindex: every ambiguity operator is a proper rotation", "[reindex]") CHECK(op.det_rot() == gemmi::Op::DEN * gemmi::Op::DEN * gemmi::Op::DEN); } } + +TEST_CASE("Reindex: a near-symmetry law moves the cell with the indices", "[reindex]") { + // Pseudo-merohedry: b within 2% of c in a primitive orthorhombic cell, and beta a degree off 90 in + // a monoclinic one. Each admits an indexing law that is not a symmetry of the cell, so relabelling + // the reflections by it permutes the axes (or flips beta), and the cell they are written with has + // to follow - every reflection then keeps its own d-spacing. + const std::vector> cases{ + {UnitCell{50.0f, 60.0f, 61.2f, 90.0f, 90.0f, 90.0f}, 16}, // P222, b ~ c + {UnitCell{70.0f, 79.0f, 83.0f, 90.0f, 91.1f, 90.0f}, 3}, // P2, beta ~ 90 + }; + for (const auto &[cell, sg] : cases) { + const auto laws = ReindexAmbiguityOperators(cell, SG(sg)); + REQUIRE_FALSE(laws.empty()); + bool moved = false; + for (const gemmi::Op &law : laws) { + const UnitCell written = CellInBasis(cell, HklOperator(law)); + moved = moved || std::abs(written.b - cell.b) > 0.5f || std::abs(written.beta - cell.beta) > 0.5f; + const gemmi::UnitCell before = cell, after = written; + for (const gemmi::Op::Miller hkl : {gemmi::Op::Miller{1, 2, 3}, gemmi::Op::Miller{-4, 0, 7}, + gemmi::Op::Miller{5, -3, 1}}) + CHECK(after.calculate_d(law.apply_to_hkl(hkl)) == Catch::Approx(before.calculate_d(hkl)).epsilon(1e-5)); + } + CHECK(moved); + } +} diff --git a/tools/battery/depdata_check.py b/tools/battery/depdata_check.py index 0a1759c6c..ab828c7fc 100644 --- a/tools/battery/depdata_check.py +++ b/tools/battery/depdata_check.py @@ -134,7 +134,7 @@ def check(mtz_path, pdb_id, cache_dir=model_check.DEFAULT_CACHE, wavelength=None os.path.join(cache_dir, pdb_id + "-sf.cif.gz")) if sf is None: return dict(res, dep_status="n/a", dep_reason="no deposited structure factors") - st = gemmi.read_structure(meta["xyz"]) + st = model_check.read_model(meta["xyz"]) sg_m, cell_m = st.find_spacegroup(), st.cell if sg_m is None: return dict(res, dep_reason=f"model space group '{st.spacegroup_hm}' not recognised") diff --git a/tools/battery/model_check.py b/tools/battery/model_check.py index 1c668b8fe..5212c5c8f 100644 --- a/tools/battery/model_check.py +++ b/tools/battery/model_check.py @@ -404,6 +404,16 @@ def refined_shared(workdir, data, dep, cell_m, sg_m, xyz, twin): return {k_: (round(v, 4) if isinstance(v, float) else v) for k_, v in out.items()} +def read_model(xyz): + """The deposited model with its strict-NCS copies generated: an entry can hold part of the + asymmetric unit and the operators that make the rest (_struct_ncs_oper 'generate', MTRIX without + the given flag) - an icosahedral virus deposits a few chains of dozens.""" + st = gemmi.read_structure(xyz) + if any(not op.given for op in st.ncs): + st.expand_ncs(gemmi.HowToNameCopiedChain.Short) + return st + + def model_for_refmac(st, workdir): """The deposited model written in a form REFMAC reads, unchanged in everything it scores; returns the file and what had to be changed ("none" normally).""" @@ -587,7 +597,7 @@ def check(mtz_path, pdb_id, workdir, cache_dir=DEFAULT_CACHE, baseline=True, ref return done("n/a", "not an X-ray structure: " + ", ".join(meta["methods"])) if meta["rfree"] is None: return done("n/a", "no deposited R-free") - st = gemmi.read_structure(meta["xyz"]) + st = read_model(meta["xyz"]) if len(st) > 1: return done("skipped", f"{len(st)} models in the entry") natoms = st[0].count_atom_sites()