docs: refinement chapter (§15), REFINE_ report keys, acknowledgements and in-source credits

- New docs/CPU_DATA_ANALYSIS_REFINEMENT.md (§15): the ML amplitude target and its FFT gradient,
  scale/robust fit/fine solvent mask, D and the variance (work/free), the embedded Monomer Library
  subset and topology, geometry/metal/KL-B restraints, the minimiser, the weight trials, refined maps,
  blobs, GPU/CPU, and a "Not yet" list; with its own References. Linked from TECHNICAL.md and
  CPU_DATA_ANALYSIS.md; §14 no longer says the structure is never refined.
- RUGNUX_REPORT.md: a Refinement (REFINE_*) section - every REFINE_ and REFINE_SCORE_ key, the
  --mode refine input keys and the _blobs.json fields.
- RUGNUX_MODEL_AND_REFERENCE.md: "What it does not do yet"; RUGNUX_ADVANCED.md: --model mentions the
  refinement, --no-refine as in the usage, a --mode refine options table; GPU_DATA_ANALYSIS.md: the
  refinement's CUDA engines and the missing Metal path.
- ACKNOWLEDGEMENT.md: Brunger, Karplus & Petsko 1989 (gradient-balanced starting weight), Nocedal 1980
  (two-loop recursion), Charbonnier et al. 1997 (pseudo-Huber), Jiang & Brunger 1994 (flat solvent
  mask), Kullback & Leibler 1951; all DOIs checked against Crossref.
- In-source credit lines: the weight balancing (replacing a "Brunger (1992)-style" line that named the
  R-free paper), Nocedal 1980, Charbonnier 1997, Jiang & Brunger 1994, Emsley et al. 2010 at
  FindDensityBlobs, Read 1986 at the refined maps. Comments only.
- docs/THIRD_PARTY_NOTICES.md regenerated from the root file with update_version.sh's own step (it
  lacked the CCP4 Monomer Library row).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit 09375dba351bae0cc73e9a6f0f79455dc92b3a10)
This commit is contained in:
leonarski_f committed 2026-10-26 11:07:43 +01:00
1 parent f960bd7a80
commit e406eb17ea
18 files changed
+496 -21

No files matched your search

