Jungfraujoch/image_analysis/structure_refinement/RefinementNcs.h

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#pragma once

#include <array>
#include <string>
#include <vector>

#include "gemmi/model.hpp"

#include "RefinementRestraints.h"

// Non-crystallographic symmetry for restrained refinement (Refinement.h): the copies of one molecule in the
// asymmetric unit, found from the model alone, and local restraints that hold each copy's interatomic distances to
// the same distances in another copy. Whether they are used, and how tightly, the weight trials decide by R-free.

// The copies of a model and the restraints between them.
struct NcsCopies {
    // Per pair of copies: a distance in one, a-b, against the same distance in the other, c-d (NcsDistanceTerm),
    // with sigma 1 A (NcsRestraintsAt sets it) and the release distance in `ideal`.
    RestraintSet distances;
    std::vector<int> record_residue;              // per record: its residue, numbered over the model in order
    // Per pair of copies: their corresponding C-alpha (P for a nucleic acid) atoms, for the rmsd between them.
    std::vector<std::vector<std::array<int, 2>>> anchors;
    int groups = 0;                               // molecules present more than once
    int copies = 0;                               // the chains of those molecules
    std::string description;                      // the groups as chain names, e.g. "A=B=C D=E"
    bool Empty() const { return distances.records.empty(); }
};

// The copies of `model` (whose entities are set up, as BuildRefinementTopology leaves them): polymer chains whose
// sequences align with more than 15 residues identical, identical in 80% of the longer chain's residues, and whose
// C-alpha (P) atoms are within 2.5 A rms of each other after superposition are copies of one molecule (REFMAC's
// criterion, the identity taken on the longer chain); the same sequence in a conformation further away - a
// quasi-equivalent subunit of a virus capsid - is another molecule. Their residues of the same name aligned to each
// other correspond, and in them the atoms of the same name - neither a hydrogen nor an alternative conformation, all
// refinable (topo.refinable). Each copy is paired with the next of its molecule (the last with the first, where there
// are three or more), and each pair restrained on the distances within a chain shorter than NCS_DISTANCE_MAX in either
// copy, those the dictionary's bonds and angles fix excepted. Atoms are numbered in model order, as ModelPositions()
// lists them.
NcsCopies FindNcsCopies(const gemmi::Model &model, const RefinementTopology &topo);

// The restraints at `sigma` (A).
RestraintSet NcsRestraintsAt(const NcsCopies &ncs, double sigma);

// The rmsd of corresponding C-alpha (P) atoms of two copies after superposition, averaged over the pairs of copies.
double NcsRmsd(const NcsCopies &ncs, const std::vector<std::array<double, 4>> &xyzb);

// The residues let go: those whose restraints are more than half beyond their release distance.
int NcsReleasedResidues(const NcsCopies &ncs, const std::vector<std::array<double, 4>> &xyzb);