leonarski_fandClaude Opus 5.5
4867eda022
A chain is a copy of a molecule only where, besides the sequence rule, its C-alpha (P) atoms lie within 2.5 A rms of the molecule's first copy after superposition (REFMAC5's second criterion): the same sequence in a conformation further away is another molecule. 4NWV (a T=3 capsid expanded from 14 strict-NCS operators) has 45 protein chains of ONE sequence: A, B and C are its quasi-equivalent subunits (A starts 14 residues earlier; A is 3.3-3.6 A from B and C, B and C 1.2 A apart). All 45 were one molecule, paired in model order A-B-C, copies 2.9 A apart, and NCS was refused. Now A x 15 and B/C x 30: NCS kept (sigma 0.1 A), --mode refine R-free 0.2565 -> 0.2450. Every other NCS set (open arm and the MD starts) has all its copies within 2.5 A, so its groups and results are unchanged. Ordering each molecule's copies as a chain of nearest neighbours (strict copies beside each other) was tried too: 4NWV 0.2401, but one MD start of a four-copy set lost 0.039 and a twelve-copy set 0.001 on three starts. Not kept. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37
52 lines
3.1 KiB
C++
52 lines
3.1 KiB
C++
// 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);
|