Computational backends#

The most time-consuming step in a dynasor calculation is the evaluation of the Fourier-transformed density

\[\rho(\boldsymbol{q}) = \sum_{i=1}^{N} e^{i \boldsymbol{q} \cdot \boldsymbol{r}_i}\]

and the corresponding current \(\boldsymbol{j}(\boldsymbol{q})\) for every frame in the trajectory. dynasor exposes three interchangeable backends for this kernel, selectable via set_compute_config, which applies to every subsequent call to compute_dynamic_structure_factors and compute_static_structure_factors:

`backend`	Library	Device	Notes
`'numba'`	Numba	CPU	Default. No extra dependencies.
`'torch'`	PyTorch	CUDA GPU	`pip install dynasor[gpu-torch]`. Fastest overall.
`'cupy'`	CuPy	CUDA GPU	See Installation: the CUDA-version-specific package is strongly recommended over the plain `cupy` extra. Near-identical speed to `'torch'`.

Both GPU backends are optional: importing dynasor never requires them, and they are loaded only when backend='torch' or backend='cupy' is explicitly configured.

Selecting a backend#

Call set_compute_config once, before computing anything; it takes effect for every subsequent call to compute_dynamic_structure_factors and compute_static_structure_factors in the process (there is no per-call backend/precision argument on either function):

from dynasor import Trajectory, compute_dynamic_structure_factors, set_compute_config

traj = Trajectory("trajectory.nc", trajectory_format="nc", ...)

# Default: Numba on CPU, no call to set_compute_config needed
sample = compute_dynamic_structure_factors(traj, q_points, dt=25.0,
                                           window_size=2000)

# PyTorch on GPU (CUDA)
set_compute_config(backend='torch')
sample = compute_dynamic_structure_factors(traj, q_points, dt=25.0,
                                           window_size=2000)

# CuPy on GPU (CUDA)
set_compute_config(backend='cupy')
sample = compute_dynamic_structure_factors(traj, q_points, dt=25.0,
                                           window_size=2000)

Use get_compute_config to inspect the current setting. Since the configuration is global, mixing backends within the same script means calling set_compute_config again before each call that needs a different one. There is intentionally only one way to select a backend, rather than a global default plus a per-call override.

The command-line interface exposes the same choice via --backend (numba/torch/cupy, default numba) and --precision (float64/float32, default float64):

dynasor -f trajectory.nc --trajectory-format nc --backend torch --precision float32 ...

dynasor.set_compute_config(backend=None, precision=None, batch_size=None)[source]#

Set the computational backend, precision, and/or GPU batch size used by all subsequent calls to compute_dynamic_structure_factors/compute_static_structure_factors.

Parameters:

backend (Optional[str]) – One of 'numba' (CPU, default), 'torch', or 'cupy' (both CUDA GPU, optional dependencies). Leave unset to keep the current value.
precision (Optional[str]) – One of 'float64' (default) or 'float32'. Only affects the GPU backends. Leave unset to keep the current value.
batch_size (Optional[int]) – Number of trajectory frames stacked into a single GPU kernel call (default 32). Only affects the GPU backends. Leave unset to keep the current value; see ComputeConfig.batch_size for tuning guidance.

Return type:

None

Examples

>>> from dynasor import set_compute_config
>>> set_compute_config(backend='torch', precision='float32')

dynasor.get_compute_config()[source]#

Return the current compute configuration.

Returns:: The current backend/precision/batch_size setting; ComputeConfig(backend='numba', precision='float64', batch_size=32) if set_compute_config() has never been called.
Return type:: ComputeConfig

The backend parameter is ignored for compute_spectral_energy_density, which always uses its own Numba kernel; it has no backend parameter at all. Mode projection (ModeProjector/project_modes()) is CPU-only in the same way, with no backend parameter and no GPU code path.

Incoherent scattering (calculate_incoherent=True) is a partial exception: even when backend='torch' or backend='cupy' is requested, the incoherent part of the calculation always runs on the CPU Numba kernel internally; only the coherent part (and, if requested, the currents) actually run on the GPU. This is a deliberate scope limitation, not a bug: the incoherent path recomputes a full density transform for every time lag in every window (an \(O(\text{window\_size})\) cost per window, rather than once per frame like the coherent path), and porting it to GPU is a larger, separate piece of work. The result is still correct: this only affects which device does the work, not the numbers, but a run with calculate_incoherent=True will not see the GPU speedup on that part of the calculation.

Floating-point precision#

