From 0287e13945fc5b78f6027a3f137f5d5ad6b696ef Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 13 Nov 2026 08:11:03 +0100 Subject: [PATCH 1/3] MR: bring the placed copies together in the cell The search places each further copy wherever its translation function peaks, and every symmetry image of that position (an operator of the group plus a lattice shift) peaks alike, so the copies of a multi-copy solution came out at arbitrary images, often 50-80 A apart and in different cells, although the crystal they build is right. CompactCopies replaces each copy, in order, by the image whose centroid is nearest the mean centroid of the copies before it, and moves the whole by one lattice shift that puts the mean of the centroids in [0, 1). It runs once the decision is taken, on every candidate: within its own group the crystal is unchanged, while the check that carries the chosen copies into another group reads them as they stand, so compacting earlier would change that check's input. On a three-copy P 32 2 1 solution of an in-house cytochrome C set the copies went from 52-78 A apart to 26 A apart, each still matching a reference chain to 0.3-0.4 A C-alpha rms; the search's decisions were unchanged, and the held-out R-free before the validation's rigid body was identical (0.4331). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK --- .../MolecularReplacement.cpp | 61 +++++++++++++++- .../MolecularReplacement.h | 6 ++ tests/MolecularReplacementTest.cpp | 72 +++++++++++++++++++ 3 files changed, 138 insertions(+), 1 deletion(-) diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp index 6e8b3349e..087b0e94c 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.cpp +++ b/image_analysis/structure_refinement/MolecularReplacement.cpp @@ -2168,6 +2168,56 @@ std::vector SplitCopies(const gemmi::Structure &assembly, int return body; } +void CompactCopies(gemmi::Structure &placed, int copies) { + const gemmi::UnitCell &cell = placed.cell; + const std::vector ops = placed.find_spacegroup()->operations().all_ops_sorted(); + const auto apply = [](const gemmi::Op &op, const gemmi::Fractional &f, const gemmi::Fractional &shift) { + const std::array x = op.apply_to_xyz({{f.x, f.y, f.z}}); + return gemmi::Fractional(x[0] + shift.x, x[1] + shift.y, x[2] + shift.z); + }; + // Each copy in turn: the image of it - an operator and a lattice shift - whose centroid is nearest the mean + // centroid of the copies before it. + std::vector op(copies, gemmi::Op::identity()); + std::vector shift(copies); + const std::vector body = SplitCopies(placed, copies); + gemmi::Fractional sum; + for (int k = 0; k < copies; k++) { + const gemmi::Fractional c = centroid(body[k]); + gemmi::Fractional centre = c; + if (k > 0) { + const gemmi::Fractional target(sum * (1.0 / k)); + double nearest = INFINITY; + for (const gemmi::Op &o : ops) { + const gemmi::Fractional m = apply(o, c, gemmi::Fractional(0, 0, 0)); + for (int dx = -1; dx <= 1; dx++) + for (int dy = -1; dy <= 1; dy++) + for (int dz = -1; dz <= 1; dz++) { + const gemmi::Fractional s(std::round(target.x - m.x) + dx, std::round(target.y - m.y) + dy, + std::round(target.z - m.z) + dz); + const double d = cell.orthogonalize_difference(m + s - target).length(); + if (d < nearest) { + nearest = d; + op[k] = o; + shift[k] = s; + centre = m + s; + } + } + } + } + sum = sum + centre; + } + // The whole assembly then by one lattice shift, the mean of the centroids in the unit cell. + const gemmi::Fractional mean(sum * (1.0 / copies)); + const gemmi::Fractional to_cell(-std::floor(mean.x), -std::floor(mean.y), -std::floor(mean.z)); + const size_t chains = placed.models[0].chains.size(); + for (size_t i = 0; i < chains; i++) { + const size_t k = i * copies / chains; + for (gemmi::Residue &r : placed.models[0].chains[i].residues) + for (gemmi::Atom &a : r.atoms) + a.pos = cell.orthogonalize(apply(op[k], cell.fractionalize(a.pos), shift[k] + to_cell)); + } +} + std::optional SubgroupBasis(const gemmi::UnitCell &g, const gemmi::SpaceGroup &g_sg, const gemmi::UnitCell &h, const gemmi::SpaceGroup &h_sg) { constexpr double TOLERANCE = 0.01; // relative, on lengths, and on the cosines of the angles @@ -4030,7 +4080,16 @@ MolecularReplacementResult run_molecular_replacement(const RugnuxModel &input, c MolecularReplacementResult RunMolecularReplacement(const RugnuxModel &model, const std::vector &merges, Logger &logger, size_t nthreads, size_t gpu_resident_bytes) { WaitForModelWarmUp(); - return run_molecular_replacement(model, merges, logger, nthreads, gpu_resident_bytes); + MolecularReplacementResult result = run_molecular_replacement(model, merges, logger, nthreads, gpu_resident_bytes); + // Each placement's copies brought together once the decision is taken: the same crystal in its own group, but + // the copies carried into another group (carry_copies) are read as they stand, so not before. + for (MRCandidate &c : result.candidates) + if (c.placed && c.copies >= 1) { + auto compact = std::make_shared(*c.placed); + CompactCopies(*compact, c.copies); + c.placed = compact; + } + return result; } std::optional CarryMolecularReplacement(const MolecularReplacementResult &early, diff --git a/image_analysis/structure_refinement/MolecularReplacement.h b/image_analysis/structure_refinement/MolecularReplacement.h index 1ef05b244..f645099af 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.h +++ b/image_analysis/structure_refinement/MolecularReplacement.h @@ -263,6 +263,12 @@ gemmi::Structure AssembleCopies(const std::vector &copies); // The copies of an assembly (AssembleCopies), its chains split evenly in order. std::vector SplitCopies(const gemmi::Structure &assembly, int copies); +// The copies of a placement brought together, the crystal unchanged: the search puts each copy wherever the +// translation function peaks, which is any of its symmetry images alike. Each copy, in order, replaced by the image +// - an operator of the group and a lattice shift - whose centroid lies nearest the mean centroid of the copies +// before it; then the whole by one lattice shift that puts the mean of their centroids in the unit cell [0, 1). +void CompactCopies(gemmi::Structure &placed, int copies); + // The setting of a subgroup's merge (cell `h`, group `h_sg`) in the cell `g` of a group `g_sg`: `basis` has h's // axes as columns, lattice vectors of g (with g's centring) in g's fractional coordinates, right-handed and // matching h's lengths and angles to 1%; `origin` is where g's origin lies in h's fractional coordinates, so diff --git a/tests/MolecularReplacementTest.cpp b/tests/MolecularReplacementTest.cpp index dd92c098d..e51635d5f 100644 --- a/tests/MolecularReplacementTest.cpp +++ b/tests/MolecularReplacementTest.cpp @@ -613,6 +613,78 @@ TEST_CASE("MolecularReplacement_ExpandIntoSubgroup", "[MolecularReplacement]") { } } +// Three copies side by side, two of them then moved to far symmetry images - the crystal unchanged, as the search +// may place them: CompactCopies brings them back side by side, in the unit cell, and |F| stays as it was. +TEST_CASE("MolecularReplacement_CompactCopies", "[MolecularReplacement]") { + const SearchModel model = synthetic_model(); + const gemmi::UnitCell cell(70, 70, 80, 90, 90, 120); + const gemmi::SpaceGroup &sg = gemmi::get_spacegroup_by_name("P 32 2 1"); + const gemmi::Fractional t1(0.30, 0.20, 0.10), t2(0.62, 0.30, 0.12), t3(0.45, 0.45, 0.30); + const gemmi::Structure side_by_side = AssembleCopies({ + PlaceSearchModel(model, axis_angle({1, 2, 3}, 50), t1, cell, sg), + PlaceSearchModel(model, axis_angle({-2, 1, 1}, 110), t2, cell, sg), + PlaceSearchModel(model, axis_angle({0, 1, -1}, 30), t3, cell, sg)}); + + // Copy 2 by -x+y,-x,z+1/3 and (1, 0, -1), copy 3 by x-y,-y,-z+1/3 and (0, -1, 2) + gemmi::Structure scattered = side_by_side; + const auto move = [&](gemmi::Chain &ch, const gemmi::Op &op, const gemmi::Fractional &s) { + for (gemmi::Residue &r : ch.residues) + for (gemmi::Atom &a : r.atoms) { + const gemmi::Fractional f = cell.fractionalize(a.pos); + const std::array x = op.apply_to_xyz({{f.x, f.y, f.z}}); + a.pos = cell.orthogonalize(gemmi::Fractional(x[0] + s.x, x[1] + s.y, x[2] + s.z)); + } + }; + move(scattered.models[0].chains[1], gemmi::parse_triplet("-x+y,-x,z+1/3"), gemmi::Fractional(1, 0, -1)); + move(scattered.models[0].chains[2], gemmi::parse_triplet("x-y,-y,-z+1/3"), gemmi::Fractional(0, -1, 2)); + + const auto centres = [&](const gemmi::Structure &st) { + std::vector out; + for (const gemmi::Chain &ch : st.models[0].chains) { + gemmi::Vec3 c; + for (const gemmi::Atom &a : ch.residues[0].atoms) + c += a.pos; + out.push_back(cell.fractionalize(gemmi::Position(c * (1.0 / ch.residues[0].atoms.size())))); + } + return out; + }; + const auto widest = [&](const gemmi::Structure &st) { + const auto c = centres(st); + double d = 0; + for (size_t i = 0; i < c.size(); i++) + for (size_t j = i + 1; j < c.size(); j++) + d = std::max(d, cell.orthogonalize_difference(c[i] - c[j]).length()); + return d; + }; + REQUIRE(widest(side_by_side) < 30.0); + REQUIRE(widest(scattered) > 60.0); + + gemmi::Structure compact = scattered; + CompactCopies(compact, 3); + CHECK(widest(compact) < 30.0); + gemmi::Fractional mean; + for (const gemmi::Fractional &c : centres(compact)) + mean = mean + c; + for (double m : {mean.x / 3, mean.y / 3, mean.z / 3}) { + CHECK(m >= 0.0); + CHECK(m < 1.0); + } + + // The same crystal: the same amplitudes at every reflection + const auto before = observed(scattered); + const auto after = observed(compact); + REQUIRE(before.size() == after.size()); + double diff = 0, sum = 0; + for (size_t i = 0; i < before.size(); i++) { + REQUIRE(before[i].h == after[i].h); + REQUIRE(before[i].k == after[i].k); + REQUIRE(before[i].l == after[i].l); + diff += std::fabs(before[i].F - after[i].F); + sum += before[i].F; + } + CHECK(diff / sum < 1e-4); +} + TEST_CASE("MolecularReplacement_MatthewsProbabilities", "[MolecularReplacement]") { // Vm 2.5 A^3/Da for one copy: one copy is by far the most probable, four impossible (Vm 0.6). const auto p = MatthewsProbabilities(2.5 * 20000, 20000); From 564565ef45f65f5f0b4a146067e0dd6e8a637d89 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 13 Nov 2026 09:41:34 +0100 Subject: [PATCH 2/3] Difference-map blobs: a blob found again from its symmetry mate is dropped gemmi's flood fill marks the asymmetric-unit mates of a blob's points as it goes, but a difference map from a Fourier transform is symmetric only to ~0.06 sigma: an edge point of the mate that is just above the cutoff while its partner is just below it stays unmarked, seeds a second flood fill, and the whole mate is reported as a blob of its own. Moved beside the model, the two copies coincide. On a cytochrome C crystal with an AlphaFold model one of the three Fe sites was listed twice (21.5 sigma, 0.05 A apart); 11 of 19 runs of that data set carried one to three such pairs. A blob whose placed centroid is within 1.5 A of a stronger one is now left out. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK --- .../structure_refinement/DensityBlobs.cpp | 4 + tests/CMakeLists.txt | 1 + tests/DensityBlobsTest.cpp | 75 +++++++++++++++++++ 3 files changed, 80 insertions(+) create mode 100644 tests/DensityBlobsTest.cpp diff --git a/image_analysis/structure_refinement/DensityBlobs.cpp b/image_analysis/structure_refinement/DensityBlobs.cpp index b2a9b7fb9..81063d300 100644 --- a/image_analysis/structure_refinement/DensityBlobs.cpp +++ b/image_analysis/structure_refinement/DensityBlobs.cpp @@ -105,6 +105,10 @@ std::vector FindDensityBlobs(const gemmi::Grid &difference, cra.residue->seqid.str(); } else d.nearest_atom_distance = INFINITY; + // The flood fill can find a blob that crosses the asymmetric unit's boundary again from a symmetry mate; + // beside the model the copies coincide, and the first (strongest) is kept. + if (std::any_of(out.begin(), out.end(), [&](const DensityBlob &o) { return o.centroid.dist(d.centroid) < 1.5; })) + continue; ns.for_each(d.centroid, '\0', 3.4, [&](gemmi::NeighborSearch::Mark &m, double dist_sq) { const double r = std::sqrt(dist_sq); if (Polar(m.element) && r >= 2.4) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 668a5e715..0bc1bafd5 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -196,6 +196,7 @@ ADD_EXECUTABLE(jfjoch_test RefinementTest.cpp RefineOnlyTest.cpp SolventSitesTest.cpp + DensityBlobsTest.cpp SolventBuildTest.cpp SearchSpaceGroupTest.cpp SearchSpaceGroupTwinTest.cpp diff --git a/tests/DensityBlobsTest.cpp b/tests/DensityBlobsTest.cpp new file mode 100644 index 000000000..8019a1c6a --- /dev/null +++ b/tests/DensityBlobsTest.cpp @@ -0,0 +1,75 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include + +#include +#include +#include + +#include "../image_analysis/structure_refinement/DensityBlobs.h" + +namespace { + +// A P 32 2 1 map in sigma units with one Gaussian peak at `f` and its symmetry mates, the mates equal to within +// 0.1 sigma as in a map from a Fourier transform (the live cytochrome C difference map: up to 0.06 sigma) +gemmi::Grid PeakMap(const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg, const gemmi::Fractional &f, + std::mt19937 &rng) { + gemmi::Grid grid; + grid.unit_cell = cell; + grid.spacegroup = sg; + grid.set_size(72, 72, 84); + size_t idx = 0; + for (int w = 0; w < grid.nw; w++) + for (int v = 0; v < grid.nv; v++) + for (int u = 0; u < grid.nu; u++, idx++) { + gemmi::Fractional d(double(u) / grid.nu - f.x, double(v) / grid.nv - f.y, double(w) / grid.nw - f.z); + d = gemmi::Fractional(d.x - std::round(d.x), d.y - std::round(d.y), d.z - std::round(d.z)); + const double r2 = cell.orthogonalize_difference(d).length_sq(); + grid.data[idx] = static_cast(10.0 * std::exp(-r2 / 2.0)); + } + grid.symmetrize_max(); + std::normal_distribution noise(0.0f, 0.1f); + for (float &v : grid.data) + v += noise(rng); + return grid; +} + +} + +TEST_CASE("FindDensityBlobs_SymmetryMatesOnce", "[DensityBlobs]") { + const gemmi::UnitCell cell(60, 60, 70, 90, 90, 120); + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 32 2 1"); + REQUIRE(sg != nullptr); + gemmi::BlobCriteria criteria; + criteria.cutoff = 3.0; + criteria.min_volume = 0.0; + criteria.min_score = 0.0; + criteria.min_peak = 4.0; + + // Peaks over the whole cell, the asymmetric unit's faces among them: gemmi's flood fill finds some of those + // that cross a face twice - an edge point above the cutoff whose mate is just below it seeds the mate's copy. + std::mt19937 rng(1); + int found_twice = 0; + for (int i = 0; i < 8; i++) + for (int j = 0; j < 8; j++) + for (int k = 0; k < 8; k++) { + const gemmi::Fractional f(i / 8.0 + 0.01, j / 8.0 + 0.02, k / 8.0 + 0.005); + const gemmi::Grid map = PeakMap(cell, sg, f, rng); + const gemmi::Position p = cell.orthogonalize(f); + const std::string pdb = fmt::format( + "CRYST1 60.000 60.000 70.000 90.00 90.00 120.00 P 32 2 1\n" + "HETATM 1 O HOH A 1 {:8.3f}{:8.3f}{:8.3f} 1.00 20.00 O\nEND\n", + p.x + 2.5, p.y, p.z); + const gemmi::Structure st = gemmi::read_pdb_string(pdb, "one water"); + if (gemmi::find_blobs_by_flood_fill(map, criteria).size() > 1) + found_twice++; + const std::vector blobs = FindDensityBlobs(map, 1.0, st, 2.0); + INFO(fmt::format("peak at {:.3f} {:.3f} {:.3f}", f.x, f.y, f.z)); + CHECK(blobs.size() == 1); + } + CHECK(found_twice > 0); +} From 6573f48444084dd580536360d960909fe58dca4e Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 13 Nov 2026 09:41:24 +0100 Subject: [PATCH 3/3] Rugnux live structure view: step through the Fo-Fc peaks A bar over the Mol* stage carries previous/next arrows through the refined model's difference-map blobs (already in the live result, highest peak first), with the peak's index, height and reading, and a button back to the whole model; the Fo-Fc toggle moves into the same bar. At a peak the camera looks at it from the side with the fewest atoms in front, 26 A away, with the clip planes 5 A either side of it; the maps are drawn uncarved within 10 A of it (the carved maps about the model would hide unmodelled density such as a missing cofactor) and the model as sticks over the cartoon. Smaller view fixes: stick carbons take the chain colours instead of Mol*'s default greens, which read as Fo-Fc density; riding hydrogens are not drawn; the 2Fo-Fc mesh is dimmed (alpha 0.35 over the whole model, 0.5 at a peak) so the cartoon and the difference density show through it. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK --- frontend/src/components/RugnuxPage.tsx | 2 +- .../src/components/RugnuxStructureViewer.tsx | 188 ++++++++++++++---- frontend/src/components/structure/loadMaps.ts | 19 +- .../src/components/structure/loadModel.ts | 9 + frontend/src/components/structure/peakView.ts | 60 ++++++ 5 files changed, 234 insertions(+), 44 deletions(-) create mode 100644 frontend/src/components/structure/peakView.ts diff --git a/frontend/src/components/RugnuxPage.tsx b/frontend/src/components/RugnuxPage.tsx index 38f9d1cef..9e685047e 100644 --- a/frontend/src/components/RugnuxPage.tsx +++ b/frontend/src/components/RugnuxPage.tsx @@ -141,7 +141,7 @@ function Centre({s, pal}: { s: rugnux_live_status, pal: LivePalette }) { }> - diff --git a/frontend/src/components/RugnuxStructureViewer.tsx b/frontend/src/components/RugnuxStructureViewer.tsx index 093381aa4..ba8ebc6d5 100644 --- a/frontend/src/components/RugnuxStructureViewer.tsx +++ b/frontend/src/components/RugnuxStructureViewer.tsx @@ -1,13 +1,18 @@ import React, {useCallback, useEffect, useRef, useState} from "react"; -import {Alert, Box, Tooltip, Typography} from "@mui/material"; +import {Alert, Box, IconButton, Tooltip, Typography} from "@mui/material"; +import ChevronLeftIcon from "@mui/icons-material/ChevronLeft"; +import ChevronRightIcon from "@mui/icons-material/ChevronRight"; +import ZoomOutMapIcon from "@mui/icons-material/ZoomOutMap"; import type {PluginContext} from "molstar/lib/mol-plugin/context"; import {authHeaders, noteResponse} from "../auth"; +import type {rugnux_live_blob} from "../client"; import MolstarCanvas from "./structure/MolstarCanvas"; import {applyModelDisplay, clearScene, getStructureAtoms, getStructureBox, LoadedModel, loadModel, ReprMode} from "./structure/loadModel"; -import {loadMapLayer, MAP_COLORS, MapLayer, setLayerVisible, updateLayerSigma} from "./structure/loadMaps"; -import {fullCellFromCcp4, prepareMap} from "./structure/densityCrop"; +import {loadMapLayer, MAP_COLORS, MapLayer, removeMapLayer, setLayerVisible} from "./structure/loadMaps"; +import {FullCellMap, fullCellFromCcp4, prepareMap} from "./structure/densityCrop"; import {revealCamera} from "./structure/reveal"; +import {PEAK_MAP_RADIUS_A, viewPeak} from "./structure/peakView"; // A glance at the model rugnux live placed - refined, where the run refined it - in its density: cartoon, @@ -15,11 +20,25 @@ import {revealCamera} from "./structure/reveal"; // map at 4 sigma, in Mol*. Set up as aaredb's structure viewer (gitea.psi.ch/mx/aaredb, pages/StructureView.tsx // and components/Structure/, the PSI MX group's own code): its plugin, Coot mouse controls, model presets, map // meshes, crop and carve - on a black stage. It is not a viewer to work in: what it shows is gone with the next -// collection, so the only control is a toggle for the difference map; every map is carved within the same -// radius of the atoms (aaredb carves 2Fo-Fc only), and the first model makes an entrance (reveal.ts). Mol* is -// only loaded with this component (React.lazy in RugnuxPage), and every Mol* import lives in ./structure. +// collection, so the controls are a toggle for the difference map and arrows through the refined model's +// difference-map peaks (the blobs of the live result, highest first), in a bar over the stage. Every map is carved +// within the same radius of the atoms (aaredb carves 2Fo-Fc only); at a peak the maps are drawn uncarved within +// a sphere about it instead, and the camera frames it in a slab (peakView.ts). The first model makes an entrance +// (reveal.ts). Mol* is only loaded with this component (React.lazy in RugnuxPage), and every Mol* import lives in +// ./structure. type Kind = 'twofofc' | 'fofc' | 'anom'; +const MAP_NAME: Record = {twofofc: '2fofc', fofc: 'fofc', anom: 'anom'}; +// The 2Fo-Fc mesh dimmed, so that the cartoon (whole model) and the difference density (at a peak) show through it +const TWOFOFC_ALPHA = {model: 0.35, peak: 0.5}; +const KIND_WORD: Record = { + ION: 'ion', WATER: 'water', LIGAND: 'ligand', MODEL_ERROR: 'model error', UNASSIGNED: 'unassigned'}; + +// ["peak 1/20 · 23.4 σ", "ion · NE2 HIS A19 1.9 Å"] +function peakLabel(p: rugnux_live_blob, i: number, n: number): [string, string] { + const near = p.nearest_atom_distance_A != null ? `${p.nearest_atom} ${p.nearest_atom_distance_A.toFixed(1)} Å` : ''; + return [`peak ${i + 1}/${n} · ${p.peak_sigma.toFixed(1)} σ`, [KIND_WORD[p.kind], near].filter(Boolean).join(' · ')]; +} const hex = (c: number) => '#' + c.toString(16).padStart(6, '0'); @@ -39,9 +58,11 @@ const nextFrames = () => new Promise(r => requestAnimationFrame(() => requestAni // resultKey: changes when a new result is there to load; anomalous: it has an anomalous difference map; // anomalousSignal: the data carry usable anomalous signal, so that map is worth showing from the start; -// fill: take the parent's height (the live page's stage) instead of the viewer's own. -function RugnuxStructureViewer({resultKey, anomalous, anomalousSignal, fill}: - { resultKey: string, anomalous: boolean, anomalousSignal: boolean, fill?: boolean }) { +// fill: take the parent's height (the live page's stage) instead of the viewer's own; peaks: the blobs of the +// refined model's difference map, in the model's frame. +function RugnuxStructureViewer({resultKey, anomalous, anomalousSignal, fill, peaks}: + { resultKey: string, anomalous: boolean, anomalousSignal: boolean, fill?: boolean, + peaks?: rugnux_live_blob[] }) { const pluginRef = useRef(null); const [pluginReady, setPluginReady] = useState(0); const [error, setError] = useState(); @@ -49,7 +70,19 @@ function RugnuxStructureViewer({resultKey, anomalous, anomalousSignal, fill}: const [mapsLoading, setMapsLoading] = useState(false); const layersRef = useRef>>({}); const modelRefsRef = useRef({}); - const sigmaTimer = useRef | null>(null); + // The full-cell maps and the atoms, for the maps about a peak; the layers drawn about the peak shown + const fullRef = useRef>>({}); + const atomsRef = useRef(new Float32Array(0)); + const regionRef = useRef>>({}); + // The peak shown (index into `sorted`), null for the whole model. Steps run one after another; a step + // overtaken by a later click is skipped. + const [peak, setPeak] = useState(null); + const peakRef = useRef(null); + const stepQueue = useRef>(Promise.resolve()); + const stepToken = useRef(0); + const sorted = [...(peaks ?? [])].sort((a, b) => b.peak_sigma - a.peak_sigma); + const sortedRef = useRef(sorted); + sortedRef.current = sorted; // First load frames the molecule; later results keep the user's camera and clip (aaredb StructureView.tsx:117-121) const sceneLoadedRef = useRef(false); const [hasMap, setHasMap] = useState>({twofofc: false, fofc: false, anom: false}); @@ -80,8 +113,13 @@ function RugnuxStructureViewer({resultKey, anomalous, anomalousSignal, fill}: setMapsLoading(false); setError(undefined); try { + stepToken.current++; + peakRef.current = null; + setPeak(null); await clearScene(plugin); layersRef.current = {}; + regionRef.current = {}; + fullRef.current = {}; setHasMap({twofofc: false, fofc: false, anom: false}); const [cif, map2, map1, mapA] = await Promise.all([ fetchBytes("/rugnux/model.cif", abort.signal), @@ -102,6 +140,7 @@ function RugnuxStructureViewer({resultKey, anomalous, anomalousSignal, fill}: const box = getStructureBox(loaded.structure); const crop = box ? {...box, marginA: 6} : null; const atoms = getStructureAtoms(loaded.structure); + atomsRef.current = atoms; setModelLoading(false); setMapsLoading(true); const maps: [Kind, '2fofc' | 'fofc' | 'anom', ArrayBuffer | undefined][] = @@ -109,8 +148,10 @@ function RugnuxStructureViewer({resultKey, anomalous, anomalousSignal, fill}: for (const [kind, name, bytes] of maps) { if (!bytes || abort.signal.aborted) continue; - const computed = prepareMap(name, fullCellFromCcp4(bytes), crop, atoms.length ? atoms : null); - const layer = await loadMapLayer(plugin, computed, name, SIGMA[kind]); + const full = fullCellFromCcp4(bytes); + fullRef.current[kind] = full; + const computed = prepareMap(name, full, crop, atoms.length ? atoms : null); + const layer = await loadMapLayer(plugin, computed, name, SIGMA[kind], TWOFOFC_ALPHA.model); layersRef.current[kind] = layer; if (!showRef.current[kind]) setLayerVisible(plugin, layer, false); @@ -138,38 +179,111 @@ function RugnuxStructureViewer({resultKey, anomalous, anomalousSignal, fill}: const toggleLayer = useCallback((kind: Kind, visible: boolean) => { setShow(s => ({...s, [kind]: visible})); const plugin = pluginRef.current; - const layer = layersRef.current[kind]; + const layer = peakRef.current === null ? layersRef.current[kind] : regionRef.current[kind]; if (plugin && layer) setLayerVisible(plugin, layer, visible); }, []); + // Show peak i (null: the whole model): the maps about the model give way to the maps about the peak, and the + // camera frames it. + const showPeak = useCallback(async (i: number | null) => { + const plugin = pluginRef.current; + if (!plugin) + return; + for (const layer of Object.values(regionRef.current)) + await removeMapLayer(plugin, layer); + regionRef.current = {}; + for (const [kind, layer] of Object.entries(layersRef.current) as [Kind, MapLayer][]) + setLayerVisible(plugin, layer, i === null && showRef.current[kind]); + // at a peak the side chains about it matter: sticks over the cartoon + applyModelDisplay(plugin, modelRefsRef.current, i === null ? REPR : 'both', false); + if (i === null) { + plugin.managers.camera.reset(undefined, reducedMotion() ? 0 : 400); + return; + } + const p = sortedRef.current[i]; + const centre = [p.x, p.y, p.z]; + viewPeak(plugin, centre, atomsRef.current, reducedMotion() ? 0 : 400); + for (const [kind, full] of Object.entries(fullRef.current) as [Kind, FullCellMap][]) { + const map = prepareMap(MAP_NAME[kind], full, {min: centre, max: centre, marginA: PEAK_MAP_RADIUS_A}, + Float32Array.from(centre), PEAK_MAP_RADIUS_A); + const layer = await loadMapLayer(plugin, map, MAP_NAME[kind], SIGMA[kind], TWOFOFC_ALPHA.peak); + regionRef.current[kind] = layer; + if (!showRef.current[kind]) + setLayerVisible(plugin, layer, false); + } + }, []); + + const step = useCallback((i: number | null) => { + const token = ++stepToken.current; + peakRef.current = i; + setPeak(i); + stepQueue.current = stepQueue.current + .then(() => token === stepToken.current ? showPeak(i) : undefined) + .catch(e => setError(e instanceof Error ? e.message : String(e))); + }, [showPeak]); + + const n = sorted.length; + const busy = modelLoading || mapsLoading; + const [label, detail] = peak === null + ? [`${n} Fo−Fc peak${n === 1 ? '' : 's'}`, `highest ${sorted[0]?.peak_sigma.toFixed(1)} σ`] + : peakLabel(sorted[peak], peak, n); + const arrow = {color: '#c9d3e6', '&.Mui-disabled': {color: 'rgba(154,165,196,.3)'}}; + return ( - - - - - {(modelLoading || (mapsLoading && !revealed)) && !error && - - {modelLoading ? 'loading model…' : 'computing density maps…'} - } - {error && {error}} - {/* the one control: the difference map on or off */} - {hasMap.fofc && ( - - toggleLayer('fofc', !show.fofc)} disabled={mapsLoading} - sx={{position: 'absolute', top: 8, right: 8, display: 'inline-flex', alignItems: 'center', gap: 0.75, - px: 1, height: 26, borderRadius: '13px', border: '1px solid', cursor: 'pointer', - borderColor: show.fofc ? 'primary.main' : 'rgba(154,165,196,.3)', bgcolor: 'rgba(8,12,26,.8)', - color: show.fofc ? '#e9edf7' : '#9aa5c4', fontSize: '0.75rem', fontFamily: 'inherit'}}> - - Fo−Fc + {/* the bar over the stage: the difference-map peaks, highest first, and the difference map on or off */} + + + {n > 0 && <> + step(peak! - 1)}> + + + + {label} + {detail} - )} - {!modelLoading && !mapsLoading && !hasMap.twofofc && !error && - No maps} + step(peak === null ? 0 : peak + 1)}> + + + {peak !== null && + + step(null)}> + + + } + } + + {hasMap.fofc && ( + + toggleLayer('fofc', !show.fofc)} disabled={mapsLoading} + sx={{display: 'inline-flex', alignItems: 'center', gap: 0.75, flexShrink: 0, ml: 'auto', mr: 0.5, my: 0.5, + px: 1, height: 26, borderRadius: '13px', border: '1px solid', cursor: 'pointer', + borderColor: show.fofc ? 'primary.main' : 'rgba(154,165,196,.3)', bgcolor: 'rgba(8,12,26,.8)', + color: show.fofc ? '#e9edf7' : '#9aa5c4', fontSize: '0.75rem', fontFamily: 'inherit'}}> + + Fo−Fc + + )} + + + + + + {(modelLoading || (mapsLoading && !revealed)) && !error && + + {modelLoading ? 'loading model…' : 'computing density maps…'} + } + {error && {error}} + {!modelLoading && !mapsLoading && !hasMap.twofofc && !error && + No maps} + ); } diff --git a/frontend/src/components/structure/loadMaps.ts b/frontend/src/components/structure/loadMaps.ts index b15c64bf0..6d1ebed33 100644 --- a/frontend/src/components/structure/loadMaps.ts +++ b/frontend/src/components/structure/loadMaps.ts @@ -21,6 +21,7 @@ export const MAP_COLORS = { } as const; export interface MapLayer { + dataRef: string; // the raw bytes the volume is parsed from: removing it removes the layer volumeRef: string; reprRefs: string[]; // one for 2fofc and anom, two (+/−) for fofc rmsd: number; @@ -36,6 +37,7 @@ async function addIso( map: ComputedMap, sigma: number, color: number, + alpha: number, ): Promise { const volumeCell = plugin.state.data.select(volumeRef)[0]; const repr = plugin @@ -45,7 +47,7 @@ async function addIso( StateTransforms.Representation.VolumeRepresentation3D, createVolumeRepresentationParams(plugin, volumeCell?.obj?.data, { type: 'isosurface', - typeParams: { isoValue: abs(map, sigma), visuals: ['wireframe'], sizeFactor: 1 }, + typeParams: { isoValue: abs(map, sigma), visuals: ['wireframe'], sizeFactor: 1, alpha }, color: 'uniform', colorParams: { value: Color(color) }, }), @@ -68,6 +70,7 @@ export async function loadMapLayer( map: ComputedMap, kind: '2fofc' | 'fofc' | 'anom', sigma: number, + alpha = 1, // of the 2Fo-Fc mesh, which can be dimmed so the model and the difference map show through it ): Promise { const data = await plugin.builders.data.rawData({ data: new Uint8Array(map.ccp4), @@ -78,14 +81,14 @@ export async function loadMapLayer( const reprRefs: string[] = []; if (kind === '2fofc') { - reprRefs.push(await addIso(plugin, volumeRef, map, sigma, MAP_COLORS.twofofc)); + reprRefs.push(await addIso(plugin, volumeRef, map, sigma, MAP_COLORS.twofofc, alpha)); } else if (kind === 'anom') { - reprRefs.push(await addIso(plugin, volumeRef, map, sigma, MAP_COLORS.anom)); + reprRefs.push(await addIso(plugin, volumeRef, map, sigma, MAP_COLORS.anom, 1)); } else { - reprRefs.push(await addIso(plugin, volumeRef, map, +sigma, MAP_COLORS.fofcPos)); - reprRefs.push(await addIso(plugin, volumeRef, map, -sigma, MAP_COLORS.fofcNeg)); + reprRefs.push(await addIso(plugin, volumeRef, map, +sigma, MAP_COLORS.fofcPos, 1)); + reprRefs.push(await addIso(plugin, volumeRef, map, -sigma, MAP_COLORS.fofcNeg, 1)); } - return { volumeRef, reprRefs, rmsd: map.rmsd, mean: map.mean }; + return { dataRef: data.ref, volumeRef, reprRefs, rmsd: map.rmsd, mean: map.mean }; } /** σ slider update: param-only change, re-mesh without recompute. */ @@ -107,3 +110,7 @@ export async function updateLayerSigma( export function setLayerVisible(plugin: PluginContext, layer: MapLayer, visible: boolean): void { setSubtreeVisibility(plugin.state.data, layer.volumeRef, !visible); // arg = isHidden } + +export async function removeMapLayer(plugin: PluginContext, layer: MapLayer): Promise { + await plugin.build().delete(layer.dataRef).commit(); +} diff --git a/frontend/src/components/structure/loadModel.ts b/frontend/src/components/structure/loadModel.ts index 1d10d28da..8779a12bb 100644 --- a/frontend/src/components/structure/loadModel.ts +++ b/frontend/src/components/structure/loadModel.ts @@ -15,6 +15,10 @@ const CHAIN_PALETTE = { }, }; +// Sticks: carbons in the chain's colour, not Mol*'s default palette, whose greens read as Fo-Fc density; the +// riding hydrogens a refined model carries are left out. +const STICK_COLORS = { carbonColor: { name: 'chain-id' as const, params: { ...CHAIN_PALETTE, asymId: 'auth' as const } } }; + // Load a PDB model from bytes into the plugin's state tree. All representations // are built up front with stable refs; "representation mode" and waters are // pure visibility flips (setSubtreeVisibility), so switching is instant. @@ -55,7 +59,9 @@ export async function loadModel( refs.polymerCartoon = cartoon.ref; const sticks = await plugin.builders.structure.representation.addRepresentation(polymer, { type: 'ball-and-stick', + typeParams: { ignoreHydrogens: true }, color: 'element-symbol', + colorParams: STICK_COLORS, }); refs.polymerSticks = sticks.ref; } @@ -63,6 +69,9 @@ export async function loadModel( if (ligand) { const r = await plugin.builders.structure.representation.addRepresentation(ligand, { type: 'ball-and-stick', + typeParams: { ignoreHydrogens: true }, + color: 'element-symbol', + colorParams: STICK_COLORS, }); refs.ligand = r.ref; } diff --git a/frontend/src/components/structure/peakView.ts b/frontend/src/components/structure/peakView.ts new file mode 100644 index 000000000..143ab4876 --- /dev/null +++ b/frontend/src/components/structure/peakView.ts @@ -0,0 +1,60 @@ +import { PluginContext } from 'molstar/lib/mol-plugin/context'; +import { Vec3 } from 'molstar/lib/mol-math/linear-algebra'; + +// Framing one difference-map peak: the camera looks at it from the side where the least of the model lies in +// front of it, close enough for the blob and what surrounds it (a haem is ~12 A across) to fill the view, and the +// clip planes cut a slab through it so that nothing further in front hides it and nothing further behind clutters it. +export const PEAK_DISTANCE_A = 26; // camera to peak: a 45 degree field of view shows ~11 A above and below it +export const PEAK_SLAB_A = 5; // the near and far clip planes this far in front of and behind the peak +export const PEAK_MAP_RADIUS_A = 10; // the maps are drawn within this radius of the peak, uncarved + +// The 26 directions to the corners, edges and faces of a cube +const DIRECTIONS: Vec3[] = []; +for (let x = -1; x <= 1; x++) + for (let y = -1; y <= 1; y++) + for (let z = -1; z <= 1; z++) + if (x || y || z) + DIRECTIONS.push(Vec3.normalize(Vec3(), Vec3.create(x, y, z))); + +// Atoms between the peak and the near clip plane, within 3 A of the line of sight +function atomsInFront(atoms: Float32Array, c: Vec3, dir: Vec3): number { + let n = 0; + for (let i = 0; i < atoms.length; i += 3) { + const dx = atoms[i] - c[0], dy = atoms[i + 1] - c[1], dz = atoms[i + 2] - c[2]; + const t = dx * dir[0] + dy * dir[1] + dz * dir[2]; + if (t > 1.5 && t < PEAK_SLAB_A && dx * dx + dy * dy + dz * dz - t * t < 9) + n++; + } + return n; +} + +export function viewPeak(plugin: PluginContext, center: number[], atoms: Float32Array, durationMs: number): void { + const camera = plugin.canvas3d?.camera; + if (!camera) return; + const c = Vec3.create(center[0], center[1], center[2]); + // the current direction first, so that it is kept where no other is clearer + const current = Vec3.normalize(Vec3(), Vec3.sub(Vec3(), camera.state.position, camera.state.target)); + let best = current; + let bestCount = atomsInFront(atoms, c, current); + let bestDot = 1; + for (const d of DIRECTIONS) { + const n = atomsInFront(atoms, c, d); + const dot = Vec3.dot(d, current); + if (n < bestCount || (n === bestCount && dot > bestDot)) { + best = d; + bestCount = n; + bestDot = dot; + } + } + // "up" as close to the current one as the new direction allows + const up = Vec3.scaleAndAdd(Vec3(), camera.state.up, best, -Vec3.dot(camera.state.up, best)); + if (Vec3.magnitude(up) < 1e-3) + Vec3.cross(up, best, Math.abs(best[0]) < 0.9 ? Vec3.unitX : Vec3.unitY); + Vec3.normalize(up, up); + camera.setState({ + target: c, + position: Vec3.scaleAndAdd(Vec3(), c, best, PEAK_DISTANCE_A), + up, + radius: PEAK_SLAB_A, + }, durationMs); +}