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); +} 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/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/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); +} 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);