leonarski_fandClaude Opus 5.5
8d17d05227
Where the molecular transform's FFT runs on the CPU - a CPU-only build, --no-gpu, a card too small for it - it was one single-threaded FFTW 3D r2c: 12.5 s of 13 s of the transform on 4nwv's 768^3 box, one core busy while fifteen waited. Now MapToFPhiInParallel does it as two passes over the threads: the r2c along w of every (u, v) column, then the 2D transform of each plane l the search keeps (l <= reach), cropping as CropHalfL does. Each column and each plane is done by one plan, so the result does not depend on the thread count; it agrees with FFTW's 3D plan to float rounding, not bit for bit. The GPU path is unchanged. CPU-only (--no-gpu), timing benchmark: 4nwv's transform 13.1 -> 2.4 s (provisional, load average 13-25). What the rounding moved: 4nwv and myob came out byte-identical; on cytc (three copies, tNCS) the null's replicates, placed by the translation function, read the fit at +22.8 sd instead of +18.7 (MR_HAND_FIT_SIGMA), and the placement - same R-work, R-free and every MR candidate line - came out at an origin-equivalent position 43 A from the previous one, so its model file and maps differ while the answer (ENANTIOMORPH, P 32 2 1, placement CC 0.744) does not. GPU runs byte-identical on all three. [MolecularReplacement] passes. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
42 lines
2.8 KiB
C++
42 lines
2.8 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#pragma once
|
|
|
|
#include "gemmi/grid.hpp"
|
|
#include "gemmi/recgrid.hpp"
|
|
|
|
// The structure factors of a real-space map on the model path (Fcalc from the model density, Fmask
|
|
// from the bulk-solvent mask): what gemmi::transform_map_to_f_phi(map, true) returns - the
|
|
// conjugated forward transform (F = sum rho exp(+2 pi i h.x)) scaled by V/N, kept for l >= 0 on an (nu, nv, nw/2+1) grid in the map's own XYZ order,
|
|
// ready for prepare_asu_data() - computed with FFTW rather than gemmi's single-threaded pocketfft.
|
|
// The same transform, not the same bits: the two libraries order their arithmetic differently.
|
|
//
|
|
// Plans are FFTW_ESTIMATE only - MEASURE would time candidate algorithms and could pick a different
|
|
// one from run to run - and made once per grid size and kept, so the rigid body's hundred-odd
|
|
// evaluations on one grid plan once. Safe to call from several threads at once.
|
|
gemmi::FPhiGrid<float> MapToFPhi(const gemmi::Grid<float> &map);
|
|
|
|
// The inverse: the real-space map of an (nu, nv, nw/2+1) half-l coefficient grid in XYZ order, as
|
|
// get_f_phi_on_grid() fills it - what gemmi::transform_f_phi_grid_to_map() returns (rho = 1/V sum F
|
|
// exp(-2 pi i h.x) on an (nu, nv, 2*(nw-1)) grid, NaN coefficients read as zero), computed with FFTW
|
|
// under the same planning rules.
|
|
gemmi::Grid<float> MapFromFPhi(const gemmi::FPhiGrid<float> &hkl);
|
|
|
|
// MapFromFPhi's transform itself, for a caller that fills its input: `in`, the (nu, nv, nw/2+1) half-l coefficients
|
|
// u fastest as MapFromFPhi hands them to FFTW (conjugated, a NaN as zero), is overwritten, and `out` receives the
|
|
// (nu, nv, nw) map before MapFromFPhi's 1/V. Both from fftwf_malloc, as a plan requires.
|
|
void MapFromFPhiExecute(int nu, int nv, int nw, std::complex<float> *in, float *out);
|
|
|
|
// CropHalfL(MapToFPhi(map), reach) on `nthreads` threads, for a map of hundreds of millions of points (the
|
|
// molecular transform's box), whose single-threaded transform is minutes of CPU: the r2c along w of every
|
|
// (u, v) column, then the 2D transform of every plane l <= reach - each column and each plane by one plan, so
|
|
// the result does not depend on the thread count. The same transform as MapToFPhi to float rounding, not bit
|
|
// for bit: FFTW's 3D plan orders its arithmetic its own way.
|
|
gemmi::FPhiGrid<float> MapToFPhiInParallel(const gemmi::Grid<float> &map, int reach, size_t nthreads);
|
|
|
|
// The half-l grid `full` with only the indices |k|, l <= reach kept, on a grid of full.nu x (2 reach + 1) x
|
|
// (reach + 1) that an index wraps into as it wraps into the whole one (the molecular transform's, which is
|
|
// read no further out).
|
|
gemmi::FPhiGrid<float> CropHalfL(const gemmi::FPhiGrid<float> &full, int reach);
|