leonarski_fandClaude Opus 5.5
6472ce3ff1
S1 - the model's displacements kept: an atom's ANISOU is no longer made isotropic; where the model has TLS
groups whose atoms bring no U, their U is made from the groups and the B column read as the total or as the
residual to them, whichever (or the isotropic B) scores the lowest R-free as given. Refined models are written
with ANISOU and TLS groups (mmCIF _pdbx_refine_tls*, and PDB REMARK 3, which gemmi's writer leaves out).
S2/S3 - both engines grid an atom with its full U (U = dev + B/(8 pi^2) I, so the B machinery works on it with
its shape held) and give the target's derivatives in U by the FFT method, CPU and GPU (one warp per atom, no
atomics). After the weights are chosen, the refinement tries, as trials side by side: the B factors refined on;
TLS groups with a residual B (the model's own, one per chain, and chains cut into 2/4/8/16 TLSMD segments by
dynamic programming over a quadratic B(r) fit); and a U per atom with rigid-bond (RIGU form), similarity and
isotropy restraints at three weights. A trial displaces another only within the geometry bounds and with both
R-free (beyond the tie) and the free -ln L lower. The water-build rounds hold the chosen U shapes.
An anisotropic atom's density is now cut at the radius of its largest eigenvalue rather than its largest
diagonal element (also in the validation's density and the GPU gather), which left skewed atoms truncated.
Refine arm + open-arm subset (11 sets, vs 51cbe0ef7), refined R-free: 1exr 0.2064 -> 0.1377 (SHELXL U kept),
193l 0.1689 -> 0.1322 (aniso), 2r8s 0.2262 -> 0.2054 (TLS 47), 1hr2 0.2463 -> 0.2316 (TLS 30), 1k4c 0.2110 ->
0.1881 (TLS 42), 36gk 0.2445 -> 0.2221, 7l6j 0.1766 -> 0.1553, 7os3 0.2040 -> 0.2014, 8agq 0.1948 -> 0.1595,
8sa8 0.1830 -> 0.1548, 6f3p 0.1821 -> 0.1453. Every output U positive definite. The first refinement takes
about 2-3x as long (the hypothesis trials; measured on a shared GPU).
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37
79 lines
5.2 KiB
C++
79 lines
5.2 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/math.hpp"
|
|
#include "gemmi/metadata.hpp" // TlsGroup
|
|
#include "gemmi/model.hpp"
|
|
|
|
// The anisotropic side of restrained refinement (Refinement.h): TLS groups, and the restraints of atoms refined
|
|
// with a U of their own. Atoms are numbered in model order, as everywhere in the refinement.
|
|
|
|
// A group of atoms whose displacement is in part that of a rigid body: U_TLS(r) = T + A L A^T + A S + S^T A^T at
|
|
// r - origin = (x, y, z), A = [[0, z, -y], [-z, 0, x], [y, -x, 0]]. T in A^2, L in rad^2, S in rad A (the files'
|
|
// L and S are in degrees). The trace of S does not enter U: S33 is -(S11 + S22).
|
|
// Following Schomaker & Trueblood (1968) Acta Cryst. B24, 63-76, and Winn, Isupov & Murshudov (2001) Acta Cryst. D57, 122-133
|
|
struct TlsBody {
|
|
std::string id;
|
|
gemmi::Position origin;
|
|
gemmi::SMat33<double> T{0, 0, 0, 0, 0, 0}, L{0, 0, 0, 0, 0, 0};
|
|
gemmi::Mat33 S{0};
|
|
std::vector<int> atoms; // in model order
|
|
|
|
gemmi::SMat33<double> U(const gemmi::Position &pos) const;
|
|
// The 20 parameters (T11 T22 T33 T12 T13 T23, L likewise, S11 S22 S12 S13 S21 S23 S31 S32) and back.
|
|
static constexpr int PARAMETERS = 20;
|
|
std::array<double, PARAMETERS> Get() const;
|
|
void Set(const std::array<double, PARAMETERS> &p);
|
|
// dU/dp at `pos`: per parameter the six U11 U22 U33 U12 U13 U23 (U is linear in the parameters).
|
|
std::array<std::array<double, 6>, PARAMETERS> Derivatives(const gemmi::Position &pos) const;
|
|
// The derivatives of a target in the parameters, from its derivatives in U11 U22 U33 U12 U13 U23 of every atom
|
|
// (`grad_u`, in model order as `atoms`): sum over the body's atoms of grad_u . dU/dp, by the chain rule in matrix
|
|
// form - with G the symmetric matrix of an atom's derivatives (an off-diagonal element half the U_ab one),
|
|
// dT = tr(G dU) gives G for T, A^T G A for L and 2 A^T G for S.
|
|
std::array<double, PARAMETERS> Gradient(const std::vector<gemmi::Atom *> &atoms,
|
|
const std::vector<std::array<double, 6>> &grad_u) const;
|
|
};
|
|
|
|
// The TLS groups of a model's metadata (mmCIF _pdbx_refine_tls*, PDB REMARK 3) with their atoms: the atoms of every
|
|
// residue of a selection's chain from its first residue to its last, waters left out; a group without numbers or
|
|
// without atoms is left out.
|
|
std::vector<TlsBody> TlsBodiesFromModel(const gemmi::Model &model, const std::vector<gemmi::TlsGroup> &groups);
|
|
|
|
// The bodies as the metadata holds them (degrees), each selection one run of consecutive residues of a chain.
|
|
std::vector<gemmi::TlsGroup> TlsBodiesToGroups(const gemmi::Model &model, const std::vector<TlsBody> &bodies);
|
|
|
|
// The restraints of atoms refined with a U of their own (in A^2, U11 U22 U33 U12 U13 U23 per atom):
|
|
// - rigid bond, for each pair of bonded or angle-related atoms (the pairs of the B similarity): the difference of
|
|
// their U along the line between them, |(U_i - U_j) n|^2 / sigma^2, sigma 0.004 A^2 - the components U33, U13 and
|
|
// U23 of the difference in a frame whose z is n, as SHELXL's RIGU restrains them;
|
|
// - similarity, for the same pairs: |U_i - U_j|^2 (Frobenius) / sigma^2, sigma 0.04 A^2 (SHELXL's SIMU);
|
|
// - near isotropy, for an atom in no pair (a water, an ion): |U - tr(U)/3 I|^2 / sigma^2, sigma 0.1 A^2 (ISOR).
|
|
// Following Hirshfeld (1976) Acta Cryst. A32, 239-244, and Thorn, Dittrich & Sheldrick (2012) Acta Cryst. A68, 448-451
|
|
struct URestraints {
|
|
static constexpr double SIGMA_RIGID = 0.004, SIGMA_SIMILAR = 0.04, SIGMA_ISOTROPIC = 0.1;
|
|
std::vector<std::pair<int, int>> pairs;
|
|
std::vector<int> isotropic;
|
|
};
|
|
// The pairs of `b_similarity_pairs` whose atoms both refine a U, and the refined atoms of no pair.
|
|
URestraints MakeURestraints(const std::vector<std::pair<int, int>> &b_similarity_pairs, const std::vector<char> &refines_u);
|
|
// The sum of E at the positions `pos` and the U of every atom `u` (atoms that refine none are not read) and, with
|
|
// `grad` (sized as u), dE/dU per atom added to it.
|
|
double EvaluateURestraints(const URestraints &r, const std::vector<gemmi::Position> &pos,
|
|
const std::vector<std::array<double, 6>> &u, std::vector<std::array<double, 6>> *grad);
|
|
|
|
// Contiguous segments of a chain (TLSMD's idea): residues cut into `k` runs so that the B factors of each run's atoms
|
|
// are best fitted, in least squares, by an isotropic TLS model of its own - B(r) a general quadratic in r (tr(A L A^T)
|
|
// is tr(L)|r|^2 - r^T L r and tr(A S) linear in r, ten parameters as Painter & Merritt's isotropic TLS), by dynamic
|
|
// programming over where the runs end. `residue_atoms`: per residue (in chain order) its atoms, `pos` and `b` per
|
|
// atom. Returns per segment its first residue (the first is 0); fewer than k where the chain has too few residues for
|
|
// `min_residues` each.
|
|
// Following Painter & Merritt (2006) Acta Cryst. D62, 439-450
|
|
std::vector<int> TlsSegments(const std::vector<std::vector<int>> &residue_atoms, const std::vector<gemmi::Position> &pos,
|
|
const std::vector<double> &b, int k, int min_residues);
|