Merge branch 'rc187-molstar-peaks' into rc186-merge2
This commit is contained in:
leonarski_f committed
commit
e14f3a0761
11 files changed
+452
-45
No files matched your search
@@ -141,7 +141,7 @@ function Centre({s, pal}: { s: rugnux_live_status, pal: LivePalette }) {
|
||||
<Division sx={{flexGrow: 1}} label={s.refined_model ? "Model · refined" : "Model · placed"}>
|
||||
<Box sx={{flexGrow: 1, minHeight: 0}}>
|
||||
<Suspense fallback={<Later text="Loading viewer..."/>}>
|
||||
<RugnuxStructureViewer resultKey={resultKey} anomalous={s.anomalous_map} fill
|
||||
<RugnuxStructureViewer resultKey={resultKey} anomalous={s.anomalous_map} fill peaks={m.refinement?.blobs}
|
||||
anomalousSignal={s.anomalous_signal}/>
|
||||
</Suspense>
|
||||
</Box>
|
||||
|
||||
@@ -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<Kind, '2fofc' | 'fofc' | 'anom'> = {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<rugnux_live_blob['kind'], string> = {
|
||||
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<PluginContext | null>(null);
|
||||
const [pluginReady, setPluginReady] = useState(0);
|
||||
const [error, setError] = useState<string>();
|
||||
@@ -49,7 +70,19 @@ function RugnuxStructureViewer({resultKey, anomalous, anomalousSignal, fill}:
|
||||
const [mapsLoading, setMapsLoading] = useState(false);
|
||||
const layersRef = useRef<Partial<Record<Kind, MapLayer>>>({});
|
||||
const modelRefsRef = useRef<LoadedModel['refs']>({});
|
||||
const sigmaTimer = useRef<ReturnType<typeof setTimeout> | 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<Partial<Record<Kind, FullCellMap>>>({});
|
||||
const atomsRef = useRef<Float32Array>(new Float32Array(0));
|
||||
const regionRef = useRef<Partial<Record<Kind, MapLayer>>>({});
|
||||
// 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<number | null>(null);
|
||||
const peakRef = useRef<number | null>(null);
|
||||
const stepQueue = useRef<Promise<void>>(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<Record<Kind, boolean>>({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 (
|
||||
<Box sx={{position: 'relative', bgcolor: '#000', borderRadius: 1, overflow: 'hidden',
|
||||
<Box sx={{display: 'flex', flexDirection: 'column', bgcolor: '#000', borderRadius: 1, overflow: 'hidden',
|
||||
height: fill ? '100%' : {xs: '55vh', md: 'clamp(520px, 64vh, 820px)'}, minHeight: fill ? 0 : undefined}}>
|
||||
<Box sx={{position: 'absolute', inset: 0, opacity: revealed ? 1 : 0,
|
||||
transition: 'opacity 900ms ease-out', '@media (prefers-reduced-motion: reduce)': {transition: 'none'}}}>
|
||||
<MolstarCanvas onReady={onPluginReady}/>
|
||||
</Box>
|
||||
{(modelLoading || (mapsLoading && !revealed)) && !error && <Box sx={{position: 'absolute', inset: 0, display: 'flex', alignItems: 'center',
|
||||
justifyContent: 'center', pointerEvents: 'none'}}>
|
||||
<Typography sx={{fontSize: '0.875rem', color: '#8aa0b5'}}>
|
||||
{modelLoading ? 'loading model…' : 'computing density maps…'}</Typography>
|
||||
</Box>}
|
||||
{error && <Alert severity="warning" sx={{position: 'absolute', top: 8, left: 8}}>{error}</Alert>}
|
||||
{/* the one control: the difference map on or off */}
|
||||
{hasMap.fofc && (
|
||||
<Tooltip title={`mFo−DFc at ±${SIGMA.fofc.toFixed(0)} σ (green positive, red negative); 2mFo−DFc at ${SIGMA.twofofc.toFixed(0)} σ in blue`}>
|
||||
<Box component="button" type="button" onClick={() => 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'}}>
|
||||
<Box sx={{width: 8, height: 8, borderRadius: '2px', bgcolor: hex(MAP_COLORS.fofcPos)}}/>
|
||||
Fo−Fc
|
||||
{/* the bar over the stage: the difference-map peaks, highest first, and the difference map on or off */}
|
||||
<Box sx={{display: 'flex', flexWrap: 'wrap', alignItems: 'center', columnGap: 0.5, px: 0.5, minHeight: 40,
|
||||
flexShrink: 0, bgcolor: '#0b1020', borderBottom: '1px solid rgba(154,165,196,.2)'}}>
|
||||
<Box sx={{display: 'flex', alignItems: 'center', gap: 0.5, flex: '1 1 250px', minWidth: 0}}>
|
||||
{n > 0 && <>
|
||||
<IconButton size="small" aria-label="previous peak" sx={arrow}
|
||||
disabled={busy || peak === null || peak === 0} onClick={() => step(peak! - 1)}>
|
||||
<ChevronLeftIcon fontSize="small"/>
|
||||
</IconButton>
|
||||
<Box title={`${label} · ${detail}`} sx={{minWidth: 0, lineHeight: 1.25, fontVariantNumeric: 'tabular-nums'}}>
|
||||
<Typography noWrap sx={{fontSize: '0.8rem', fontWeight: 600, color: '#e9edf7'}}>{label}</Typography>
|
||||
<Typography noWrap sx={{fontSize: '0.72rem', color: '#9aa5c4'}}>{detail}</Typography>
|
||||
</Box>
|
||||
</Tooltip>)}
|
||||
{!modelLoading && !mapsLoading && !hasMap.twofofc && !error &&
|
||||
<Typography sx={{position: 'absolute', bottom: 8, right: 12, fontSize: '0.8rem', color: '#8aa0b5'}}>No maps</Typography>}
|
||||
<IconButton size="small" aria-label="next peak" sx={arrow}
|
||||
disabled={busy || peak === n - 1} onClick={() => step(peak === null ? 0 : peak + 1)}>
|
||||
<ChevronRightIcon fontSize="small"/>
|
||||
</IconButton>
|
||||
{peak !== null &&
|
||||
<Tooltip title="Whole model">
|
||||
<IconButton size="small" aria-label="whole model" sx={arrow} onClick={() => step(null)}>
|
||||
<ZoomOutMapIcon fontSize="small"/>
|
||||
</IconButton>
|
||||
</Tooltip>}
|
||||
</>}
|
||||
</Box>
|
||||
{hasMap.fofc && (
|
||||
<Tooltip title={`mFo−DFc at ±${SIGMA.fofc.toFixed(0)} σ (green positive, red negative); 2mFo−DFc at ${SIGMA.twofofc.toFixed(0)} σ in blue`}>
|
||||
<Box component="button" type="button" onClick={() => 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'}}>
|
||||
<Box sx={{width: 8, height: 8, borderRadius: '2px', bgcolor: hex(MAP_COLORS.fofcPos)}}/>
|
||||
Fo−Fc
|
||||
</Box>
|
||||
</Tooltip>)}
|
||||
</Box>
|
||||
<Box sx={{position: 'relative', flexGrow: 1, minHeight: 0}}>
|
||||
<Box sx={{position: 'absolute', inset: 0, opacity: revealed ? 1 : 0,
|
||||
transition: 'opacity 900ms ease-out', '@media (prefers-reduced-motion: reduce)': {transition: 'none'}}}>
|
||||
<MolstarCanvas onReady={onPluginReady}/>
|
||||
</Box>
|
||||
{(modelLoading || (mapsLoading && !revealed)) && !error && <Box sx={{position: 'absolute', inset: 0, display: 'flex', alignItems: 'center',
|
||||
justifyContent: 'center', pointerEvents: 'none'}}>
|
||||
<Typography sx={{fontSize: '0.875rem', color: '#8aa0b5'}}>
|
||||
{modelLoading ? 'loading model…' : 'computing density maps…'}</Typography>
|
||||
</Box>}
|
||||
{error && <Alert severity="warning" sx={{position: 'absolute', top: 8, left: 8}}>{error}</Alert>}
|
||||
{!modelLoading && !mapsLoading && !hasMap.twofofc && !error &&
|
||||
<Typography sx={{position: 'absolute', bottom: 8, right: 12, fontSize: '0.8rem', color: '#8aa0b5'}}>No maps</Typography>}
|
||||
</Box>
|
||||
</Box>
|
||||
);
|
||||
}
|
||||
|
||||
@@ -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<string> {
|
||||
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<MapLayer> {
|
||||
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<void> {
|
||||
await plugin.build().delete(layer.dataRef).commit();
|
||||
}
|
||||
@@ -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;
|
||||
}
|
||||
|
||||
@@ -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);
|
||||
}
|
||||
@@ -105,6 +105,10 @@ std::vector<DensityBlob> FindDensityBlobs(const gemmi::Grid<float> &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)
|
||||
|
||||
@@ -2168,6 +2168,56 @@ std::vector<gemmi::Structure> SplitCopies(const gemmi::Structure &assembly, int
|
||||
return body;
|
||||
}
|
||||
|
||||
void CompactCopies(gemmi::Structure &placed, int copies) {
|
||||
const gemmi::UnitCell &cell = placed.cell;
|
||||
const std::vector<gemmi::Op> 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<double, 3> 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<gemmi::Op> op(copies, gemmi::Op::identity());
|
||||
std::vector<gemmi::Fractional> shift(copies);
|
||||
const std::vector<gemmi::Structure> 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<SubgroupSetting> 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<MRMerge> &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<gemmi::Structure>(*c.placed);
|
||||
CompactCopies(*compact, c.copies);
|
||||
c.placed = compact;
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
std::optional<MolecularReplacementResult> CarryMolecularReplacement(const MolecularReplacementResult &early,
|
||||
|
||||
@@ -263,6 +263,12 @@ gemmi::Structure AssembleCopies(const std::vector<gemmi::Structure> &copies);
|
||||
// The copies of an assembly (AssembleCopies), its chains split evenly in order.
|
||||
std::vector<gemmi::Structure> 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
|
||||
|
||||
@@ -196,6 +196,7 @@ ADD_EXECUTABLE(jfjoch_test
|
||||
RefinementTest.cpp
|
||||
RefineOnlyTest.cpp
|
||||
SolventSitesTest.cpp
|
||||
DensityBlobsTest.cpp
|
||||
SolventBuildTest.cpp
|
||||
SearchSpaceGroupTest.cpp
|
||||
SearchSpaceGroupTwinTest.cpp
|
||||
|
||||
@@ -0,0 +1,75 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <catch2/catch_all.hpp>
|
||||
|
||||
#include <cmath>
|
||||
#include <random>
|
||||
|
||||
#include <spdlog/fmt/fmt.h>
|
||||
#include <gemmi/blob.hpp>
|
||||
#include <gemmi/pdb.hpp>
|
||||
|
||||
#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<float> PeakMap(const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg, const gemmi::Fractional &f,
|
||||
std::mt19937 &rng) {
|
||||
gemmi::Grid<float> 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<float>(10.0 * std::exp(-r2 / 2.0));
|
||||
}
|
||||
grid.symmetrize_max();
|
||||
std::normal_distribution<float> 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<float> 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<DensityBlob> 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);
|
||||
}
|
||||
@@ -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<double, 3> 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<gemmi::Fractional> 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);
|
||||
|
||||
Reference in new issue
Block a user