GPU backends support an optional precision parameter:

`precision`	dtype	Notes
`'float64'`	double	Default. Full IEEE-754 double precision throughout.
`'float32'`	single	Several times faster on GPU, see Performance. Relative error on \(\rho(\boldsymbol{q})\) depends strongly on the magnitude of the atomic coordinates; see Accuracy below.

# Moves half the memory and uses the faster single-precision units;
# check the accuracy caveat below
set_compute_config(backend='torch', precision='float32')
sample = compute_dynamic_structure_factors(traj, q_points, dt=25.0,
                                           window_size=2000)

precision is silently ignored when backend='numba' (which always operates in float64).

Performance#

The kernel timings quoted here were measured on an NVIDIA GeForce RTX 3080 Ti with 20 000 atoms, 200 \(\boldsymbol{q}\)-points and a batch of eight frames, against the Numba backend on the same machine. Treat them as indicative: every ratio below depends on the GPU, the CPU it is compared against, the number of atoms, and the number of \(\boldsymbol{q}\)-points. Measure your own case before drawing conclusions from a factor.

The end-to-end speedup is smaller than the kernel speedup because trajectory I/O and the correlation-function accumulation (which run on CPU) are unaffected by the GPU backend. For larger systems or more \(\boldsymbol{q}\)-points the kernel fraction grows and the GPU advantage increases.

The 'torch' backend outperforms 'cupy' at the same precision because PyTorch’s reduction kernel (summing \(N_\text{atoms}\) contributions per \(\boldsymbol{q}\)-point) is significantly faster for the typical \((N_\text{atoms}, N_q)\) shapes encountered in practice.

Memory usage and batch size#

GPU backends process frames in fixed-size batches (32 frames per batch by default) so that each kernel call is large enough to make good use of GPU memory bandwidth. A background thread reads the next batch from disk while the GPU processes the current one. The batch size is configurable via set_compute_config:

set_compute_config(backend='torch', batch_size=8)

Memory usage scales with \(\text{batch\_size} \times N_\text{atoms} \times N_q\), and it is not adaptive: very large systems or q-point grids can exhaust GPU memory (an out-of-memory error from PyTorch/CuPy, not a graceful fallback). If you hit this, reduce batch_size (the default of 32 is tuned for small-to-medium systems; large systems may need single digits, or even batch_size=1), reduce \(N_q\) (e.g. fewer q-points per call, run in multiple passes), or switch to precision='float32' (halves the memory footprint). Conversely, for small systems where disk I/O rather than GPU compute is the bottleneck, a larger batch_size than the default may improve throughput by giving the GPU more work per kernel call relative to the fixed per-call overhead.

Accuracy#

precision='float64' is indistinguishable from the Numba reference to machine precision (\({\sim}10^{-12}\) relative error), regardless of the magnitude of the atomic coordinates.

precision='float32' accuracy is dominated by rounding in the phase argument \(\boldsymbol{q} \cdot \boldsymbol{r}_i\), not by the summation over atoms that follows it. dynasor never wraps or re-references atomic positions: there is no periodic-boundary wrapping anywhere in the trajectory or core code, and positions are used exactly as read from the trajectory file. As a result, \(\boldsymbol{r}_i\) can be far larger than the simulation cell for long or diffusive trajectories with unwrapped coordinates: hundreds to thousands of Ångström is common. Since float32 has a fixed relative precision (\(\varepsilon_{32} \approx 1.2\times10^{-7}\)), the absolute rounding error on a coordinate grows with its magnitude, and so does the resulting phase error and the relative error on \(\rho(\boldsymbol{q})\), roughly

\[\text{relative error} \sim q_\text{max} \, |\boldsymbol{r}|_\text{max} \, \varepsilon_{32}\]

For small, box-local coordinates (\(|\boldsymbol{r}| \lesssim 10\) Å) this stays around \(10^{-6}\) to \(10^{-5}\). For unwrapped coordinates reaching \(10^3\) to \(10^4\) Å, which is realistic for long diffusive trajectories, the same rule of thumb gives \(10^{-3}\) to \(10^{-2}\) at typical \(q_\text{max}\), and this has been confirmed empirically: with \(q_\text{max}=40\) rad/Å and \(|\boldsymbol{r}|=10^4\) Å the measured relative error on \(\rho(\boldsymbol{q})\) is \({\sim}1.5\times10^{-2}\), roughly 1000 times worse than the small-coordinate case above (see tests/test_gpu.py’s test_float32_relative_error_grows_with_unwrapped_position_magnitude). Structure factors \(S(\boldsymbol{q},\omega)\) derived from ratios of such quantities inherit a similar relative error. This is specific to precision='float32': the float64 path stays accurate regardless of coordinate magnitude, since \(\varepsilon_{64} \approx 2\times 10^{-16}\) keeps the phase error negligible even at \(|\boldsymbol{r}| \sim 10^4\).

