diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 024cb6fe2..4ce120f94 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -242,6 +242,15 @@ reading: A. Thorn and G. M. Sheldrick, "ANODE: anomalous and heavy-atom density J. Appl. Cryst. 44, 1285-1287 [doi:10.1107/S0021889811041768](https://doi.org/10.1107/S0021889811041768). +**sigma_A map coefficients** — the maps written by `--model` are weighted by a maximum-likelihood +sigma_A estimated per resolution shell, giving 2mFo-DFc and mFo-DFc rather than 2Fo-Fc and Fo-Fc. +What is taken is the formalism itself: the Rice and Woolfson likelihoods of |Fo| given |Fc| and +sigma_A, the figure of merit m and the scale D that follow from it, and the result that the +bias-corrected coefficient is 2mFo-DFc for an acentric reflection and mFo for a centric one. +R. J. Read, "Improved Fourier coefficients for maps using phases from partial structures with +errors" (1986), Acta Cryst. A42, 140-149 +[doi:10.1107/S0108767386099622](https://doi.org/10.1107/S0108767386099622). + **Data-quality statistics** follow the established conventions rather than any one program: R_meas and R_pim, CC1/2 and CC\*, and the reporting of I/sigma(I). K. Diederichs and P. A. Karplus, "Improved R-factors for diffraction data analysis in macromolecular crystallography" (1997), Nat. Struct. Biol. diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 4cb751bbf..58f73f86a 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.166 +* `rugnux --model` places the model against the data as a rigid body before scoring it, keeping the move only if R-free improves, and writes sigma_A-weighted 2mFo-DFc and mFo-DFc maps in place of the unweighted 2Fo-Fc and Fo-Fc. * `rugnux` reads PILATUS miniCBF sweeps natively, and opens masters written at other facilities, including Eiger 1.x and third-party NXmx. * `rugnux` measures the beam centre on every run, and indexes with it when the file's value indexes nothing. * `rugnux` places a detector swung out on a 2theta arm where the file says it stands. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index e68fbfaef..86bc57eb2 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -17,7 +17,7 @@ This document describes the crystallographic algorithms implemented in Jungfrauj 11. space-group determination from the merged intensities (Laue group, screw axes, centering) and the twinning check, 12. auxiliary statistics (Wilson plot, ⟨I/σ(I)⟩, CC1/2, CCref), 13. amplitude estimation (French–Wilson) and R-free test-set flagging, -14. optional model-based validation: R-free against a supplied model, 2Fo−Fc / Fo−Fc electron-density maps, and an anomalous difference map with the strongest anomalous sites named. +14. optional model-based validation: rigid-body placement of a supplied model, R-free against it, sigma_A-weighted 2mFo−DFc / mFo−DFc electron-density maps, and an anomalous difference map with the strongest anomalous sites named. The reference is split into four parts, in pipeline order; the section numbers run continuously @@ -53,6 +53,7 @@ The methods draw on, and in places reimplement, solutions from: - M. Wojdyr, "GEMMI: A library for structural biology", *J. Open Source Softw.* **7** (2022), 4200 (model / structure-factor / map machinery used in §14). - 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. J. Read, "Improved Fourier coefficients for maps using phases from partial structures with errors", *Acta Cryst.* **A42** (1986), 140-149 (the sigma_A formalism and the m, D weighting of the map coefficients of §14.4). - Z. Otwinowski & W. Minor, "Processing of X-ray diffraction data collected in oscillation mode", *Methods Enzymol.* **276** (1997), 307-326 (reweighted, de-biased profile-fit variances). - G. Winter et al., "DIALS: implementation and evaluation of a new integration package", *Acta Cryst.* **D74** (2018), 85-97, and J. Beilsten-Edmands et al., *Acta Cryst.* **D76** (2020), 385-399 (CC1/2 resolution cutoff, merge outlier rejection, scaling error model). - P. Evans, "Scaling and assessment of data quality", *Acta Cryst.* **D62** (2006), 72-82, and P. R. Evans, *Acta Cryst.* **D67** (2011), 282-292 (POINTLESS: operator-by-operator point-group scoring, and the axial-zone screw-absence test). diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index a146c1863..bab5c6649 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -89,7 +89,7 @@ Everything lands in `_report.txt` section 9 (`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 **initial** 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 only re-fractionalized into the data unit cell (a rigid cell adjustment, so a deposited model with a slightly different cell still lines up), 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**; the model is re-fractionalized into the data unit cell and then **placed as one rigid body** (§14.7), 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 @@ -117,7 +117,11 @@ with R-free the same sum restricted to the free set. Note that the scaling of § ### 14.4 Electron-density maps -Two maps are formed with the model phases $\varphi_\mathrm{model}$: a $2F_o-F_c$ map, coefficients $(2|F_o|-|F_\mathrm{model}|)\,e^{i\varphi_\mathrm{model}}$, and an $F_o-F_c$ difference map, $(|F_o|-|F_\mathrm{model}|)\,e^{i\varphi_\mathrm{model}}$, each inverse-Fourier-transformed to a real-space CCP4 map (`_2fofc.ccp4`, `_fofc.ccp4`). A map-coefficient MTZ (`_maps.mtz`: `FP`, `FC`, `PHIC`, `FWT`/`PHWT`, `DELFWT`/`PHDELWT`, `FREE`) is written alongside so the maps can be reopened or rebuilt in Coot / PyMOL. These are unweighted difference coefficients (no $\sigma_A$ / figure-of-merit weighting), which is why they are described as *initial* maps. +Two maps are formed with the model phases $\varphi_\mathrm{model}$, both $\sigma_A$ weighted: a $2mF_o-DF_c$ map and an $mF_o-DF_c$ difference map, each inverse-Fourier-transformed to a real-space CCP4 map (`_2fofc.ccp4`, `_fofc.ccp4`). A map-coefficient MTZ (`_maps.mtz`: `FP`, `FC`, `PHIC`, `FWT`/`PHWT`, `DELFWT`/`PHDELWT`, `FOM`, `FREE`) is written alongside so the maps can be reopened or rebuilt in Coot / PyMOL. + +$\sigma_A$ is estimated by maximum likelihood **per resolution shell**, on the **free reflections only** — on the working set the model has been fitted to the data, so $\sigma_A$ would come out too high and the weighting would understate exactly the model error the map is meant to reveal. The shells are cut by equal reflection count, and it is the *number of shells* that is chosen from the size of the free set (about 50 free reflections to a shell, at most 20 shells), so no shell is thin by construction; a shell that still ends up with fewer than 10 free reflections takes the estimate made over the whole free set instead. Acentric and centric reflections enter with their own likelihoods (Rice and Woolfson respectively), and the epsilon factor is divided out before normalising both amplitudes to $\langle|E|^2\rangle = 1$ within the shell. The figure of merit is then $m = I_1(X)/I_0(X)$ with $X = 2\sigma_A|E_o||E_c|/(1-\sigma_A^2)$ for an acentric reflection and $m = \tanh(X/2)$ for a centric one, and $D = \sigma_A\sqrt{\Sigma_o/\Sigma_c}$ carries $F_c$ onto the observed amplitudes' scale. The $2F_o-F_c$ combination becomes $2mF_o - DF_c$ for acentric reflections and $mF_o$ for centric ones — a centric reflection's phase is either exactly right or 180° wrong, never in between — while the difference coefficient is $mF_o - DF_c$ throughout. + +$m$ and $D$ are estimated on each dataset separately, so two datasets of one crystal form get slightly different weights and their maps are to that extent no longer scaled identically — the same property §14.2 deliberately protects by refusing a free-form per-shell rescale. The two are not the same thing: the weighting never rescales $F_o$ and never touches the R-factors of §14.3, and the difference between two datasets' $\sigma_A$ curves is the difference in how well the model explains each of them, which is what a screening campaign is looking for. A PanDDA-style analysis consumes $2mF_o-DF_c$ maps and normalises each dataset's map against the ensemble before comparing them. The per-reflection `FOM` is written to the MTZ so the weighting can be read off and undone. ### 14.5 Aligning the data to the model: enantiomorph and indexing ambiguity @@ -141,3 +145,15 @@ Rather than searching the map for blobs and leaving a list of coordinates, the m Because the map is built on the model's phases and the data's own indexing, it is also the only test of whether the two agree about the **hand** (§14.5). A model and a dataset in opposite hands turn every anomalous peak into a trough, so a map whose deepest hole at an atom is both deeper than $5\sigma$ and deeper than its highest peak says so, and the run reports it as a warning naming that atom. It is not repaired by reindexing: that would make the two agree by construction and destroy the evidence for which of the model and the data is in the wrong hand. Note that R-free cannot see this at all — a mirrored model gives R-free to four decimal places unchanged, and an exactly inverted anomalous map. The list is always ten entries long, so it is their height that carries the information: on a sulfur-SAD dataset the sulfurs fill the top of the list and are followed by a clear drop to the couple of sigma that is the map's noise, while a dataset with no anomalous signal has no such separation and lists ten unrelated atoms at noise level. A scatterer the **model does not contain** — a bound ion, a soaked heavy atom — is by construction invisible in the list, and is what the map file is for. + +### 14.7 Rigid-body placement of the model + +Re-fractionalizing a model into the data cell puts it in the right box but not necessarily in the right place: a non-isomorphous cell squeezes the box without moving the body inside it, and the body's own position in the cell differs from crystal to crystal. Six parameters recover that — an angle-axis rotation about the model's own centroid, then a translation. Rotating about the centroid rather than the cell origin is what keeps the rotation from moving the body bodily, so the two triplets are close to independent. Parameters are carried as six lengths in ångström (the rotation vector multiplied by the model's r.m.s. radius), so a unit of each moves a typical atom by the same amount. + +It is **one** rigid body. A fragment-screening model arrives already solved and isomorphous, and what is being recovered is the crystal's movement, not the molecule's; splitting it into domains, or giving a bound ligand six parameters of its own, would refine against evidence these data do not separately carry — and the ligand is what the difference map is there to show, not to model away. + +The refinement walks a coarse-to-fine ladder, 6 Å → 4.5 Å → 3.5 Å, each zone starting from the previous one's answer. It stops at 3.5 Å because that is where rigid-body refinement is conventionally run and because the movement being recovered is a few tenths of an ångström, a tenth of that resolution; a finer zone costs $(1/d)^3$ in grid points and reflections for a placement it cannot meaningfully sharpen. Each evaluation recomputes $F_\mathrm{calc}$ and the solvent mask for the moved model **and re-fits the scale of §14.2** — otherwise the target would measure the scale as much as the placement, and the body would translate to repair a scale error instead of moving to where the density is. The minimiser is Levenberg–Marquardt on the amplitude residuals with a forward-difference Jacobian (seven evaluations per Jacobian; GEMMI's structure-factor path carries no analytic $\partial F/\partial p$, and at six parameters it does not need one), the step being a fixed fraction of the zone's resolution. + +The refinement sees **only the working reflections**. The step is then committed only if it lowers **R-free**, computed at full resolution on the free set it never saw; otherwise the model is put back exactly where it was read and the maps are the ones it would have given. The whole addition costs about two seconds. Note that the origin is a gauge in some space groups — free in all three directions in $P1$, and along the unique axis in a polar group — so those components of the translation are undetermined; nothing is done about that beyond the Levenberg–Marquardt damping and the R-free gate, which between them make an undetermined direction harmless rather than unstable. + +The rotation and translation actually taken are reported (`RIGID_BODY_ROTATION_DEG`, `RIGID_BODY_SHIFT_A`), together with the R-free before it (`R_FREE_BEFORE_RIGID_BODY`), so what the placement bought is visible. The **placed coordinates are not written out** — the deliverable is the map — so a model overlaid on these maps is the file as it was supplied, out by the reported shift. diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index fe8418a5a..458254fe9 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -10,6 +10,10 @@ ADD_LIBRARY(Rugnux STATIC RugnuxCommandLine.h ModelValidation.cpp ModelValidation.h + RigidBodyRefine.cpp + RigidBodyRefine.h + SigmaA.cpp + SigmaA.h ResultReport.cpp ResultReport.h SpotWidth.cpp diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index e00576174..361747a77 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -24,6 +24,8 @@ #include "../common/JFJochMath.h" // PI (M_PI is not standard, and MSVC does not define it) #include "../common/Logger.h" #include "../image_analysis/scale_merge/ReindexAmbiguity.h" // ReindexReflections +#include "RigidBodyRefine.h" +#include "SigmaA.h" namespace { @@ -159,26 +161,29 @@ ModelValidationResult ValidateAgainstModel(const std::vector & logger.Info("Model validation: {} atoms, cell a={:.2f} b={:.2f} c={:.2f}, sg {}, to {:.2f} A", gemmi::count_atom_sites(st.models[0]), ucell.a, ucell.b, ucell.c, sg->hm, d_min); - // --- Fcalc (atomic) via electron density on a grid + FFT --- - gemmi::DensityCalculator dc; - dc.d_min = d_min; - dc.rate = 1.5; - dc.set_grid_cell_and_spacegroup(st); - dc.set_refmac_compatible_blur(st.models[0]); - dc.put_model_density_on_grid(st.models[0]); - gemmi::AsuData> fcalc = - gemmi::transform_map_to_f_phi(dc.grid, true).prepare_asu_data(dc.d_min, dc.blur, false, false, false); + // --- Fcalc (atomic) via electron density on a grid + FFT, plus a flat bulk-solvent mask -> Fmask. + // A lambda because the rigid-body step below moves the model and then needs both again. --- + gemmi::AsuData> fcalc, fmask; + auto compute_model_factors = [&]() { + gemmi::DensityCalculator dc; + dc.d_min = d_min; + dc.rate = 1.5; + dc.set_grid_cell_and_spacegroup(st); + dc.set_refmac_compatible_blur(st.models[0]); + dc.put_model_density_on_grid(st.models[0]); + fcalc = gemmi::transform_map_to_f_phi(dc.grid, true) + .prepare_asu_data(dc.d_min, dc.blur, false, false, false); - // --- flat bulk-solvent mask -> Fmask --- - // Refmac radii give a slightly lower R than the Cctbx set on our test cases, at no cost. - gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac); - gemmi::Grid mask_grid; - mask_grid.unit_cell = dc.grid.unit_cell; - mask_grid.spacegroup = dc.grid.spacegroup; - mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up); - masker.put_mask_on_grid(mask_grid, st.models[0]); - gemmi::AsuData> fmask = - gemmi::transform_map_to_f_phi(mask_grid, true).prepare_asu_data(dc.d_min, 0); + // Refmac radii give a slightly lower R than the Cctbx set on our test cases, at no cost. + gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac); + gemmi::Grid mask_grid; + mask_grid.unit_cell = dc.grid.unit_cell; + mask_grid.spacegroup = dc.grid.spacegroup; + mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up); + masker.put_mask_on_grid(mask_grid, st.models[0]); + fmask = gemmi::transform_map_to_f_phi(mask_grid, true).prepare_asu_data(dc.d_min, 0); + }; + compute_model_factors(); gemmi::GroupOps gops = sg->operations(); gemmi::ReciprocalAsu asu(sg); @@ -187,7 +192,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // Factored into a lambda so we can probe indexing (merohedral) ambiguities: run the same scale + // R computation on each reindexing of the observed reflections and keep the lowest-R-free one. struct Fit { - gemmi::AsuData> fmodel, map2fofc, mapfofc; + gemmi::AsuData> fmodel; + gemmi::AsuData> fobs; // what it was fitted to std::unordered_map> obs_by_hkl; // hkl -> (Fobs, is_free) double r_work = 1, r_free = 1, k_sol = 0, b_sol = 0, k_overall = 0; int n_w = 0, n_f = 0; @@ -200,7 +206,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // Observed amplitudes are the French-Wilson |F| already computed at the end of the merge // (MergedReflection.F), so the model R-free / maps use exactly the same amplitudes as the // written reflection file. - gemmi::AsuData> fobs; + gemmi::AsuData> &fobs = out.fobs; fobs.unit_cell_ = ucell; fobs.spacegroup_ = sg; for (const MergedReflection &r : obs_in) { @@ -230,22 +236,19 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // it reshapes each map's radial amplitude profile differently, so a batch of maps would no // longer be directly comparable. For a fragment-screening / PanDDA campaign, comparable maps // across datasets matter more than the last bit of per-dataset R, so it is deliberately omitted. + // (The sigma_A weighting further down is a different thing and does not reopen this: it never + // rescales Fobs, and it leaves the R-factors below untouched. It does weight the map + // coefficients per shell and per dataset - see the note where it is computed.) - // --- R-work / R-free and map coefficients (2Fo-Fc and Fo-Fc, model phases) --- - out.map2fofc.unit_cell_ = ucell; out.map2fofc.spacegroup_ = sg; - out.mapfofc.unit_cell_ = ucell; out.mapfofc.spacegroup_ = sg; + // --- R-work / R-free --- double num_w = 0, den_w = 0, num_f = 0, den_f = 0; for (const auto &hv : out.fmodel.v) { auto it = out.obs_by_hkl.find(hkl_key(hv.hkl)); if (it == out.obs_by_hkl.end()) continue; double Fo = it->second.first; double Fc = std::abs(hv.value); - double phi = std::arg(hv.value); if (it->second.second) { num_f += std::fabs(Fo - Fc); den_f += Fo; ++out.n_f; } else { num_w += std::fabs(Fo - Fc); den_w += Fo; ++out.n_w; } - std::complex ph = std::polar(1.0f, static_cast(phi)); - out.map2fofc.v.push_back({hv.hkl, static_cast(2 * Fo - Fc) * ph}); - out.mapfofc.v.push_back({hv.hkl, static_cast(Fo - Fc) * ph}); } out.r_work = den_w > 0 ? num_w / den_w : 1; out.r_free = den_f > 0 ? num_f / den_f : 1; @@ -261,13 +264,21 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // enantiomorph/screw ambiguity is never probed by R-free: |Fcalc| is the same for both hands, so // it cannot distinguish them - that is taken from the model hand above. Fit best = fit_model(obs); + std::vector reindexed; // kept only where a reindexing wins + const std::vector *obs_best = &obs; if (probe_indexing_ambiguity) { const auto reindex_ops = ReindexAmbiguityOperators(cell, *sg); std::vector candidate_r_free{best.r_free}; // identity first, then the twin laws for (const auto &op : reindex_ops) { - Fit cand = fit_model(ReindexReflections(obs, op)); + std::vector cand_obs = ReindexReflections(obs, op); + Fit cand = fit_model(cand_obs); candidate_r_free.push_back(cand.r_free); - if (cand.r_free < best.r_free) { best = std::move(cand); result.indexing_op = op; } + if (cand.r_free < best.r_free) { + best = std::move(cand); + result.indexing_op = op; + reindexed = std::move(cand_obs); + obs_best = &reindexed; + } } if (!reindex_ops.empty()) { // The runner-up as well as the winner: the margin between them is what says whether the @@ -281,9 +292,46 @@ ModelValidationResult ValidateAgainstModel(const std::vector & } } + // --- rigid-body placement of the model in the data cell --- + // Re-fractionalizing the model above puts it in the right box but not in the right place: a + // non-isomorphous cell squeezes the box without moving the body inside it, and the body's own + // position differs from crystal to crystal anyway. Six parameters recover that - a rotation about + // the model's centroid and a translation - which is all a fragment-screening model needs, since it + // arrives already solved. The step is committed only if the FREE reflections, which the refinement + // never saw, agree that it helped; on data the model cannot be placed against, the model stays + // exactly where it was read. + { + gemmi::AsuData> fobs_work; + fobs_work.unit_cell_ = ucell; + fobs_work.spacegroup_ = sg; + for (const auto &hv : best.fobs.v) { + const auto it = best.obs_by_hkl.find(hkl_key(hv.hkl)); + if (it != best.obs_by_hkl.end() && !it->second.second) + fobs_work.v.push_back(hv); + } + const std::vector before = ModelPositions(st.models[0]); + const RigidBodyRefineResult rb = + RefineRigidBody(st.models[0], ucell, *sg, fobs_work, d_min, logger); + if (rb.converged) { + compute_model_factors(); + Fit moved = fit_model(*obs_best); + const bool commit = moved.r_free < best.r_free; + logger.Info("Model validation: rigid body held-out R-free {:.4f} -> {:.4f} => {}", + best.r_free, moved.r_free, commit ? "committed" : "rejected, model put back"); + if (commit) { + result.rigid_body_applied = true; + result.rigid_body_angle_deg = rb.angle_deg; + result.rigid_body_shift_A = rb.shift_A; + result.r_free_before_rigid_body = best.r_free; + best = std::move(moved); + } else { + SetModelPositions(st.models[0], before); + compute_model_factors(); + } + } + } + gemmi::AsuData> &fmodel = best.fmodel; - gemmi::AsuData> &map2fofc = best.map2fofc; - gemmi::AsuData> &mapfofc = best.mapfofc; std::unordered_map> &obs_by_hkl = best.obs_by_hkl; result.r_work = best.r_work; @@ -294,7 +342,53 @@ ModelValidationResult ValidateAgainstModel(const std::vector & result.b_sol = best.b_sol; result.k_overall = best.k_overall; - // --- write the maps and score the 2Fo-Fc map at atom centres (a real map peaks there) --- + // --- sigma_A weighting: the maps are 2mFo-DFc and mFo-DFc, not 2Fo-Fc and Fo-Fc --- + // m and D come from a maximum-likelihood sigma_A per resolution shell, so a shell the model + // describes badly is damped rather than carried into the map at full weight, and the difference + // map is correspondingly less biased towards the model that made its phases. + // + // m and D are estimated on THIS dataset, so two datasets of one crystal form get slightly + // different weights, and to that extent their maps are no longer scaled identically - the same + // property the scaling above deliberately protects. It is kept anyway: the difference between two + // datasets' sigma_A curves is the difference in how well the model explains each of them, which is + // real and is what a screening campaign is looking for, and a PanDDA-style analysis consumes + // 2mFo-DFc maps and normalizes each dataset's map against the ensemble before comparing them. The + // per-reflection FOM is written to the MTZ so the weighting can be undone. + // Following Read (1986) Acta Cryst. A42, 140-149 + struct MapTerm { gemmi::Miller hkl; double fo, fc, phi; bool free, centric; }; + std::vector terms; + std::vector sa_input; + for (const auto &hv : fmodel.v) { + const auto it = obs_by_hkl.find(hkl_key(hv.hkl)); + if (it == obs_by_hkl.end()) continue; + const double Fo = it->second.first; + const double Fc = std::abs(hv.value); + const bool centric = gops.is_reflection_centric(hv.hkl); + terms.push_back({hv.hkl, Fo, Fc, std::arg(hv.value), it->second.second, centric}); + sa_input.push_back({Fo, Fc, ucell.calculate_1_d2(hv.hkl), gops.epsilon_factor(hv.hkl), + centric, it->second.second}); + } + const SigmaAResult sigma_a = EstimateSigmaA(sa_input, ucell); + result.mean_fom = sigma_a.mean_fom; + result.sigma_a_shells = sigma_a.shells; + + gemmi::AsuData> map2fofc, mapfofc; + map2fofc.unit_cell_ = ucell; map2fofc.spacegroup_ = sg; + mapfofc.unit_cell_ = ucell; mapfofc.spacegroup_ = sg; + std::vector fwt(terms.size()), delfwt(terms.size()); + for (size_t i = 0; i < terms.size(); i++) { + const double m = sigma_a.weight[i].m, D = sigma_a.weight[i].d; + // A centric reflection's phase is either exactly right or 180 degrees wrong, never in + // between, so its bias-free coefficient is mFo and not 2mFo - DFc. + fwt[i] = static_cast(terms[i].centric ? m * terms[i].fo + : 2 * m * terms[i].fo - D * terms[i].fc); + delfwt[i] = static_cast(m * terms[i].fo - D * terms[i].fc); + const std::complex ph = std::polar(1.0f, static_cast(terms[i].phi)); + map2fofc.v.push_back({terms[i].hkl, fwt[i] * ph}); + mapfofc.v.push_back({terms[i].hkl, delfwt[i] * ph}); + } + + // --- write the maps and score the 2mFo-DFc map at atom centres (a real map peaks there) --- const gemmi::Grid grid2fofc = map_from_coefficients(map2fofc); const double rms2 = write_ccp4(grid2fofc, output_prefix + "_2fofc.ccp4"); write_ccp4(map_from_coefficients(mapfofc), output_prefix + "_fofc.ccp4"); @@ -400,25 +494,23 @@ ModelValidationResult ValidateAgainstModel(const std::vector & mtz.add_column("PHWT", 'P', -1, -1, false); mtz.add_column("DELFWT", 'F', -1, -1, false); mtz.add_column("PHDELWT", 'P', -1, -1, false); + // The figure of merit the coefficients carry, so the weighting can be read off - and undone - + // from the file rather than having to be taken on trust. + mtz.add_column("FOM", 'W', -1, -1, false); mtz.add_column("FREE", 'I', -1, -1, false); std::vector data; - int nref = 0; - for (const auto &hv : fmodel.v) { - auto it = obs_by_hkl.find(hkl_key(hv.hkl)); - if (it == obs_by_hkl.end()) continue; - double Fo = it->second.first; - double Fc = std::abs(hv.value); - double phi_deg = gemmi::phase_in_angles(hv.value); - data.insert(data.end(), {static_cast(hv.hkl[0]), static_cast(hv.hkl[1]), - static_cast(hv.hkl[2]), - static_cast(Fo), static_cast(Fc), - static_cast(phi_deg), - static_cast(2 * Fo - Fc), static_cast(phi_deg), - static_cast(Fo - Fc), static_cast(phi_deg), - it->second.second ? 0.0f : 1.0f}); - ++nref; + for (size_t i = 0; i < terms.size(); i++) { + const auto phi_deg = static_cast(terms[i].phi * 180.0 / PI); + data.insert(data.end(), {static_cast(terms[i].hkl[0]), static_cast(terms[i].hkl[1]), + static_cast(terms[i].hkl[2]), + static_cast(terms[i].fo), static_cast(terms[i].fc), + phi_deg, + fwt[i], phi_deg, + delfwt[i], phi_deg, + static_cast(sigma_a.weight[i].m), + terms[i].free ? 0.0f : 1.0f}); } - mtz.nreflections = nref; + mtz.nreflections = static_cast(terms.size()); mtz.data = std::move(data); mtz.write_to_file(output_prefix + "_maps.mtz"); } catch (const std::exception &e) { @@ -432,7 +524,11 @@ ModelValidationResult ValidateAgainstModel(const std::vector & result.r_work, result.n_work, result.r_free, result.n_free); logger.Info("Model validation: bulk solvent k_sol={:.3f} b_sol={:.1f}, k_overall={:.3f}", result.k_sol, result.b_sol, result.k_overall); - logger.Info("Model validation: mean 2Fo-Fc density at atom centres = {:.2f} sigma", result.mean_atom_density_sigma); + logger.Info("Model validation: mean 2mFo-DFc density at atom centres = {:.2f} sigma", result.mean_atom_density_sigma); + logger.Info("Model validation: sigma_A weighting over {} resolution shell(s) (estimated on the {} " + "free reflections): sigma_A {:.2f} at low resolution, {:.2f} at high, mean FOM {:.3f}", + sigma_a.shells, sigma_a.free_reflections, sigma_a.sigma_a_lowest_shell, + sigma_a.sigma_a_highest_shell, sigma_a.mean_fom); if (!result.anomalous_sites.empty()) { std::string sites; for (const auto &s : result.anomalous_sites) diff --git a/rugnux/ModelValidation.h b/rugnux/ModelValidation.h index 24f5073ec..fb76b0391 100644 --- a/rugnux/ModelValidation.h +++ b/rugnux/ModelValidation.h @@ -29,9 +29,21 @@ struct ModelValidationResult { double k_overall = 0.0; int n_work = 0; int n_free = 0; - double mean_atom_density_sigma = 0.0; // mean 2Fo-Fc value at atom centres, in sigma + double mean_atom_density_sigma = 0.0; // mean 2mFo-DFc value at atom centres, in sigma std::string maps_prefix; // where the .ccp4 / _maps.mtz were written + // The rigid-body placement of the model against these data, where it was committed. The step is + // taken only when it lowers R-free, so an unmoved model means the data said no, not that nothing + // was tried; r_free_before_rigid_body is what it would have been without it. + bool rigid_body_applied = false; + double rigid_body_angle_deg = 0.0; + double rigid_body_shift_A = 0.0; + double r_free_before_rigid_body = 0.0; + + // The sigma_A weighting the map coefficients carry. + double mean_fom = 1.0; + int sigma_a_shells = 0; + // The anomalous difference map read at the model's own atoms: the strongest sites, highest // first. Empty when the merge kept no Bijvoet split, and so had nothing to make the map from. struct AnomalousSite { @@ -65,9 +77,12 @@ struct ModelValidationResult { // amplitudes, then report R-work / R-free and write 2Fo-Fc and Fo-Fc electron-density maps // (CCP4) plus an MTZ of map coefficients next to output_prefix. // -// No refinement of the structure itself is done. The model is only re-fractionalized into the -// data unit cell (a cheap rigid cell adjustment) so a deposited model with a slightly different -// cell still lines up with the processed data. Returns ok=false (and logs) on any failure. +// The structure itself is not refined. The model is re-fractionalized into the data unit cell and +// then placed as ONE RIGID BODY - an angle-axis rotation about its centroid plus a translation, +// refined over a coarse-to-fine resolution ladder and kept only where it lowers R-free - so a +// deposited model from a non-isomorphous crystal still sits where the density is. The map +// coefficients are sigma_A weighted (2mFo-DFc and mFo-DFc). Returns ok=false (and logs) on any +// failure. // data_space_group is the space group the data was merged in. If it is the enantiomorph of // the model's space group (e.g. data P4(1)2(1)2, model P4(3)2(1)2 - indistinguishable from merged // intensities), the model's group is adopted as a label, without touching the reflections. diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 66dd9e8b3..e0c36b45c 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -584,6 +584,17 @@ std::string RenderResultReport(const std::string &output_prefix, Key(os, "BULK_SOLVENT_K_SOL", fmt::format("{:.3f}", mv.k_sol)); Key(os, "BULK_SOLVENT_B_SOL", fmt::format("{:.1f}", mv.b_sol)); Key(os, "SCALE_OVERALL", fmt::format("{:.4f}", mv.k_overall)); + // The model is placed against the data as one rigid body, and only where the free + // reflections say it helped - so R_FREE_BEFORE_RIGID_BODY says what the placement bought. + Key(os, "RIGID_BODY", mv.rigid_body_applied ? "APPLIED" : "NOT_APPLIED"); + if (mv.rigid_body_applied) { + Key(os, "RIGID_BODY_ROTATION_DEG", fmt::format("{:.3f}", mv.rigid_body_angle_deg)); + Key(os, "RIGID_BODY_SHIFT_A", fmt::format("{:.3f}", mv.rigid_body_shift_A)); + Key(os, "R_FREE_BEFORE_RIGID_BODY", fmt::format("{:.4f}", mv.r_free_before_rigid_body)); + } + Key(os, "MAP_COEFFICIENTS", "2mFo-DFc / mFo-DFc"); + Key(os, "MAP_SIGMA_A_SHELLS", mv.sigma_a_shells); + Key(os, "MAP_MEAN_FOM", fmt::format("{:.3f}", mv.mean_fom)); Key(os, "MEAN_ATOM_DENSITY_SIGMA", fmt::format("{:.2f}", mv.mean_atom_density_sigma)); // The anomalous scatterers the data themselves found, named by the model's atoms. if (!mv.anomalous_sites.empty()) { diff --git a/rugnux/RigidBodyRefine.cpp b/rugnux/RigidBodyRefine.cpp new file mode 100644 index 000000000..29f4fcc99 --- /dev/null +++ b/rugnux/RigidBodyRefine.cpp @@ -0,0 +1,295 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "RigidBodyRefine.h" + +#include +#include +#include +#include +#include + +#include +#include + +#include "gemmi/dencalc.hpp" // DensityCalculator +#include "gemmi/fourier.hpp" // transform_map_to_f_phi +#include "gemmi/it92.hpp" // IT92 x-ray form factors +#include "gemmi/scaling.hpp" // Scaling (bulk solvent + anisotropic B) +#include "gemmi/solmask.hpp" // SolventMasker + +#include "../common/JFJochMath.h" // PI +#include "../common/Logger.h" + +namespace { + +using Table = gemmi::IT92; + +// The ladder the placement is walked down. It starts coarse because the model arrives already placed +// but out by a cell's worth of non-isomorphism: at 6 A a few hundred reflections see the body as a +// blob and the target has one broad minimum, and each finer zone starts from the previous one's +// answer. It stops at 3.5 A, which is where rigid-body refinement is conventionally run (it is +// REFMAC's own default through dimple) - the movement being recovered is a few tenths of an +// angstrom, a tenth of that resolution, so it is well determined there, while a finer zone costs +// (1/d)^3 in grid points and reflections for a placement it cannot meaningfully sharpen. +constexpr double LADDER[] = {6.0, 4.5, 3.5}; + +// The step of the forward-difference Jacobian, as a fraction of the zone's resolution - so it is +// 0.06 A of atom displacement at 6 A and 0.035 A at 3.5 A. A step fixed in angstroms instead is far +// too small for the coarse zones, where a structure factor barely notices it and the derivative is +// swallowed by the jitter of the scale re-fit: measured, a fixed 0.02 A left the 6 A zone at 0.35 +// degrees where this rule takes it to 2.79, which is most of the way to the answer. +constexpr double JACOBIAN_STEP_FRACTION = 0.01; + +// Parameters are carried as six lengths in angstroms - the first three are the angle-axis rotation +// vector multiplied by the model's rms radius, so a unit of each of the six moves a typical atom by +// the same amount. That makes the Jacobian step and the trust region isotropic in something +// physical, rather than mixing radians with angstroms. +struct Placement { + gemmi::Position centre; // the model centroid: rotating about it decorrelates R from t + double rms_radius = 1.0; // rms distance of the atoms from the centroid + + void Apply(const double q[6], const std::vector &base, gemmi::Model &model) const { + const double aa[3] = {q[0] / rms_radius, q[1] / rms_radius, q[2] / rms_radius}; + size_t i = 0; + for (gemmi::Chain &ch : model.chains) + for (gemmi::Residue &r : ch.residues) + for (gemmi::Atom &a : r.atoms) { + const double p[3] = {base[i].x - centre.x, base[i].y - centre.y, base[i].z - centre.z}; + double rp[3]; + ceres::AngleAxisRotatePoint(aa, p, rp); + a.pos = gemmi::Position(rp[0] + centre.x + q[3], + rp[1] + centre.y + q[4], + rp[2] + centre.z + q[5]); + ++i; + } + } +}; + +// One target evaluation: place the model, recompute Fcalc and the bulk-solvent mask to the zone's +// resolution, re-fit the scale, and hand back the amplitude residuals. +class Evaluator { +public: + Evaluator(gemmi::Model &model, const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, + const std::vector &base, const Placement &placement) + : model_(model), cell_(cell), sg_(sg), base_(base), placement_(placement) {} + + // The zone's observations, and the scale the residuals are expressed in. + void SetZone(const gemmi::AsuData> &fobs, double d_min) { + fobs_ = fobs; + d_min_ = d_min; + double sum = 0; + for (const auto &hv : fobs_.v) + sum += hv.value.value; + f_mean_ = fobs_.v.empty() ? 1.0 : sum / static_cast(fobs_.v.size()); + } + + size_t NumObservations() const { return fobs_.v.size(); } + double JacobianStep() const { return JACOBIAN_STEP_FRACTION * d_min_; } + int evaluations = 0; + + bool Residuals(const double q[6], double *residuals) { + ++evaluations; + placement_.Apply(q, base_, model_); + + gemmi::DensityCalculator dc; + dc.d_min = d_min_; + dc.rate = 1.5; + dc.grid.unit_cell = cell_; + dc.grid.spacegroup = &sg_; + dc.set_refmac_compatible_blur(model_); + dc.put_model_density_on_grid(model_); + gemmi::AsuData> fcalc = + gemmi::transform_map_to_f_phi(dc.grid, true).prepare_asu_data(dc.d_min, dc.blur, false, false, false); + + gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac); + gemmi::Grid mask_grid; + mask_grid.unit_cell = cell_; + mask_grid.spacegroup = &sg_; + mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up); + masker.put_mask_on_grid(mask_grid, model_); + gemmi::AsuData> fmask = + gemmi::transform_map_to_f_phi(mask_grid, true).prepare_asu_data(dc.d_min, 0); + if (fmask.size() != fcalc.size()) + return false; + + // Re-fitted at every evaluation: with the scale held at the starting placement's value the + // target would measure the scale as much as the placement, and the body would translate to + // repair a scale error instead of moving where the density is. + gemmi::Scaling scaling(cell_, &sg_); + scaling.use_solvent = true; + scaling.prepare_points(fcalc, fobs_, &fmask); + if (scaling.points.empty()) + return false; + scaling.fit_isotropic_b_approximately(); + scaling.fit_parameters(); + scaling.scale_data(fcalc, &fmask); + + // Both are sorted and in the same ASU, so one merge pass matches them. + auto c = fcalc.v.begin(); + for (size_t i = 0; i < fobs_.v.size(); ++i) { + const gemmi::Miller &h = fobs_.v[i].hkl; + while (c != fcalc.v.end() && c->hkl < h) + ++c; + residuals[i] = (c != fcalc.v.end() && c->hkl == h) + ? (fobs_.v[i].value.value - std::abs(c->value)) / f_mean_ + : 0.0; + } + return true; + } + +private: + gemmi::Model &model_; + const gemmi::UnitCell &cell_; + const gemmi::SpaceGroup &sg_; + const std::vector &base_; + Placement placement_; + gemmi::AsuData> fobs_; + double d_min_ = 0; + double f_mean_ = 1; +}; + +// Ceres' own numeric differentiation steps by |x| * relative_step_size, which is zero at the start of +// every zone (the placement begins at no shift), so the Jacobian is supplied here instead, by +// forward differences at a step chosen in the parameters' units. Analytic dF/dp would need +// derivatives GEMMI's structure-factor path does not have, and at six parameters it is not worth it: +// a Jacobian costs seven evaluations, and the evaluations at 6-3.5 A are cheap. +class RigidBodyCost : public ceres::CostFunction { +public: + explicit RigidBodyCost(Evaluator &ev) : ev_(ev) { + set_num_residuals(static_cast(ev.NumObservations())); + mutable_parameter_block_sizes()->push_back(6); + } + + bool Evaluate(double const *const *parameters, double *residuals, double **jacobians) const override { + const int n = num_residuals(); + if (!ev_.Residuals(parameters[0], residuals)) + return false; + if (jacobians != nullptr && jacobians[0] != nullptr) { + std::vector shifted(n); + for (int j = 0; j < 6; j++) { + double q[6]; + std::copy(parameters[0], parameters[0] + 6, q); + const double step = ev_.JacobianStep(); + q[j] += step; + if (!ev_.Residuals(q, shifted.data())) + return false; + for (int i = 0; i < n; i++) + jacobians[0][i * 6 + j] = (shifted[i] - residuals[i]) / step; + } + } + return true; + } + +private: + Evaluator &ev_; +}; + +} // namespace + +std::vector ModelPositions(const gemmi::Model &model) { + std::vector pos; + for (const gemmi::Chain &ch : model.chains) + for (const gemmi::Residue &r : ch.residues) + for (const gemmi::Atom &a : r.atoms) + pos.push_back(a.pos); + return pos; +} + +void SetModelPositions(gemmi::Model &model, const std::vector &pos) { + size_t i = 0; + for (gemmi::Chain &ch : model.chains) + for (gemmi::Residue &r : ch.residues) + for (gemmi::Atom &a : r.atoms) + a.pos = pos[i++]; +} + +// One rigid body, not groups: a fragment-screening model arrives already solved and isomorphous, and +// the movement to recover is the crystal's, not the molecule's. Splitting it into domains or giving a +// bound ligand its own six parameters would refine against evidence this data does not separately +// carry, and the ligand is what the difference map is meant to show rather than model away. +// +// Some of the translation can be a gauge rather than a quantity - the origin is free in all three +// directions in P1 and along the unique axis in a polar group, and |F| does not change when the whole +// content moves along it. Nothing is done about that beyond the Levenberg-Marquardt damping, which +// leaves an undetermined direction where it started, and the caller's R-free gate, which throws away +// a step that wandered along one. +RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, + const gemmi::UnitCell &cell, + const gemmi::SpaceGroup &sg, + const gemmi::AsuData> &fobs, + double d_min, + Logger &logger) { + const auto t0 = std::chrono::steady_clock::now(); + RigidBodyRefineResult result; + + const std::vector base = ModelPositions(model); + if (base.empty() || fobs.v.empty()) + return result; + + Placement placement; + for (const gemmi::Position &p : base) + placement.centre += p; + placement.centre *= 1.0 / static_cast(base.size()); + double r2 = 0; + for (const gemmi::Position &p : base) + r2 += placement.centre.dist_sq(p); + placement.rms_radius = std::sqrt(r2 / static_cast(base.size())); + if (!(placement.rms_radius > 0)) + return result; + + for (double zone : LADDER) + if (zone >= d_min) + result.zones.push_back(zone); + if (result.zones.empty()) + result.zones.push_back(d_min); + + Evaluator ev(model, cell, sg, base, placement); + double q[6] = {0, 0, 0, 0, 0, 0}; + bool any_zone_solved = false; + for (double zone : result.zones) { + gemmi::AsuData> zone_obs; + zone_obs.unit_cell_ = fobs.unit_cell_; + zone_obs.spacegroup_ = fobs.spacegroup_; + for (const auto &hv : fobs.v) + if (cell.calculate_d(hv.hkl) >= zone) + zone_obs.v.push_back(hv); + if (zone_obs.v.size() < 50) + continue; + ev.SetZone(zone_obs, zone); + + ceres::Problem problem; + problem.AddResidualBlock(new RigidBodyCost(ev), nullptr, q); + ceres::Solver::Options options; + options.linear_solver_type = ceres::DENSE_QR; + options.max_num_iterations = 15; + options.function_tolerance = 1e-4; + options.parameter_tolerance = 1e-4; + options.logging_type = ceres::LoggingType::SILENT; + ceres::Solver::Summary summary; + const int evaluations_before = ev.evaluations; + const auto zone_t0 = std::chrono::steady_clock::now(); + ceres::Solve(options, &problem, &summary); + any_zone_solved = any_zone_solved || summary.IsSolutionUsable(); + logger.Debug("Rigid body zone {:.1f} A: {} reflections, {} iterations, {} evaluations, {:.2f} s, " + "rotation {:.3f} deg, translation {:.3f} A", zone, zone_obs.v.size(), + summary.iterations.size() - 1, ev.evaluations - evaluations_before, + std::chrono::duration(std::chrono::steady_clock::now() - zone_t0).count(), + std::sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2]) / placement.rms_radius * 180.0 / PI, + std::sqrt(q[3]*q[3] + q[4]*q[4] + q[5]*q[5])); + } + + placement.Apply(q, base, model); // Ceres left the model at a Jacobian probe; put it at the answer + result.evaluations = ev.evaluations; + result.converged = any_zone_solved; + const double aa = std::sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2]) / placement.rms_radius; + result.angle_deg = aa * 180.0 / PI; + result.shift_A = std::sqrt(q[3] * q[3] + q[4] * q[4] + q[5] * q[5]); + result.seconds = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); + + logger.Info("Model validation: rigid body over {} resolution zone(s) down to {:.1f} A, " + "{} evaluations in {:.2f} s: rotation {:.3f} deg, translation {:.3f} A", + result.zones.size(), result.zones.back(), result.evaluations, result.seconds, + result.angle_deg, result.shift_A); + return result; +} diff --git a/rugnux/RigidBodyRefine.h b/rugnux/RigidBodyRefine.h new file mode 100644 index 000000000..3bd51283f --- /dev/null +++ b/rugnux/RigidBodyRefine.h @@ -0,0 +1,42 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "gemmi/asudata.hpp" +#include "gemmi/model.hpp" +#include "gemmi/symmetry.hpp" +#include "gemmi/unitcell.hpp" + +class Logger; + +struct RigidBodyRefineResult { + bool converged = false; + double angle_deg = 0.0; // magnitude of the rotation about the model centroid + double shift_A = 0.0; // magnitude of the translation + std::vector zones; // the resolution ladder actually walked, coarsest first + int evaluations = 0; // structure-factor evaluations (each one re-fits the scale) + double seconds = 0.0; +}; + +// Every atom position of `model`, in the order the model iterates them. +std::vector ModelPositions(const gemmi::Model &model); +void SetModelPositions(gemmi::Model &model, const std::vector &pos); + +// Refine the placement of `model` as one rigid body against the observed amplitudes: an angle-axis +// rotation about the model's own centroid followed by a translation, six parameters, over a +// coarse-to-fine resolution ladder. Each evaluation recomputes Fcalc and the bulk-solvent mask for +// the moved model and re-fits the scale (k_overall, anisotropic B, k_sol, b_sol), so the target +// measures the placement and not the scale. +// +// `fobs` should be the working set only - the free reflections are what the caller decides on. The +// model is left MOVED (whether or not the refinement helped): the caller scores it and puts it back +// with SetModelPositions() if it did not. +RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, + const gemmi::UnitCell &cell, + const gemmi::SpaceGroup &sg, + const gemmi::AsuData> &fobs, + double d_min, + Logger &logger); diff --git a/rugnux/SigmaA.cpp b/rugnux/SigmaA.cpp new file mode 100644 index 000000000..fcf00a9d9 --- /dev/null +++ b/rugnux/SigmaA.cpp @@ -0,0 +1,151 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "SigmaA.h" + +#include +#include +#include + +#include "gemmi/bessel.hpp" // bessel_i1_over_i0, log_bessel_i0 +#include "gemmi/binner.hpp" // Binner +#include "gemmi/math.hpp" // log_cosh + +namespace { + +// How many free reflections a shell is given. sigma_A is one number per shell, so a few tens of +// reflections already pin it; sizing the SHELLS from the free count (rather than fixing the number of +// shells and then patching the thin ones) is what keeps every shell usable on a small dataset, where +// the free set is 5 % of a few thousand reflections. +constexpr int FREE_PER_SHELL = 50; +constexpr int MAX_SHELLS = 20; + +// Below this a shell's own maximum likelihood is noise, and the estimate over the whole free set is +// used for it instead. Equal-count shells are cut on all the reflections, so a shell's free count +// still fluctuates about FREE_PER_SHELL and the odd one comes out short. +constexpr int MIN_FREE_FOR_SHELL = 10; + +constexpr double SIGMA_A_MIN = 0.01; +constexpr double SIGMA_A_MAX = 0.99; + +struct Normalized { + double eo = 0, ec = 0; + bool centric = false; +}; + +// Rice (acentric) and Woolfson (centric) log likelihoods of |Eo| given |Ec| and sigma_A, with the +// terms that do not depend on sigma_A dropped. +double LogLikelihood(const std::vector &e, const std::vector &idx, double sigma_a) { + const double u = 1 - sigma_a * sigma_a; + double ll = 0; + for (int i : idx) { + const double eo = e[i].eo, ec = e[i].ec; + const double quad = eo * eo + sigma_a * sigma_a * ec * ec; + if (e[i].centric) + ll += -0.5 * std::log(u) - quad / (2 * u) + gemmi::log_cosh(sigma_a * eo * ec / u); + else + ll += -std::log(u) - quad / u + gemmi::log_bessel_i0(2 * sigma_a * eo * ec / u); + } + return ll; +} + +// The likelihood is unimodal in sigma_A, so a golden-section search finds its maximum without +// derivatives and without a starting guess. +double MaximizeSigmaA(const std::vector &e, const std::vector &idx) { + constexpr double GOLDEN = 0.6180339887498949; + double lo = SIGMA_A_MIN, hi = SIGMA_A_MAX; + double x1 = hi - GOLDEN * (hi - lo), x2 = lo + GOLDEN * (hi - lo); + double f1 = LogLikelihood(e, idx, x1), f2 = LogLikelihood(e, idx, x2); + for (int it = 0; it < 60; it++) { + if (f1 > f2) { + hi = x2; x2 = x1; f2 = f1; + x1 = hi - GOLDEN * (hi - lo); f1 = LogLikelihood(e, idx, x1); + } else { + lo = x1; x1 = x2; f1 = f2; + x2 = lo + GOLDEN * (hi - lo); f2 = LogLikelihood(e, idx, x2); + } + } + return 0.5 * (lo + hi); +} + +} // namespace + +SigmaAResult EstimateSigmaA(const std::vector &refl, const gemmi::UnitCell &cell) { + SigmaAResult result; + result.weight.assign(refl.size(), SigmaAWeight{}); + if (refl.empty()) + return result; + + for (const SigmaAReflection &r : refl) + if (r.free) + ++result.free_reflections; + if (result.free_reflections == 0) + return result; // nothing to estimate on: leave the coefficients unweighted + + const int nbins = std::clamp(result.free_reflections / FREE_PER_SHELL, 1, MAX_SHELLS); + std::vector inv_d2; + inv_d2.reserve(refl.size()); + for (const SigmaAReflection &r : refl) + inv_d2.push_back(r.inv_d2); + gemmi::Binner binner; + binner.setup_from_1_d2(nbins, gemmi::Binner::Method::EqualCount, std::move(inv_d2), &cell); + std::vector bin(refl.size()); + for (size_t i = 0; i < refl.size(); i++) + bin[i] = binner.get_bin_from_1_d2(refl[i].inv_d2); + result.shells = nbins; + + // Normalize both amplitudes to <|E|^2> = 1 within the shell, which is the scale the likelihood + // above is written in. The epsilon factor is the symmetry enhancement of a reflection's expected + // intensity, and dividing it out is what makes the reflections of one shell comparable. + std::vector sum_o(nbins, 0), sum_c(nbins, 0); + std::vector count(nbins, 0); + for (size_t i = 0; i < refl.size(); i++) { + const double eps = std::max(1, refl[i].epsilon); + sum_o[bin[i]] += refl[i].f_obs * refl[i].f_obs / eps; + sum_c[bin[i]] += refl[i].f_calc * refl[i].f_calc / eps; + ++count[bin[i]]; + } + std::vector sigma_o(nbins), sigma_c(nbins); + for (int b = 0; b < nbins; b++) { + sigma_o[b] = count[b] > 0 ? sum_o[b] / count[b] : 0; + sigma_c[b] = count[b] > 0 ? sum_c[b] / count[b] : 0; + } + + std::vector e(refl.size()); + for (size_t i = 0; i < refl.size(); i++) { + const int b = bin[i]; + const double eps = std::max(1, refl[i].epsilon); + e[i].eo = sigma_o[b] > 0 ? refl[i].f_obs / std::sqrt(eps * sigma_o[b]) : 0; + e[i].ec = sigma_c[b] > 0 ? refl[i].f_calc / std::sqrt(eps * sigma_c[b]) : 0; + e[i].centric = refl[i].centric; + } + + std::vector> free_in_bin(nbins); + std::vector free_all; + for (size_t i = 0; i < refl.size(); i++) + if (refl[i].free) { + free_in_bin[bin[i]].push_back(static_cast(i)); + free_all.push_back(static_cast(i)); + } + const double sigma_a_overall = MaximizeSigmaA(e, free_all); + std::vector sigma_a(nbins, sigma_a_overall); + for (int b = 0; b < nbins; b++) + if (static_cast(free_in_bin[b].size()) >= MIN_FREE_FOR_SHELL) + sigma_a[b] = MaximizeSigmaA(e, free_in_bin[b]); + result.sigma_a_lowest_shell = sigma_a.front(); + result.sigma_a_highest_shell = sigma_a.back(); + + double fom_sum = 0; + for (size_t i = 0; i < refl.size(); i++) { + const int b = bin[i]; + const double s = sigma_a[b]; + const double x = s * e[i].eo * e[i].ec / (1 - s * s); + result.weight[i].m = refl[i].centric ? std::tanh(x) : gemmi::bessel_i1_over_i0(2 * x); + // D takes Fc from its own shell scale to the observed one; the two are already close, because + // Fcalc reaches here scaled to Fobs overall and with a bulk-solvent and anisotropic-B model. + result.weight[i].d = sigma_c[b] > 0 ? s * std::sqrt(sigma_o[b] / sigma_c[b]) : s; + fom_sum += result.weight[i].m; + } + result.mean_fom = fom_sum / static_cast(refl.size()); + return result; +} diff --git a/rugnux/SigmaA.h b/rugnux/SigmaA.h new file mode 100644 index 000000000..afc403590 --- /dev/null +++ b/rugnux/SigmaA.h @@ -0,0 +1,42 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "gemmi/unitcell.hpp" + +// One reflection as the sigma_A estimate sees it: an observed and a calculated amplitude, already on +// a common scale, plus what the likelihood needs to know about the reflection itself. +struct SigmaAReflection { + double f_obs = 0; + double f_calc = 0; + double inv_d2 = 0; + int epsilon = 1; + bool centric = false; + bool free = false; +}; + +// The weights of one reflection's map coefficients: 2mFo - DFc and mFo - DFc. +struct SigmaAWeight { + double m = 1.0; // figure of merit + double d = 1.0; // D, carried on the observed amplitudes' scale +}; + +struct SigmaAResult { + std::vector weight; // one per input reflection, in the input order + int shells = 0; + int free_reflections = 0; + double mean_fom = 1.0; + double sigma_a_lowest_shell = 0.0; // for a log line: sigma_A at the low- and high-resolution ends + double sigma_a_highest_shell = 0.0; +}; + +// Per-resolution-shell sigma_A by maximum likelihood, and from it the figure of merit m and the +// scale D of each reflection. Estimated on the FREE reflections only - on the working set the model +// has been fitted to the data, so sigma_A would come out too high and the weighting would under-state +// the model's error where the map is meant to reveal it. +// +// Following Read (1986) Acta Cryst. A42, 140-149 +SigmaAResult EstimateSigmaA(const std::vector &refl, const gemmi::UnitCell &cell); diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index e7b8df034..1c96726c7 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -3,11 +3,18 @@ #include +#include +#include #include #include +#include + +#include #include "../common/Logger.h" #include "../rugnux/ModelValidation.h" +#include "../rugnux/RigidBodyRefine.h" +#include "../rugnux/SigmaA.h" namespace { // A synthetic P1 cell with two carbon atoms - enough for a reader to produce a Structure with @@ -45,6 +52,24 @@ ATOM 2 C . CB GLY A 1 12.000 14.000 16.000 1.00 20.00 "ATOM 2 CB GLY A 1 12.000 14.000 16.000 1.00 20.00 C\n" "END\n"; + // A synthetic "protein": carbons filling one asymmetric unit of a small P2(1)2(1)2(1) cell. The + // space group matters - in P1 the origin is free in all three directions, so |F| does not change + // when the whole content is translated and there is no translation to recover. No specimen is + // involved; the positions come out of a fixed seed. + std::string ClusterPdb() { + std::string pdb = "CRYST1 30.000 34.000 38.000 90.00 90.00 90.00 P 21 21 21 4\n"; + std::mt19937 rng(20260902); + std::uniform_real_distribution x(2, 14), y(2, 16), z(2, 18); + char line[96]; + for (int i = 1; i <= 150; i++) { + std::snprintf(line, sizeof line, + "ATOM %5d C UNK A 1 %8.3f%8.3f%8.3f 1.00 20.00 C\n", + i, x(rng), y(rng), z(rng)); + pdb += line; + } + return pdb + "END\n"; + } + std::string WriteTemp(const std::string &name, const char *content) { std::ofstream f(name); f << content; @@ -102,3 +127,87 @@ TEST_CASE("ModelValidation_UnusableModelGivesAReason", "[ModelValidation]") { std::filesystem::remove(empty_file); } + +// A rigid-body step is only worth taking if it can find a shift it was not told about. The check is +// closed: the "observed" amplitudes are the model's own, so the answer is known to be zero shift. +TEST_CASE("ModelValidation_RigidBodyRecoversASmallShift", "[ModelValidation]") { + Logger logger("ModelValidation_RigidBodyRecoversASmallShift"); + + const auto path = WriteTemp("rigid_body_test.pdb", ClusterPdb().c_str()); + gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect); + const gemmi::SpaceGroup *sg = st.find_spacegroup(); + REQUIRE(sg != nullptr); + st.setup_cell_images(); + + // "Observed" amplitudes: the model's own structure factors, so the target's minimum is exactly + // where the model started. + const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger); + REQUIRE_FALSE(ref.empty()); + gemmi::AsuData> fobs; + fobs.unit_cell_ = st.cell; + fobs.spacegroup_ = sg; + for (const auto &r : ref) + fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}}); + fobs.ensure_sorted(); + + const std::vector original = ModelPositions(st.models[0]); + std::vector displaced; + for (const gemmi::Position &p : original) + displaced.emplace_back(p.x + 0.40, p.y - 0.30, p.z + 0.20); // 0.54 A off + SetModelPositions(st.models[0], displaced); + + const RigidBodyRefineResult result = + RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger); + CHECK(result.converged); + + const std::vector refined = ModelPositions(st.models[0]); + double before = 0, after = 0; + for (size_t i = 0; i < original.size(); i++) { + before += original[i].dist_sq(displaced[i]); + after += original[i].dist_sq(refined[i]); + } + before = std::sqrt(before / original.size()); + after = std::sqrt(after / original.size()); + logger.Info("Rigid-body test: rmsd from the truth {:.3f} A -> {:.3f} A", before, after); + CHECK(after < 0.2 * before); + + std::filesystem::remove(path); +} + +// sigma_A is what says how much of the model to believe, so the two ends of its range are what the +// weighting has to get right: a model that explains the data completely, and one that explains none +// of it. +TEST_CASE("ModelValidation_SigmaAWeightsFollowTheModelsAgreement", "[ModelValidation]") { + gemmi::UnitCell cell(40, 50, 60, 90, 90, 90); + std::mt19937 rng(12345); + std::normal_distribution normal(0.0, 1.0); + + auto weights = [&](bool agreeing) { + std::vector refl; + for (int i = 0; i < 2000; i++) { + SigmaAReflection r; + r.f_calc = std::fabs(normal(rng)) * 100; + r.f_obs = agreeing ? r.f_calc : std::fabs(normal(rng)) * 100; + r.inv_d2 = 0.01 + 0.2 * (i / 2000.0); + r.free = (i % 20) == 0; // the usual 5 % + refl.push_back(r); + } + return EstimateSigmaA(refl, cell); + }; + + const SigmaAResult perfect = weights(true); + const SigmaAResult useless = weights(false); + CHECK(perfect.mean_fom > 0.85); + CHECK(useless.mean_fom < 0.2); + CHECK(perfect.shells == 2); // 100 free reflections, 50 to a shell + + // No free reflections to estimate on: the coefficients are left alone rather than weighted by a + // number that was never measured. + std::vector no_free; + for (int i = 0; i < 100; i++) + no_free.push_back({100.0, 100.0, 0.05, 1, false, false}); + const SigmaAResult unweighted = EstimateSigmaA(no_free, cell); + CHECK(unweighted.weight.size() == no_free.size()); + CHECK(unweighted.weight[0].m == 1.0); + CHECK(unweighted.weight[0].d == 1.0); +}