+19 -6
View File
@@ -418,7 +418,11 @@ observed scale with an overall factor, an anisotropic B and a flat bulk-solvent
model of A. Fokine and A. Urzhumtsev, "Flat bulk-solvent model: obtaining optimal parameters"
(2002), Acta Cryst. D58, 1387-1392
[doi:10.1107/S0907444902010284](https://doi.org/10.1107/S0907444902010284), which is also the source
of the starting values and of the range those two parameters are physically meaningful over. The
of the starting values and of the range those two parameters are physically meaningful over. The mask itself -
the region a probe cannot reach from the atoms, shrunk back - is J.-S. Jiang and A. T. Brünger's, "Protein
hydration observed by X-ray diffraction: solvation properties of penicillopepsin and neuraminidase crystal
structures" (1994), J. Mol. Biol. 243, 100-115
[doi:10.1006/jmbi.1994.1633](https://doi.org/10.1006/jmbi.1994.1633), drawn with the radii REFMAC5 uses. The
procedure that fits them — a grid search over that range for the solvent pair, with the overall
scale and the anisotropic B refitted at every grid point — follows P. V. Afonine,
R. W. Grosse-Kunstleve and P. D. Adams, "A robust bulk-solvent correction and anisotropic scaling
@@ -540,7 +544,9 @@ R. C. Agarwal, "A new least-squares refinement technique based on the fast Fouri
The positions and B factors are minimised by L-BFGS (rugnux's own, its memory kept from one macro-cycle to
the next, the B factors' steps scaled by their estimated curvature): D. C. Liu and J. Nocedal, "On the limited
memory BFGS method for large scale optimization" (1989), Math. Program. 45, 503-528
[doi:10.1007/BF01589116](https://doi.org/10.1007/BF01589116). The geometry restraints are those of the CCP4
[doi:10.1007/BF01589116](https://doi.org/10.1007/BF01589116), whose two-loop recursion is J. Nocedal's, "Updating
quasi-Newton matrices with limited storage" (1980), Math. Comp. 35, 773-782
[doi:10.1090/S0025-5718-1980-0572855-7](https://doi.org/10.1090/S0025-5718-1980-0572855-7). The geometry restraints are those of the CCP4
Monomer Library, a subset of which rugnux carries (see THIRD_PARTY_NOTICES.md): A. A. Vagin,
R. A. Steiner, A. A. Lebedev, L. Potterton, S. McNicholas, F. Long and G. N. Murshudov, "REFMAC5
dictionary: organization of prior chemical knowledge and guidelines for its use" (2004), Acta Cryst. D60,
@@ -552,11 +558,15 @@ A. Blondel and M. Karplus, "New formulation for derivatives of torsion angles an
in molecular mechanics: elimination of singularities" (1996), J. Comput. Chem. 17, 1132-1141
[doi:10.1002/(SICI)1096-987X(19960715)17:9<1132::AID-JCC5>3.0.CO;2-T](https://doi.org/10.1002/(SICI)1096-987X(19960715)17:9%3C1132::AID-JCC5%3E3.0.CO;2-T).
The B factors of bonded atoms, and of atoms in contact, are restrained to each other in the
Kullback-Leibler form REFMAC5 uses, whose tolerance grows with B (Murshudov et al. 2011, above), and the
Kullback-Leibler form REFMAC5 uses, whose tolerance grows with B (Murshudov et al. 2011, above; the divergence itself
is S. Kullback and R. A. Leibler's, "On information and sufficiency" (1951), Ann. Math. Statist. 22, 79-86
[doi:10.1214/aoms/1177729694](https://doi.org/10.1214/aoms/1177729694)), and the
overall isotropic B of the scale is moved into the atoms as phenix.refine moves it (Afonine et al. 2012,
below). The weights of the X-ray term against the geometry and against the B restraints start from the
ratios of the gradients and are chosen among short refinements by R-free with the geometry held within
bounds, as phenix.refine chooses them: P. V. Afonine, R. W. Grosse-Kunstleve, N. Echols, J. J. Headd,
ratios of the gradients, as X-PLOR's crystallographic weight was first set: A. T. Brünger, M. Karplus and
G. A. Petsko, "Crystallographic refinement by simulated annealing: application to crambin" (1989), Acta Cryst.
A45, 50-61 [doi:10.1107/S0108767388009195](https://doi.org/10.1107/S0108767388009195); and they are chosen
among short refinements by R-free with the geometry held within bounds, as phenix.refine chooses them: P. V. Afonine, R. W. Grosse-Kunstleve, N. Echols, J. J. Headd,
N. W. Moriarty, M. Mustyakimov, T. C. Terwilliger, A. Urzhumtsev, P. H. Zwart and P. D. Adams,
"Towards automated crystallographic structure refinement with phenix.refine" (2012), Acta Cryst. D68,
352-367 [doi:10.1107/S0907444912001308](https://doi.org/10.1107/S0907444912001308). R-free itself:
@@ -566,7 +576,10 @@ The model the refinement scores is scaled to the data with the overall scale and
robust least squares - iteratively reweighted, with the pseudo-Huber weight of each reflection's residual -
so that a few reflections the model cannot explain do not set the scale of the rest: P. J. Huber, "Robust
estimation of a location parameter" (1964), Ann. Math. Statist. 35, 73-101
[doi:10.1214/aoms/1177703732](https://doi.org/10.1214/aoms/1177703732).
[doi:10.1214/aoms/1177703732](https://doi.org/10.1214/aoms/1177703732), in the smooth (pseudo-Huber) form of
P. Charbonnier, L. Blanc-Féraud, G. Aubert and M. Barlaud, "Deterministic edge-preserving regularization in
computed imaging" (1997), IEEE Trans. Image Process. 6, 298-311
[doi:10.1109/83.551699](https://doi.org/10.1109/83.551699).
The bonds of metal ions to their donor atoms that the model's LINK records name are restrained at the
distances of the Monomer Library's own statistics of metal sites (its `metals.json`, embedded as
`monomers/metals.cif`; per metal, donor element and coordination number), in line with the surveys of
+14 -2
View File
@@ -21,7 +21,7 @@ This document describes the crystallographic algorithms implemented in Jungfrauj
15. restrained refinement of the placed model (wherever a model is placed): rigid bodies and group B, then positions and B against a maximum-likelihood amplitude target with the dictionary's geometry, the X-ray weight chosen by R-free within geometry bounds, and the unexplained blobs of the refined difference map.
The reference is split into four parts, in pipeline order; the section numbers run continuously
The reference is split into five parts, in pipeline order; the section numbers run continuously
across them and are the ones the rest of the documentation cites.
- [From images to spots (§0–§3)](CPU_DATA_ANALYSIS_IMAGE.md) — device-side decoding, geometry and
@@ -34,6 +34,10 @@ across them and are the ones the rest of the documentation cites.
- [Space group and validation (§13–§14)](CPU_DATA_ANALYSIS_DECISIONS.md) — the space-group search,
twinning and translational pseudo-symmetry, the resolution cutoff, diffraction anisotropy,
model-based validation and molecular replacement.
- [Refinement of a placed model (§15)](CPU_DATA_ANALYSIS_REFINEMENT.md) — the maximum-likelihood target and
its FFT gradient, the scale, bulk solvent and likelihood parameters, the dictionary and its restraints
(metal sites and B included), the minimiser, the weights chosen by trial, the refined maps and their
unexplained blobs, and what is not implemented yet.
## References
@@ -66,7 +70,15 @@ The methods draw on, and in places reimplement, solutions from:
- A. A. Vagin et al., "REFMAC5 dictionary: organization of prior chemical knowledge and guidelines for its use", *Acta Cryst.* **D60** (2004), 2184-2195 (the restraint dictionary, §15).
- P. V. Afonine et al., "Towards automated crystallographic structure refinement with phenix.refine", *Acta Cryst.* **D68** (2012), 352-367 (the weight chosen among refinements by R-free, §15).
- D. C. Liu & J. Nocedal, "On the limited memory BFGS method for large scale optimization", *Math. Program.* **45** (1989), 503-528 (the minimiser, §15).
- M. Wojdyr, "GEMMI: A library for structural biology", *J. Open Source Softw.* **7** (2022), 4200 (model / structure-factor / map machinery used in §14).
- J. Nocedal, "Updating quasi-Newton matrices with limited storage", *Math. Comp.* **35** (1980), 773-782 (the two-loop recursion of the minimiser, §15.6).
- G. N. Murshudov, P. Skubák, A. A. Lebedev et al., "REFMAC5 for the refinement of macromolecular crystal structures", *Acta Cryst.* **D67** (2011), 355-367 (the likelihood's D from the working set and its variance from the free, and the Kullback-Leibler B restraints, §15.3, §15.5).
- A. T. Brünger, M. Karplus & G. A. Petsko, "Crystallographic refinement by simulated annealing: application to crambin", *Acta Cryst.* **A45** (1989), 50-61 (the starting weight from the ratio of the gradients, §15.7).
- P. J. Huber, "Robust estimation of a location parameter", *Ann. Math. Statist.* **35** (1964), 73-101, and P. Charbonnier, L. Blanc-Féraud, G. Aubert & M. Barlaud, *IEEE Trans. Image Process.* **6** (1997), 298-311 (the robust, pseudo-Huber scale fit of the refinement's score, §15.3).
- J.-S. Jiang & A. T. Brünger, "Protein hydration observed by X-ray diffraction", *J. Mol. Biol.* **243** (1994), 100-115 (the flat bulk-solvent mask, §14.2, §15.3).
- A. Blondel & M. Karplus, *J. Comput. Chem.* **17** (1996), 1132-1141 (the torsion derivatives, §15.5).
- M. M. Harding, *Acta Cryst.* **D57** (2001), 401-411 and **D62** (2006), 678-682, and F. Leonarski, L. D'Ascenzo & P. Auffinger, *Nucleic Acids Res.* **45** (2017), 987-1004 (the metal-site distances, §15.5).
- P. Emsley, B. Lohkamp, W. G. Scott & K. Cowtan, "Features and development of Coot", *Acta Cryst.* **D66** (2010), 486-501 (unmodelled blobs, §15.9).
- M. Wojdyr, "GEMMI: A library for structural biology", *J. Open Source Softw.* **7** (2022), 4200 (model / structure-factor / map machinery used in §14; the topology, riding hydrogens and blob search of §15).
- J. P. Wright, "Experiences with GPU decompression for bitshuffle + LZ4 data", HDF5 User Group meeting (2021), and [github.com/jonwright/bslz4decoders](https://github.com/jonwright/bslz4decoders) (device-side decoding of bitshuffle+LZ4 images, §0).
- A. Thorn & G. M. Sheldrick, "ANODE: anomalous and heavy-atom density calculation", *J. Appl. Cryst.* **44** (2011), 1285-1287 (anomalous difference density read at the model's sites).
- R. Kahn, R. Fourme, A. Gadet, J. Janin, C. Dumas & D. Andre, "Macromolecular crystallography with synchrotron radiation: photographic data collection and polarization correction", *J. Appl. Cryst.* **15** (1982), 330-337 (the azimuthal polarization factor of §2.2, applied to the azimuthal profile, the Bragg intensities and the ring background the beam-stop shadow test compares against).
+2 -2
View File
@@ -1,6 +1,6 @@
# Data analysis: space group and validation (§13–§14)
Part of the [CPU/GPU data-analysis reference](CPU_DATA_ANALYSIS.md); the section numbers are continuous across its four parts.
Part of the [CPU/GPU data-analysis reference](CPU_DATA_ANALYSIS.md); the section numbers are continuous across its parts.
```{contents} On this page
:local:
@@ -124,7 +124,7 @@ Everything lands in `<prefix>_report.txt` section 4 (`ANISOTROPY_*` keys), in th
## 14. Model-based validation: R-free against a model and electron-density maps
Offline (`rugnux --model model.pdb`) the merged data can be scored against a supplied atomic model and electron-density maps computed — enough to confirm that a model fits the data and to inspect the density, not a substitute for refinement. **The structure itself is not refined**; the model is re-fractionalized into the data unit cell and then **placed as one rigid body** (§14.8) - or, where the rigid body does not place it and for a predicted model, by molecular replacement (§14.10) - and the observed amplitudes are the French–Wilson $|F|$ from §10.8, so the R-free and the maps use exactly the same amplitudes as the written reflection file. The model, structure-factor, bulk-solvent and FFT machinery is provided by GEMMI.
Offline (`rugnux --model model.pdb`) the merged data can be scored against a supplied atomic model and electron-density maps computed — enough to confirm that a model fits the data and to inspect the density, not a substitute for refinement. **The structure itself is not refined here** - a placed model is refined afterwards, on the side and deciding nothing (§15, [Refinement of a placed model](CPU_DATA_ANALYSIS_REFINEMENT.md)); the model is re-fractionalized into the data unit cell and then **placed as one rigid body** (§14.8) - or, where the rigid body does not place it and for a predicted model, by molecular replacement (§14.10) - and the observed amplitudes are the French–Wilson $|F|$ from §10.8, so the R-free and the maps use exactly the same amplitudes as the written reflection file. The model, structure-factor, bulk-solvent and FFT machinery is provided by GEMMI.
### 14.1 Model structure factors
+1 -1
View File
@@ -1,6 +1,6 @@
# Data analysis: from images to spots (§0–§3)
Part of the [CPU/GPU data-analysis reference](CPU_DATA_ANALYSIS.md); the section numbers are continuous across its four parts.
Part of the [CPU/GPU data-analysis reference](CPU_DATA_ANALYSIS.md); the section numbers are continuous across its parts.
```{contents} On this page
:local:
+1 -1
View File
@@ -1,6 +1,6 @@
# Data analysis: indexing and geometry (§4–§7)
Part of the [CPU/GPU data-analysis reference](CPU_DATA_ANALYSIS.md); the section numbers are continuous across its four parts.
Part of the [CPU/GPU data-analysis reference](CPU_DATA_ANALYSIS.md); the section numbers are continuous across its parts.
```{contents} On this page
:local:
+1 -1
View File
@@ -1,6 +1,6 @@
# Data analysis: integration, scaling and merging (§8–§12)
Part of the [CPU/GPU data-analysis reference](CPU_DATA_ANALYSIS.md); the section numbers are continuous across its four parts.
Part of the [CPU/GPU data-analysis reference](CPU_DATA_ANALYSIS.md); the section numbers are continuous across its parts.
```{contents} On this page
:local:
+330
View File
@@ -0,0 +1,330 @@
# Data analysis: refinement of a placed model (§15)
Part of the [CPU/GPU data-analysis reference](CPU_DATA_ANALYSIS.md); the section numbers are continuous across its parts.
```{contents} On this page
:local:
:depth: 2
```
## 15. Restrained refinement of a placed model
Wherever a model is placed - it fits where it is (§14.8), or molecular replacement placed it (§14.10) -
rugnux refines it against the merge it was placed on (`--no-refine` skips this). `rugnux --mode refine`
runs the same refinement on a model and merged structure factors that already exist. How to use it, and
what it writes, is in [Models and reference data](RUGNUX_MODEL_AND_REFERENCE.md#refinement); its report
keys are in [The results report](RUGNUX_REPORT.md#refinement-refine_).
The refinement comes **after every processing decision and changes none of them**. Nothing it finds is
read by the space-group search, the indexing choice, the resolution cut or the files `<prefix>.mtz`,
`<prefix>_model.cif` and `<prefix>_maps.mtz`; `R_WORK` and `R_FREE` stay those of the model as placed
(§14.3). What it produces is written beside them: the refined model, its map coefficients, the blobs of its
difference map and the `REFINE_` keys.
The library reads and writes no file. The model and the reflections come in, the refined model, its
statistics and its map coefficients go out; the restraint dictionary is compiled into the binary (§15.4).
### 15.1 What is refined, and in what order
The model is taken as the refinement can use it: the first model of the file only, its hydrogens removed,
every anisotropic B replaced by its isotropic equivalent, and every B clamped to the range the refinement
allows (2 Ų up to twice the model's largest B, within 100-300 Ų). Occupancies are kept as they are, and
so are alternative conformations, each refined with its own restraints.
Two stages follow.
1. **Placement.** Each chain of at least 50 non-water atoms is moved as a rigid body - an angle-axis
rotation about its centroid, scaled by its rms radius so that a unit of rotation moves a typical atom by
about an angstrom, and a translation - against the X-ray target (§15.2), by Ceres' L-BFGS for 30
iterations. Then one B shift per residue. Each of the two is **kept only where it lowers R-free by more
than 0.001**: a deposited model was refined atom by atom already, and moving its chains or shifting its
B factors residue by residue then only makes it worse. The placement stage needs no dictionary.
2. **Restrained refinement.** Every atom's position and isotropic B, against the X-ray target and the
geometry of the dictionary (§15.4-§15.5), in weight trials chosen by R-free (§15.7). Residues the
dictionary does not describe, and residues with an atom it does not know, are **held fixed**: they stay
in the structure factors, in the bulk-solvent mask and in the repulsions the refined atoms see, but
they do not move. So do atoms of zero occupancy.
### 15.2 The X-ray target
The target is the maximum-likelihood amplitude target: for an acentric reflection the Rice distribution of
the observed amplitude $F_o$ given the model's $|F_\mathrm{model}|$, for a centric one its centric
counterpart,
$$
-\ln L_\mathrm{acentric} = \ln\Sigma + \frac{F_o^2 + D^2|F_\mathrm{model}|^2}{\Sigma} - \ln I_0\!\left(\frac{2F_o D|F_\mathrm{model}|}{\Sigma}\right),\qquad
-\ln L_\mathrm{centric} = \tfrac12\ln\Sigma + \frac{F_o^2 + D^2|F_\mathrm{model}|^2}{2\Sigma} - \ln\cosh\frac{F_o D|F_\mathrm{model}|}{\Sigma},
$$
(the terms that depend on neither the model nor the parameters left out) with $D$ the fraction of the model that is right and $\Sigma$ the variance that the rest of the structure and
the measurement leave (§15.3), summed over the working reflections. $\ln I_0$ and $I_1/I_0$ are evaluated
with the polynomial approximations of Abramowitz and Stegun, through the exponentially scaled form for large
arguments so that nothing overflows. The amplitudes are the French-Wilson $|F|$ of §10.8 - the ones the
written reflection file carries - and the free reflections are those of §10.7.
$F_\mathrm{model}$ is the model of §14.2,
$$
F_\mathrm{model} = k\,e^{-\mathbf{s}^\top\mathbf{B}\mathbf{s}/4}\left(F_\mathrm{calc} + k_\mathrm{sol}\,e^{-B_\mathrm{sol}s^2/4}\,F_\mathrm{mask}\right),
$$
with the atoms' $F_\mathrm{calc}$ computed with the dictionary's riding hydrogens included (§15.4). The
scale, the anisotropic B, the bulk solvent and the likelihood's parameters are fitted at the start of each
**macro-cycle** (§15.3) and held fixed while the minimiser runs; the mask is drawn for the model as it is
at that point.
**The gradient by the FFT method.** The atoms' density is placed on a grid with a Gaussian blur (the one
REFMAC-compatible density calculation uses, chosen for the lowest B a macro-cycle can reach) and transformed
once; the crystal's symmetry is composed in reciprocal space from the transform of one copy, as for the rigid
body (§14.8). The derivative of the target with respect to every structure factor, taken back to real space
by one inverse transform through the adjoint of that composition, is a map; each atom's derivatives with
respect to $x, y, z$ and $B$ are that map summed against the derivatives of its own blurred density. One
evaluation of the target and of the gradient of every atom thus costs a density gather, two FFTs and a
gradient gather, whatever the number of atoms. The derivative of a riding hydrogen is passed to its parent.
### 15.3 Scale, bulk solvent and the likelihood parameters
**Scale and solvent.** The flat bulk-solvent mask is drawn with the REFMAC radii (probe 1.0 Å, shrink 0.8 Å)
on a grid of its own, with a spacing of $d_\mathrm{min}/3$ or 0.4 Å, whichever is finer (coarser only where
a cell would need more than $2^{26}$ points). On the atoms' density grid the shrink is lost at low
resolution: at 3 Å its 1 Å step leaves no grid point within the shrink radius. $k_\mathrm{sol}$ and
$B_\mathrm{sol}$ are searched over their physical range on the working reflections, and the overall scale
and anisotropic B are fitted at each point of that search (§14.2) - at the points of its fine pass by
**robust least squares**:
iteratively reweighted (four reweightings), each reflection's residual $r = F_o - |F_\mathrm{model}|$
weighted by $1/\sqrt{1 + (r/c)^2}$ with $c$ the median $|r|$ of the unweighted fit. The fit then minimises
the pseudo-Huber loss $\sum 2c^2\left(\sqrt{1 + (r/c)^2} - 1\right)$, quadratic for the bulk of the
reflections and linear for the few a model cannot explain - a beam-stop shadow, an ice ring, a missing part
of the model - which ordinary least squares would let set the scale of all the rest. When the B factors
are refined, the overall isotropic part of the scale's B is moved into the atoms (and into $B_\mathrm{sol}$,
so that $F_\mathrm{model}$ stays as it is); left in the scale, it and the atoms' B drift against each other
as the B factors are refined.
With these three choices - riding hydrogens, the fine mask grid and the robust scale - the R-work and
R-free of a model as given (`REFINE_R_WORK_START`, `REFINE_R_FREE_START`) agree with REFMAC 5.8's 0-cycle
values on the same model and data within 0.002 on twelve of fourteen depositions from 1.1 to 3.6 Å, and
within 0.006 on the other two (`tools/battery/refine_parity.py`).
**$D$ and $\Sigma$.** The two parameters of the likelihood are fitted by maximum likelihood in resolution
bins of the free set (100 free reflections to a bin, at most 30 bins; where the free set has fewer than
100 reflections, the working ones join it so that the likelihood still has a scale), and interpolated
linearly in $1/d^2$ between the bins' centres:
- $\beta$ on the **free** reflections, as the whole spread of $F_o$ about $D F_\mathrm{model}$ with variance
$\varepsilon\beta$ - the measurement error included;
- $D$ then on the **working** reflections of the bin, with the variance
$\Sigma = \varepsilon\beta + 2\sigma^2$ for an acentric reflection and $\varepsilon\beta + \sigma^2$ for a
centric one, $\sigma$ the amplitude's own error.
$\Sigma$ is also the variance of each reflection in the target. Where the measurement error is a large part
of the spread - the weak high-resolution shells - it is counted twice, and $D$ falls towards 0 there, with
the figure of merit. This division of the two parameters between the working and the free set is the one
REFMAC5 describes; measured on REFMAC 5.8's own output it reproduces its $D$ by bin and its overall figure of
merit (the mean over the working reflections) to about 0.01 on most depositions.
### 15.4 The dictionary and the topology
The restraints are those of the **CCP4 Monomer Library**, a subset of which is compiled into rugnux
(`image_analysis/structure_refinement/monomers/`, regenerated by `update_monomers.py`; nothing is read from
disk, and `$CLIBD` is never consulted): the 20 standard amino acids and MSE; the DNA and RNA nucleotides
(DA, DC, DG, DT, A, C, G, U); ATP, ADP, AMP, GTP, GDP and GMP; DMSO; the PEG fragments PEG, PGE, PG4 and 1PE;
glycerol, ethylene glycol, MPD, acetate, sulfate and phosphate; water; and the ions Na, K, Mg, Ca, Mn, Fe²⁺,
Fe³⁺, Co, Ni, Cu, Zn, Cd, Cl, Br and I - 62 compounds - with the library's links and modifications and the
energy types those compounds use.
The model's topology is built from it by GEMMI (`prepare_topology`): the links between residues (the peptide
and cis-peptide, the nucleic-acid backbone, disulfides, and the model's LINK records where the library
describes them), the modifications those links imply, and per residue the bonds, angles, torsions, chiral
volumes and planes. Torsions are restrained where the dictionary gives them an esd, except the main-chain
φ/ψ and the sugar puckers it lists as alternatives. The dictionary's **hydrogens** are added to every
residue it describes (waters excepted), placed by GEMMI from the dictionary's geometry and moved with their
parents: they take part in the structure factors and in no restraint, and are placed again on their
parents at every macro-cycle. The refined model is written with them.
### 15.5 Geometry, metal and B restraints
Every restraint is a squared normalised deviation, $E = ((x - x_\mathrm{ideal})/\sigma)^2$, with its
analytic derivatives:
- **bonds and angles**, ideal and esd from the dictionary;
- **torsions**, modulo the dictionary's period, with derivatives in the form of Blondel and Karplus, free of
the singularities of the arccos form;
- **chiral volumes** $(\mathbf{a}_1 - \mathbf{c})\cdot((\mathbf{a}_2 - \mathbf{c})\times(\mathbf{a}_3 - \mathbf{c}))$,
$\sigma = 0.2$ ų;
- **planes**: the distances of the atoms from their least-squares plane (the plane's normal the smallest
eigenvector of the atoms' scatter, by Jacobi rotations); the plane is optimal, so the derivative of each
distance carries no term from the plane moving;
- **repulsions** of every pair of atoms closer than its contact distance, symmetry mates included, that no
bond, angle or 1-4 relation joins and that are not alternative conformations of each other, with at least
one of the two refinable: $E = ((d_0 - d)/0.2\,\text{Å})^2$ below the contact distance $d_0$ and nothing
above it, $d_0$ the sum of the energy types' van der Waals radii, 2.6 Å where both atoms can hydrogen-bond
and 1.9 Å between a metal and N, O or S. The list is rebuilt at every macro-cycle.
**Metal sites.** A metal ion is restrained to the donor atoms (O, N or S) that the model's LINK /
`_struct_conn` records bind it to. Where the library has a link for the pair - its ZN-CYS, ZN-HISND,
ZN-HISNE, FE-CYS, FE-HISND and FE-HISNE - that link is used as it is, with the modification it makes to
the donor's residue (the proton the bond replaces). Every other record becomes a link of one bond at the
distance of the library's own statistics of metal sites (its `metals.json`, carried as `monomers/metals.cif`)
for the metal, the donor element and the site's coordination number - the number of such records the ion
has - interpolated linearly between two tabulated coordination numbers (rows resting on at least 20 sites
count; σ is the row's standard deviation, at least 0.02 Å). Only the bond is restrained, not
the angles at the metal. A contact without a record stays a repulsion, and a record to a symmetry mate is
left out, as every bond to a symmetry mate is. The rmsd of these bonds is reported apart from the covalent
bonds' (`REFINE_RMSD_METAL_BONDS`).
**B restraints.** The B factors of two atoms are restrained to each other in proportion to their size, by
the symmetric Kullback-Leibler divergence of two isotropic Gaussians of variances $B_i$ and $B_j$:
$$
E = \frac{(B_i - B_j)^2}{\kappa^2 B_i B_j},
$$
a difference restraint whose $\sigma$ is $\kappa\sqrt{B_iB_j}$, so that atoms of high B may differ more than
atoms of low B. $\kappa$ is 0.1 for bonded main-chain atoms and 0.133 for main-chain atoms bonded to a
common atom, 0.2 and 0.3 for the side chains and the rest, and 0.5 for atoms within 4 Å of each other (waters
among them, symmetry mates included) that no bond or angle relates. No atom's B goes below 2 Ų or above
the ceiling of §15.1.
### 15.6 The minimiser
Positions and B factors are minimised in turn, each by an L-BFGS of rugnux's own (`RefinementMinimiser`),
with the two-loop recursion of Nocedal and of Liu and Nocedal, a memory of 20 pairs and a backtracking line
search from the full step to sufficient decrease, its first trial moving no position by more than 0.3 Å and
no B by more than 10 Ų. What sets it apart from a textbook L-BFGS:
- **Its memory is kept from one macro-cycle to the next.** A macro-cycle is a few dozen iterations on a target
whose scale and likelihood parameters were refitted in between; starting each from what the earlier ones
learnt of the curvature, instead of from a steepest-descent step, is what lets two macro-cycles of a
weight trial get as far as many restarted ones. Each weight trial starts with an empty memory.
- **The B steps are scaled by their curvature.** The initial inverse Hessian of the B factors is a
diagonal: the X-ray target's curvature in each atom's B - its shape from the likelihood's curvature in
$|F_\mathrm{model}|$ and the atom's own scattering, summed in resolution shells, calibrated once per trial
against central differences of the gradient along a random $\pm1$ direction - plus the B restraints'
curvature, each element at least a thousandth of the median. The positions take none: their restraints'
curvature is that of stiff bonds and angles, which move atoms together, and as a diagonal it slows the soft
(torsion and X-ray) directions rather than scaling them.
The diagonal is rescaled at each iteration by the latest pair (Oren-Shanno scaling in that metric). A solve
stops after its iteration budget (30 for the positions and 20 for the B factors per macro-cycle) or where an
iteration lowers the target by less than $10^{-7}$ of it.
### 15.7 The weights, decided by trial
The refined target is $w\,T_\mathrm{X\text{-}ray} + E_\mathrm{geometry}$ for the positions and
$w_B\,T_\mathrm{X\text{-}ray} + E_B$ for the B factors. Rather than fixing $w$ and $w_B$ in advance, rugnux
runs **short refinements at several weights and keeps the one that does best on the free set** - a weight
chosen late, on the evidence, instead of early at a rule.
The starting weights balance the gradients: $w_0 = |\nabla E_\mathrm{geometry}| / |\nabla T_\mathrm{X\text{-}ray}|$,
the geometry's gradient taken on the model shaken by 0.1 Å (a refined model's geometry would give none), and
$w_{B,0}$ likewise from the B restraints on the B factors shaken by 20%. Each **trial** is two macro-cycles
(the scale, solvent and likelihood refitted, then the positions, then the B factors) from where it starts.
A trial is **within bounds** where its rmsd from the dictionary is at most 0.02 Å for bonds and 2° for
angles; one trial **displaces** another where it is within bounds and its R-free is lower by more than 0.001
(differences smaller than that do not decide, and the trial listed first - on a ladder, the tighter
geometry - is kept). The trials
are:
1. **The ladder from the placed model**: $w = w_0 \times \{1/1024, 1/128, 1/16\}$, each with the B factors
held and refined at $w_{B,0}$ - six trials.
2. **The trajectory**: up to eight steps, each a set of trials from where the last step ended - the weights
as they are, and in turn the positions' weight halved and doubled (halved because the geometry loosens at
one weight as the model improves, until it leaves the bounds) or the B factors' weight divided and
multiplied by four. A step is taken while it lowers R-free or R-work by more than 0.001 and its R-free
stays within 0.001 of the best so far, and the trajectory ends at its best step. On a deposited model,
whose free set the depositor's refinement has already fitted, R-work falls while R-free rises, and the
trajectory ends where it began.
3. **The final ladder**: the trajectory's weights times $\{0.5, 1, 2\}$, kept where one displaces it.
$w$ never exceeds $w_0/8$. The trials of one set run side by side - on one GPU, on several, or on the CPU's
threads - each from the same start and taking nothing from the one before it, so the result does not depend
on how many run at once.
The refinement can be cut short from outside (the broker's deadline or cancel, Ctrl+C): the check is made
between the minimiser's iterations and before every trial, the trials that finished are chosen among as
usual, and the report says so. Cut before the first trial, the model is the placement's.
### 15.8 Maps of the refined model
The refined model's map coefficients, `<prefix>_refined_maps.mtz`, are $2mF_o - DF_c$ and $mF_o - DF_c$ as in
§14.4, but with $m$ and $D$ those of the refinement's own likelihood (§15.3) rather than the shell
$\sigma_A$ of §14.4, so the maps carry the weights the target was minimised with; $m$ is
$I_1(X)/I_0(X)$, $X = 2F_oD|F_\mathrm{model}|/\Sigma$, for an acentric reflection and $\tanh(X/2)$ for a
centric one. Where the merge kept the Bijvoet pairs, the anomalous difference map is made on the refined
model's phases too, and read by the blob finder below.
### 15.9 Unexplained density: blobs
The $mF_o - DF_c$ map of the refined model is searched for what the model leaves unexplained: the connected
regions above 3 σ that peak above 4 σ, found by flood fill (GEMMI's `find_blobs_by_flood_fill`, the method
of Coot's "unmodelled blobs"), each moved to the copy nearest the model and measured - peak, volume, the
density integrated above the cutoff, the nearest atom, the polar (N, O) atoms at 2.4-3.4 Å and the N, O or S
at 1.8-2.6 Å, and the anomalous map at its peak. Each gets one reading, by rules on what is measured:
| Reading | Rule |
|---|---|
| `MODEL_ERROR` | within 1.6 Å of a model atom: the model's own atom out of place, or a missing alternative conformation |
| `ION` | at least 1.8 Å from the model, and an anomalous peak of at least 4 σ, or two metal-like contacts and a peak of at least 8 σ |
| `LIGAND` | larger than 1.5 times a single atom's blob at this resolution (a sphere of radius $d_\mathrm{min}/2$) and than 8 ų, at least 2.2 Šfrom the model |
| `WATER` | at least 2.2 Å from the model with at least one polar contact |
| `UNASSIGNED` | anything else |
They are listed strongest first (by the integrated density) in `<prefix>_blobs.json`, and the five strongest
in the report. A blob read as `LIGAND` is density of ligand size where the model has nothing: a buffer
molecule or an unmodelled side chain is that as well as the compound that was soaked. The readings point a
person at the density; they do not identify a compound.
### 15.10 GPU and CPU
On a CUDA build the target, the gradient and the restraints are evaluated on the card - the density gather,
cuFFT, the composition of the symmetry and the target per reflection, the adjoint through a CSR list, the inverse transform and one warp per atom for its gradient; the restraints in one kernel whose
per-atom gradient slots are summed in a fixed order; the robust scale fit and the fine solvent mask on the
card as well - with no floating-point atomics, so the result is the same from run to run. The weight trials
of a set run on as many cards as there are, each trial on an engine of its own. Where the case does not fit
the card's budget, or with `--no-gpu`, or on a build without CUDA (macOS included: the refinement has no
Metal path), the CPU engine computes the same target, to rounding, several times more slowly (seconds on a
GPU, typically 2-15 minutes on the CPU).
### 15.11 Not yet
What the refinement does not do yet - each a known gap, not a design choice:
- **Custom ligands and other monomer libraries.** Only the 62 compounds of §15.4 are restrained. Any other
residue is held fixed (it stays in $F_\mathrm{calc}$, the mask and the repulsions); no restraints are
generated for a new ligand, and neither a user's restraint CIF nor an external monomer library can be given.
- **TLS** refinement.
- **Anisotropic ADPs**: an anisotropic B is replaced by its isotropic equivalent, and only isotropic B is
refined.
- **NCS restraints**: copies are refined independently.
- **Twin refinement**: the target is not a twinned one, also where the twin hypothesis' merge was written.
- **Water building**: waters are neither picked nor removed; the blob finder only points at candidates.
- **Alternative conformations and occupancies**: existing conformers are refined in place; no occupancy is
refined and no conformer is added.
- **Rotamer and peptide-flip search**, and **Ramachandran restraints**.
- **Real-space refinement**, and **jelly-body or reference-model restraints**.
- **A Metal (Apple GPU) path**: on macOS the refinement runs on the CPU.
- **Validation beyond geometry**: the report gives the bond and angle rmsd, the metal-bond rmsd and the
number of clashes, not a MolProbity-style analysis (Ramachandran, rotamer and clashscore statistics).
## References
- R. C. Agarwal, "A new least-squares refinement technique based on the fast Fourier transform algorithm", *Acta Cryst.* **A34** (1978), 791-809 (the gradient by the FFT method, §15.2).
- N. S. Pannu & R. J. Read, "Improved structure refinement through maximum likelihood", *Acta Cryst.* **A52** (1996), 659-668, and G. N. Murshudov, A. A. Vagin & E. J. Dodson, "Refinement of macromolecular structures by the maximum-likelihood method", *Acta Cryst.* **D53** (1997), 240-255 (the maximum-likelihood amplitude target and its figure of merit, §15.2, §15.8).
- G. N. Murshudov, P. Skubák, A. A. Lebedev, N. S. Pannu, R. A. Steiner, R. A. Nicholls, M. D. Winn, F. Long & A. A. Vagin, "REFMAC5 for the refinement of macromolecular crystal structures", *Acta Cryst.* **D67** (2011), 355-367 (D from the working set and the variance from the free, §15.3; the Kullback-Leibler form of the B restraints, §15.5).
- R. J. Read, "Improved Fourier coefficients for maps using phases from partial structures with errors", *Acta Cryst.* **A42** (1986), 140-149 (the m, D weighting of the map coefficients, §15.8).
- V. Y. Lunin & T. P. Skovoroda, "R-free likelihood-based estimates of errors for phases calculated from atomic models", *Acta Cryst.* **A51** (1995), 880-887, and V. Y. Lunin, P. V. Afonine & A. G. Urzhumtsev, "Likelihood-based refinement. I. Irremovable model errors", *Acta Cryst.* **A58** (2002), 270-282 (the likelihood's parameters from the free set, §15.3).
- J.-S. Jiang & A. T. Brünger, "Protein hydration observed by X-ray diffraction: solvation properties of penicillopepsin and neuraminidase crystal structures", *J. Mol. Biol.* **243** (1994), 100-115, and A. Fokine & A. Urzhumtsev, "Flat bulk-solvent model: obtaining optimal parameters", *Acta Cryst.* **D58** (2002), 1387-1392 (the flat bulk-solvent mask and its parameters, §15.3).
- P. V. Afonine, R. W. Grosse-Kunstleve & P. D. Adams, "A robust bulk-solvent correction and anisotropic scaling procedure", *Acta Cryst.* **D61** (2005), 850-855 (the grid search over the solvent parameters, §15.3).
- P. J. Huber, "Robust estimation of a location parameter", *Ann. Math. Statist.* **35** (1964), 73-101, and P. Charbonnier, L. Blanc-Féraud, G. Aubert & M. Barlaud, "Deterministic edge-preserving regularization in computed imaging", *IEEE Trans. Image Process.* **6** (1997), 298-311 (the robust scale fit and its smooth, pseudo-Huber loss, §15.3).
- A. A. Vagin, R. A. Steiner, A. A. Lebedev, L. Potterton, S. McNicholas, F. Long & G. N. Murshudov, "REFMAC5 dictionary: organization of prior chemical knowledge and guidelines for its use", *Acta Cryst.* **D60** (2004), 2184-2195 (the restraint dictionary, §15.4).
- M. Wojdyr, "GEMMI: A library for structural biology", *J. Open Source Softw.* **7** (2022), 4200 (the topology, the riding hydrogens and the blob search, §15.4, §15.9).
- A. Blondel & M. Karplus, "New formulation for derivatives of torsion angles and improper torsion angles in molecular mechanics: elimination of singularities", *J. Comput. Chem.* **17** (1996), 1132-1141 (the torsion derivatives, §15.5).
- M. M. Harding, "Geometry of metal-ligand interactions in proteins", *Acta Cryst.* **D57** (2001), 401-411, and "Small revisions to predicted distances around metal sites in proteins", *Acta Cryst.* **D62** (2006), 678-682; F. Leonarski, L. D'Ascenzo & P. Auffinger, "Mg2+ ions: do they bind to nucleobase nitrogens?", *Nucleic Acids Res.* **45** (2017), 987-1004 (the metal-site distances, §15.5).
- J. Nocedal, "Updating quasi-Newton matrices with limited storage", *Math. Comp.* **35** (1980), 773-782, and D. C. Liu & J. Nocedal, "On the limited memory BFGS method for large scale optimization", *Math. Program.* **45** (1989), 503-528 (the minimiser, §15.6).
- A. T. Brünger, M. Karplus & G. A. Petsko, "Crystallographic refinement by simulated annealing: application to crambin", *Acta Cryst.* **A45** (1989), 50-61 (the starting weight from the ratio of the gradients, §15.7).
- P. V. Afonine et al., "Towards automated crystallographic structure refinement with phenix.refine", *Acta Cryst.* **D68** (2012), 352-367 (weights chosen among refinements by R-free within geometry bounds, and the overall B moved into the atoms, §15.3, §15.7).
- A. T. Brünger, "Free R value: a novel statistical quantity for assessing the accuracy of crystal structures", *Nature* **355** (1992), 472-475 (R-free).
- P. Emsley, B. Lohkamp, W. G. Scott & K. Cowtan, "Features and development of Coot", *Acta Cryst.* **D66** (2010), 486-501 (unmodelled blobs, §15.9).
- M. Abramowitz & I. A. Stegun, *Handbook of Mathematical Functions*, National Bureau of Standards (1964), 9.8.1-9.8.4 (the Bessel functions of the target, §15.2).
The full citations, with DOIs, are in [ACKNOWLEDGEMENT.md](ACKNOWLEDGEMENT.md).
+23
View File
@@ -206,6 +206,29 @@ structure factors above; the likelihood scoring runs on the host. The search is
[§14.10](CPU_DATA_ANALYSIS_DECISIONS.md#1410-molecular-replacement-placing-a-model-the-rigid-body-cannot), its keys and decision rules in
[Molecular replacement](RUGNUX_REPORT.md#molecular-replacement).
## Refinement
The restrained refinement of a placed model
([§15](CPU_DATA_ANALYSIS_REFINEMENT.md)) evaluates its target and gradients on the card where its engine fits
the card's budget, and on the CPU otherwise; the two compute the same target, to rounding, and the GPU's
result is the same from run to run (no floating-point atomics).
- **Target and gradient** (`RefinementGPU.cu`) - the atoms' density gathered on the grid, cuFFT, the
crystal's symmetry composed and the likelihood target evaluated per reflection, the derivative map formed
through the adjoint of that composition and transformed back, and each atom's derivatives gathered from it,
one warp per atom. [§15.2](CPU_DATA_ANALYSIS_REFINEMENT.md#152-the-x-ray-target)
- **Restraints** (`RefinementGPU.cu`) - the geometry and B restraints in one kernel, each writing one gradient
slot per atom it involves, summed per atom in a fixed order.
[§15.5](CPU_DATA_ANALYSIS_REFINEMENT.md#155-geometry-metal-and-b-restraints)
- **Robust scale and fine solvent mask** (`ModelScaleGPU.cu`, `ModelMaskGPU.cu`) - the reweighted scale fit and
the mask on its own fine grid, as model validation's engines compute them.
[§15.3](CPU_DATA_ANALYSIS_REFINEMENT.md#153-scale-bulk-solvent-and-the-likelihood-parameters)
The weight trials of one set run side by side, each on an engine of its own, spread over every card of the
machine; the outcome does not depend on how many run at once
([§15.7](CPU_DATA_ANALYSIS_REFINEMENT.md#157-the-weights-decided-by-trial)). The refinement has no Metal
engine: on macOS it runs on the CPU.
## Apple GPUs (Metal)
A macOS build has no CUDA, and three engines run on the Apple GPU instead (none under
+15 -3
View File
@@ -172,7 +172,8 @@ maps. Where the merge kept the Bijvoet split (a rotation merge always does), an
difference map is made too and read, not written: the report names its strongest sites
(`ANOMALOUS_SITE_01`…`10`) and the map's mean height at the model's anomalous scatterers - every
atom from phosphorus up - as one number for the anomalous signal (`ANOMALOUS_SCATTERER_MEAN_SIGMA`).
The structure itself is not refined; the model is re-fractionalized into
The model is not refined for these numbers (a placed model is refined afterwards, on the side - see
[Refinement](RUGNUX_MODEL_AND_REFERENCE.md#refinement)); it is re-fractionalized into
the data cell and then placed as **one rigid body**, so a deposited model from a crystal that is not
quite isomorphous still sits where the density is. The placement runs on the GPU where there is one
and on the CPU otherwise, with nothing to set; the two agree to rounding, not bit for bit.
@@ -355,6 +356,17 @@ Mode — `--mode <name>` (default `mx`):
| `calibration` | Determine the detector geometry from powder rings; writes `<prefix>.poni` and `<prefix>.json` |
| `refine` | Refine the `--model` against merged structure factors (an MTZ or a structure-factor mmCIF) and process nothing; see [Refining a model against structure factors](RUGNUX_MODEL_AND_REFERENCE.md#refining-a-model-against-structure-factors-rugnux---mode-refine) (`--refine-score-only`, `--refine-d-min`) |
Refinement (`--mode refine`): the columns are found, not named - amplitudes with their sigmas (F/SIGF,
FP/SIGFP, ...) where the file has them, else intensities (IMEAN/SIGIMEAN, I/SIGI, ...) made into amplitudes by
French-Wilson, else the mean of the anomalous pairs; the free set from a FreeR_flag-like column (0 the test set,
or the minority value where 0 is the majority, as phenix writes it), drawn by rugnux where the file has none.
The model is taken as placed in the data's cell, and refined in the data's space group.
| Option | Description |
| --- | --- |
| `--refine-score-only` | Score the model as given and refine nothing: R-work/R-free, the scale, the bulk solvent, the likelihood and the geometry, each term under `REFINE_SCORE_` (`--developer` prints them all) |
| `--refine-d-min <num>` | Leave out reflections finer than this (Å) |
Calibration (`--mode calibration`):
| Option | Description |
@@ -470,9 +482,9 @@ Scaling and merging:
| `--scaling-iterations <num>` | Cap on the per-frame scaling iterations; the loop stops earlier once the scales settle, and a merge that hits the cap is flagged `SCALING_NOT_CONVERGED` (default: 100) |
| `-z, --reference <file>` | Reference data of the same crystal form — an MTZ or a PDB structure-factor mmCIF (`-sf.cif`), gzipped or not, recognised by content; `--reference-mtz` is the same option: fixes the space group, its setting and the cell, resolves the [indexing ambiguity](#the-indexing-ambiguity), hands over the R-free set and reports CC<sub>ref</sub>. Not a scale anchor |
| `--reference-column <label>` | Reference column to use, by MTZ label — an mmCIF under the labels `gemmi cif2mtz` gives it (default: auto — F-model, else IMEAN/I/…, else FP/FOBS/F) |
| `--model <file>` | Validate the merged intensities against this atomic model (PDB or mmCIF, gzipped or not; the format is taken from the file's content) — R-work / R-free and maps (see [Validating against a model](#validating-against-a-model-rugnux---model)). An experimental model is placed as a rigid body, and by molecular replacement where the rigid body does not place it; a prediction (AlphaFold or ModelCIF) goes to molecular replacement straight away (`MODEL_PLACEMENT=`). It also settles the frame the reflections are written in: the enantiomorph, and the [indexing ambiguity](#the-indexing-ambiguity) where no `-z` did; another space group or the twin hypothesis' subgroup it chooses is written beside the answer as `<prefix>_mr.mtz`/`.cif`. For serial stills given `-C` / `-S`, the model's structure factors become the per-image reference |
| `--model <file>` | Validate the merged intensities against this atomic model (PDB or mmCIF, gzipped or not; the format is taken from the file's content) — R-work / R-free and maps (see [Validating against a model](#validating-against-a-model-rugnux---model)). An experimental model is placed as a rigid body, and by molecular replacement where the rigid body does not place it; a prediction (AlphaFold or ModelCIF) goes to molecular replacement straight away (`MODEL_PLACEMENT=`). It also settles the frame the reflections are written in: the enantiomorph, and the [indexing ambiguity](#the-indexing-ambiguity) where no `-z` did; another space group or the twin hypothesis' subgroup it chooses is written beside the answer as `<prefix>_mr.mtz`/`.cif`. For serial stills given `-C` / `-S`, the model's structure factors become the per-image reference. A placed model is then [refined](RUGNUX_MODEL_AND_REFERENCE.md#refinement) against the merge it was placed on (`<prefix>_refined_model.cif`, `<prefix>_refined_maps.mtz`, `<prefix>_blobs.json`, `REFINE_` keys) - seconds on a GPU, typically 2-15 minutes on the CPU - which decides nothing about the processing |
| `--mr <file>` | As `--model`, but placed by molecular replacement straight away, without trying the rigid body first. `--mode mx` only; not with `--model` |
| `--no-refine` | Do not refine the placed model (for debugging) |
| `--no-refine` | Do not refine the placed model (for debugging): by default a model that `--model` or `--mr` places is refined |
| `--write-process-h5` | Also write the (large) `_process.h5` when merging (default: only `.mtz`/`.cif`) |
| `--developer` | Write the full `<prefix>_report.txt`: the pipeline-internal keys and the long explanations the default report leaves out. Nothing is computed differently — the same report, rendered in full (see [The results report](RUGNUX_REPORT.md#the-developer-report)) |
| `--finalist-ledger` | Report the full-resolution evidence for each space group the search considered, not only the one it adopted, and the twin-immune zone evidence of the adopted group's 2-folds (report-only; the decision is unchanged). Both are read on a merge of the observations in P1 that only this option makes, so it costs one extra merge; rotation only, and not with a fixed centred `-S` |
+20
View File
@@ -291,6 +291,26 @@ fourteen depositions from 1.1 to 3.6 A this gives REFMAC 5.8's 0-cycle R-work an
within 0.002 on twelve, and within 0.006 on the other two (`tools/battery/refine_parity.py` runs the comparison;
`--refine-score-only` prints each term under `REFINE_SCORE_`, see [`--mode refine`](#refining-a-model-against-structure-factors-rugnux---mode-refine)).
The algorithms - the target and its gradient, the likelihood's parameters, the restraints, the minimiser and
the weight trials - are described in [§15](CPU_DATA_ANALYSIS_REFINEMENT.md), and the keys in
[The results report](RUGNUX_REPORT.md#refinement-refine_).
**What it does not do yet.** The refinement is a first, automatic pass on a placed model, not a replacement
for a refinement program. Not implemented yet:
- custom ligands or a monomer library beyond the embedded one: other residues are held fixed (they stay in
F<sub>calc</sub>, the solvent mask and the clash terms), and no restraints are generated for a new ligand;
- TLS;
- anisotropic ADPs: an ANISOU record is turned into its isotropic B, and only isotropic B is refined;
- NCS restraints;
- twin refinement (a twinned target);
- water building or removal (the blobs only point at candidates);
- occupancy refinement and new alternative conformations (existing ones are refined in place);
- rotamer and peptide-flip search, Ramachandran restraints;
- real-space refinement, jelly-body and reference-model restraints;
- a Metal (Apple GPU) path: on macOS the refinement runs on the CPU;
- validation beyond the geometry's rmsd and clash count: no MolProbity-style report.
## Refining a model against structure factors (`rugnux --mode refine`)
`--mode refine` refines a model against merged structure factors that already exist - processed elsewhere, a
+52
View File
@@ -1112,3 +1112,55 @@ GPU where there is one (`MR_SEARCH_DEVICE=`). `MR_DECISION=` is then:
With `SPACE_GROUP` and `TWIN_SUBGROUP` the placed model is written as `<prefix>_mr_model.pdb`/`.cif`.
Apart from the enantiomorph label, the answer's own files are never changed by the model.
## Refinement (`REFINE_*`)
Where a model was placed and refined ([Refinement](RUGNUX_MODEL_AND_REFERENCE.md#refinement); the algorithms in
[§15](CPU_DATA_ANALYSIS_REFINEMENT.md)), section `5. MODEL VALIDATION` ends with the refinement's keys. In
`--mode refine` they make up the report's only section, `1. REFINEMENT`, after `MODEL_FILE=`, `SPACE_GROUP=`,
`WAVELENGTH=` and what had to be made of the input. **No key outside `REFINE_` reads the refinement**:
`R_WORK=` and `R_FREE=` above stay the model as placed, and so does every decision of the run. Keys marked
*dev* are written with `--developer` only.
| key | meaning |
|---|---|
| `REFINE` | `PERFORMED` \| `NOT_PERFORMED` (with `REFINE_REASON`) |
| `REFINE_DEVICE`, `REFINE_SECONDS` | *dev*. `GPU` or `CPU`, and the time the refinement took |
| `REFINE_D_MIN` | The finest reflection refined against, Å |
| `REFINE_R_WORK_START`, `REFINE_R_FREE_START` | The model as given, scored as the refinement scores (riding hydrogens, B isotropic, the fine solvent mask, the robust scale): REFMAC's 0-cycle R-factors to about 0.002 on the same model and data. Not `R_WORK`/`R_FREE`, which are the validation's own scoring of the placed model |
| `REFINE_RIGID_BODIES` | The chains moved as rigid bodies, where that step was kept (0: not kept) |
| `REFINE_R_WORK_PLACEMENT`, `REFINE_R_FREE_PLACEMENT` | After the placement stage (rigid bodies, then one B per residue, each kept only where it lowers R-free) |
| `REFINE_R_WORK`, `REFINE_R_FREE` | After the restrained stage: the refined model, the one written as `<prefix>_refined_model.cif` |
| `REFINE_RMSD_BONDS`, `REFINE_RMSD_ANGLES` | The refined model's rmsd from the dictionary's bonds (Å) and angles (°); the weight is chosen within 0.02 Å and 2° |
| `REFINE_RMSD_METAL_BONDS` | The rmsd of the metal-donor bonds made from the model's LINK records, apart from the covalent bonds. Present where there are any |
| `REFINE_RMSD_BONDS_START`, `REFINE_RMSD_ANGLES_START` | *dev*. The same for the model as given |
| `REFINE_CLASHES` | Pairs of atoms in the refined model more than 0.4 Å closer than their contact distance (symmetry mates included) |
| `REFINE_XRAY_WEIGHT` | *dev*. The X-ray weight of the positions that was kept |
| `REFINE_ATOMS`, `REFINE_ATOMS_REFINED` | The model's atoms (hydrogens not counted) and those refined: the atoms of the residues the dictionary describes, with an occupancy above zero |
| `REFINE_RESIDUES_FIXED`, `REFINE_RESIDUES_FIXED_NAMES` | The residues the embedded dictionary does not describe, held fixed but kept in the structure factors, the solvent mask and the clash terms, and their names |
| `REFINE_BULK_SOLVENT_K_SOL`, `REFINE_BULK_SOLVENT_B_SOL` | *dev*. The flat bulk solvent of the refined model |
| `REFINE_BLOBS` | The blobs of the refined mFo-DFc map (above 3 σ, peaking above 4 σ) |
| `REFINE_BLOBS_LIGAND`, `REFINE_BLOBS_WATER`, `REFINE_BLOBS_ION`, `REFINE_BLOBS_MODEL_ERROR` | How many of them are read as each ([§15.9](CPU_DATA_ANALYSIS_REFINEMENT.md#159-unexplained-density-blobs)); the rest are unassigned |
A prose paragraph after them lists the five strongest blobs with their reading, peak, volume, position and
nearest atom; all of them are in `<prefix>_blobs.json`, strongest first, one object per blob: `kind`,
`centroid` (Å, Cartesian, the copy nearest the model), `peak_sigma`, `volume_A3`, `score_sigma_A3` (the
density above the cutoff integrated), `nearest_atom`, `nearest_atom_distance_A`, `polar_contacts`,
`metal_contacts` and `anomalous_sigma`.
**The score of the model as given** follows, every term under `REFINE_SCORE_`: `REFINE_SCORE_R_WORK`,
`REFINE_SCORE_R_FREE`, `REFINE_SCORE_N_WORK` and `REFINE_SCORE_N_FREE` by default, and with `--developer` the
overall scale (`REFINE_SCORE_K_OVERALL`), its isotropic and anisotropic B (`REFINE_SCORE_B_ISO`,
`REFINE_SCORE_B_ANISO`: B11 B22 B33 B12 B13 B23 less the isotropic part, Cartesian), the bulk solvent
(`REFINE_SCORE_K_SOL`, `REFINE_SCORE_B_SOL`) and the spacing of its mask grid (`REFINE_SCORE_MASK_SPACING`), the
riding hydrogens (`REFINE_SCORE_HYDROGENS`), the likelihood over the working and the free reflections
(`REFINE_SCORE_NLL_WORK`, `REFINE_SCORE_NLL_FREE`) and its mean figure of merit over the working reflections
(`REFINE_SCORE_FOM`), the geometry's rmsd and rmsZ (`REFINE_SCORE_RMSD_BONDS`, `REFINE_SCORE_RMSZ_BONDS`,
`REFINE_SCORE_RMSD_ANGLES`, `REFINE_SCORE_RMSZ_ANGLES`, `REFINE_SCORE_RMSD_METAL_BONDS`,
`REFINE_SCORE_RMSD_BONDS_ALL` for covalent and metal bonds together) and a table of R-work and R-free in 20
resolution shells. `rugnux --mode refine --refine-score-only` writes this score and nothing else.
In `--mode refine` the input's keys come first: `REFINE_REFLECTIONS` and `REFINE_FREE_REFLECTIONS` (those read,
after any `--refine-d-min`), `REFINE_FRENCH_WILSON` (`TRUE` where the file had intensities only and the amplitudes
were made by French-Wilson), `REFINE_FREE_FLAGS` (`FILE`, or `GENERATED` where the file had no free set and
rugnux drew its own) and `REFINE_D_MIN_LIMIT` (the `--refine-d-min` given).
+1
View File
@@ -13,6 +13,7 @@ CPU_DATA_ANALYSIS_IMAGE
CPU_DATA_ANALYSIS_INDEXING
CPU_DATA_ANALYSIS_INTEGRATION
CPU_DATA_ANALYSIS_DECISIONS
CPU_DATA_ANALYSIS_REFINEMENT
GPU_DATA_ANALYSIS
DETECTOR_GEOMETRY
RUGNUX_BENCHMARKING
+10 -2
View File
@@ -33,7 +33,6 @@ linked into the Jungfraujoch executables.
| [FFTW](https://www.fftw.org/) | 3.3.10 | Matteo Frigo; MIT | GPL-2.0-or-later | [fftw.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/fftw.txt) |
| [Ceres Solver](https://github.com/ceres-solver/ceres-solver) | (pinned) | Google Inc. and contributors | BSD-3-Clause | [ceres-solver.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/ceres-solver.txt) |
| [Abseil](https://github.com/abseil/abseil-cpp) | 20250127 (required by Ceres) | Google Inc. | Apache-2.0 | [abseil.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/abseil.txt) |
| [fast-feedback-indexer](https://github.com/paulscherrerinstitute/fast-feedback-indexer) | (pinned) | PSI | BSD-3-Clause | [fast-feedback-indexer.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/fast-feedback-indexer.txt) |
| [libjpeg-turbo](https://github.com/libjpeg-turbo/libjpeg-turbo) | 3.2.0 | D. R. Commander and others; IJG | IJG + BSD-3-Clause + Zlib | [libjpeg-turbo.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/libjpeg-turbo.txt) |
| [Mbed TLS](https://github.com/Mbed-TLS/mbedtls) (with TF-PSA-Crypto) | 4.2.0 | The Mbed TLS Contributors | Apache-2.0 OR GPL-2.0-or-later (used under Apache-2.0) | [mbedtls.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/mbedtls.txt), [tf-psa-crypto.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/tf-psa-crypto.txt) |
| [Catch2](https://github.com/catchorg/Catch2) | 3.16.0 | Catch2 Authors | BSL-1.0 | [catch2.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/catch2.txt) |
@@ -65,7 +64,9 @@ These live in the source tree (see the path) rather than being fetched; traccc a
| [Bitshuffle (h-perf)](https://github.com/kalcutter/bitshuffle) | `compression/bitshuffle_hperf/` | Kal Conley | Apache-2.0 | [bitshuffle-hperf.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/bitshuffle-hperf.txt) |
| [LZ4](https://github.com/lz4/lz4) | `compression/lz4/` | Yann Collet | BSD-2-Clause | [lz4.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/lz4.txt) |
| [HLS arbitrary-precision types](https://github.com/Xilinx/HLS_arbitrary_Precision_Types) | `fpga/include/` | Xilinx, Inc. | Apache-2.0 | [xilinx-hls-headers.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/xilinx-hls-headers.txt) |
| [fast-feedback-indexer](https://github.com/paulscherrerinstitute/fast-feedback-indexer) (commit a94fd4f) | `image_analysis/indexing/fast-feedback-indexer/` | PSI | BSD-3-Clause | [fast-feedback-indexer.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/fast-feedback-indexer.txt) |
| [GEMMI](https://github.com/project-gemmi/gemmi) | `gemmi_gph/` | Global Phasing Ltd. | MPL-2.0 | [gemmi.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/gemmi.txt) |
| [CCP4 Monomer Library](https://github.com/MonomerLibrary/monomers) (subset: restraint dictionary data) | `image_analysis/structure_refinement/monomers/` | CCP4 / the Monomer Library authors (Vagin et al.) | LGPL-3.0-only | [ccp4-monomer-library.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/ccp4-monomer-library.txt) |
| [PEGTL](https://github.com/taocpp/PEGTL) | `gemmi_gph/gemmi/third_party/tao/` | Dr. Colin Hirsch and Daniel Frey | MIT | [pegtl.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/pegtl.txt) |
| [sajson](https://github.com/chadaustin/sajson) | `gemmi_gph/gemmi/third_party/sajson.h` | Chad Austin | MIT | [sajson.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/sajson.txt) |
| [fast_float](https://github.com/fastfloat/fast_float) | `gemmi_gph/gemmi/third_party/fast_float.h` | The fast_float authors (Daniel Lemire et al.) | Apache-2.0 OR MIT OR BSL-1.0 | [fast-float.txt](https://gitea.psi.ch/mx/jungfraujoch/src/branch/main/licenses/fast-float.txt) |
@@ -130,12 +131,19 @@ served frontend, so the shipped web UI carries its own attribution.
its trust-region Levenberg-Marquardt minimiser, projected line search, polynomial step choice and
sphere manifold for the crystal refinement, following its source. Both files name the origin at the
top and are covered by `licenses/ceres-solver.txt`.
* **fast-feedback-indexer** is also fetched and linked in CUDA builds (table above); separately,
* **fast-feedback-indexer** is also vendored and linked in CUDA builds (table above); separately,
`FFBIDXIndexerCPU.cpp` ports its CUDA kernels to the CPU. The file names the origin at the top and
is covered by `licenses/fast-feedback-indexer.txt`.
* **FFTW** is GPL-2.0-or-later — compatible with, and absorbed by, this project's GPL-3.0 license.
* **Apache-2.0** components: where upstream ships a `NOTICE` file, it is reproduced in the
corresponding `licenses/` text.
* **CCP4 Monomer Library (LGPL-3.0)** is data, not code: the restraint dictionaries of the compounds
rugnux's refinement restrains, copied unmodified per compound into
`image_analysis/structure_refinement/monomers/` (`links_and_mods.cif` and `ener_lib.cif` cut down to
the entries those compounds use, and the library's metal-site distances of `metals.json` for its ions as
`metals.cif`), with the library's own `COPYING`. The build embeds them in the binary
as text, and `monomers/update_monomers.py` regenerates the subset from the upstream library, so the
copy can be replaced by a user with a version of their own.
* **Qt (LGPL-3.0)** and **NVIDIA CUDA (EULA)** carry redistribution conditions beyond a copyright
notice; see their dedicated notice files. The verbatim LGPL-3.0 and CUDA EULA texts are bundled
(`licenses/Qt6-LGPL-3.0.txt`, `licenses/NVIDIA-CUDA-EULA.txt`); the CUDA EULA is the one shipped
@@ -42,6 +42,7 @@ const char *DensityBlobKindName(DensityBlob::Kind kind) {
return "UNASSIGNED";
}
// Following Coot's unmodelled blobs, Emsley et al. (2010) Acta Cryst. D66, 486-501
std::vector<DensityBlob> FindDensityBlobs(const gemmi::Grid<float> &difference, double rms, const gemmi::Structure &model,
double d_min, const gemmi::Grid<float> *anomalous, double anomalous_rms) {
std::vector<DensityBlob> out;
@@ -147,6 +147,7 @@ struct FixedSolventFit {
// a model fits badly (a beam-stop shadow, an ice ring, a missing part of the model), which least squares would
// otherwise let set the scale of all the rest.
// Following Huber (1964) Ann. Math. Statist. 35, 73-101 for the M-estimator, as iteratively reweighted least squares
// The smooth (pseudo-Huber) loss of Charbonnier et al. (1997) IEEE Trans. Image Process. 6, 298-311
void RobustReweight(FixedSolventFit &fit, double &c) {
std::vector<double> r(fit.points.size());
for (size_t i = 0; i < fit.points.size(); i++)
@@ -756,8 +756,8 @@ void RefinePlacement(Context &ctx, bool group_b, RefinementResult &result) {
// The X-ray weight that balances the two gradients: the norm of the geometry's gradient (taken on the model
// shaken by 0.1 A, where a refined model's geometry would have none) over the norm of the X-ray target's.
// Following Brunger (1992)-style gradient balancing as phenix.refine starts its weight search (Afonine et al.
// (2012) Acta Cryst. D68, 352-367)
// Following Brunger, Karplus & Petsko (1989) Acta Cryst. A45, 50-61, as phenix.refine starts its weight search
// (Afonine et al. (2012) Acta Cryst. D68, 352-367)
double BalancedWeight(Context &ctx, const std::vector<int> &refinable) {
std::vector<std::array<double, 4>> gx, gr;
ctx.engine->Evaluate(ctx.Model(), &gx);
@@ -801,6 +801,7 @@ double BalancedADPWeight(Context &ctx, const RestraintSet &b_restraints, const s
// The model's map coefficients and its difference map, weighted by the refinement's own likelihood: m its figure
// of merit and D its D (2mFo - DFc, mFo - DFc).
// Following Read (1986) Acta Cryst. A42, 140-149
void MakeMaps(Context &ctx, const std::vector<MergedReflection> &merged, double wavelength_A, RefinementResult &result) {
const RefinementStage stage("maps");
ctx.Refit();
@@ -15,7 +15,7 @@
// Hessian is a diagonal the caller gives (SetDiagonal: an estimate of the Hessian's diagonal, so that a parameter
// with ten times the curvature of another takes a tenth of its step), scaled each iteration by the latest pair
// (Oren-Shanno in that metric). The line search backtracks from the full step to sufficient decrease.
// Following Liu & Nocedal (1989) Math. Program. 45, 503-528
// Following Nocedal (1980) Math. Comp. 35, 773-782 and Liu & Nocedal (1989) Math. Program. 45, 503-528
class RefinementMinimiser {
public:
// f(x, g): the target at x, and its gradient in g (sized as x); false where it cannot be evaluated there.
@@ -63,6 +63,7 @@ struct XrayRowModel {
// cell too large for that gets the spacing that fills it. The R-factors of low-resolution models settle at a
// step of 0.4 A (within 0.001 of 0.3 A) and move by up to 0.005 between 0.4 and 0.8 A.
// The flat solvent model of Fokine & Urzhumtsev (2002) Acta Cryst. D58, 1387-1392
// The probe-and-shrink mask of Jiang & Brunger (1994) J. Mol. Biol. 243, 100-115
class BulkSolventMask {
public:
static constexpr double MASK_SPACING_MAX = 0.4; // A