If your trajectory has large unwrapped displacements (diffusive liquids, long simulations, vacancy migration, …) and you are considering precision='float32', check the accuracy on a few frames against precision='float64' first. The \({\sim}10^{-5}\) figure often quoted for this feature assumes small, box-local coordinates and does not apply once atomic coordinates grow large.

When to use which backend#

Use the following decision guide (all via set_compute_config, see above):

No GPU available: use the 'numba' default; no call to set_compute_config needed. Numba automatically prefers TBB over OpenMP over its own workqueue fallback for its parallel CPU kernels; if TBB is not installed (pip install tbb), dynasor logs a one-time warning if it ends up on the slower workqueue layer.
GPU available, full double precision required: use set_compute_config(backend='torch') or set_compute_config(backend='cupy'). The two land close to each other, with 'torch' slightly ahead (kernel speedups of 12× and 11× under the conditions above). If PyTorch is already installed, prefer 'torch'; otherwise 'cupy' avoids the large PyTorch download.
GPU available, coordinates stay small (box-local, not unwrapped): use set_compute_config(backend='torch', precision='float32'). This is the fastest option, roughly 4× the float64 kernel throughput, and is appropriate for most analyses, with ~10⁻⁵ relative error. See Accuracy: this figure does not hold for large unwrapped/diffusive coordinates.
Very large systems (\(N_\text{atoms} \gtrsim 10^5\)): the kernel dominates and the end-to-end speedup approaches the kernel-level ratio, an order of magnitude or more at float64 on the hardware above. precision='float32' gives maximum throughput, but check the Accuracy caveat first if the trajectory has large unwrapped displacements. Large systems and large diffusive displacements often go together (e.g. long liquid-state simulations).

Installation#

The 'torch' and 'cupy' backends are not installed by default; both are optional dependencies.

PyTorch:

pip install dynasor[gpu-torch]   # or: pip install torch, following
                                 # https://pytorch.org/get-started

CuPy is more subtle, and there are two ways to install it:

pip install dynasor[gpu-cupy]   # plain "cupy" from PyPI
pip install cupy-cuda12x        # or cupy-cuda11x, etc. (recommended instead)

pip install dynasor[gpu-cupy] pulls in the plain cupy package from PyPI, which compiles its CUDA extensions from source at install time. This requires a local CUDA toolkit and a C++ compiler to already be present, and the install itself is slow and more prone to failing if your environment doesn’t have exactly what CuPy’s build expects.

The cupy-cudaXXx packages (e.g. cupy-cuda11x, cupy-cuda12x), by contrast, are prebuilt binary wheels with the matching CUDA runtime libraries already bundled in, so no compiler or local CUDA toolkit is needed, and installation is fast and reliable. If you know which CUDA version your system has, installing the matching cupy-cudaXXx package directly (instead of the generic extra) is the recommended path; see CuPy’s own installation guide for which package matches your CUDA version.

An additional optional package for the Numba CPU backend is TBB, which Numba automatically prefers over OpenMP and over its own workqueue fallback for its parallel CPU kernels; no code changes or environment variables are needed:

pip install tbb

If neither TBB nor OpenMP is available, Numba falls back to the slower workqueue layer; dynasor logs a one-time warning the first time this happens, recommending pip install tbb.

Multi-GPU#

By default both GPU backends use the first available CUDA device (cuda:0). To target a different device, set the relevant environment variable before importing dynasor:

CUDA_VISIBLE_DEVICES=1 python my_script.py

or call the library’s device-selection API before the computation:

# PyTorch
import torch
torch.cuda.set_device(1)

# CuPy
import cupy as cp
cp.cuda.Device(1).use()

Either call selects the device for the thread that runs it, so make the call in the same thread that then calls the compute function, and before that call. The device is picked up when the q-points are transferred to the GPU, at the start of the computation, and the whole calculation stays on it: the per-frame reading and the GPU kernels run on a background thread, which starts out with the first device current, but the kernels take their device from the q-point array rather than from the thread they run in.

A single calculation always runs on a single device. Neither backend splits one calculation across several GPUs.