Rimu.jl: Random integrators for many-body quantum systems

TL;DR

Rimu.jl presents a scalable, matrix-free framework for solving many-body quantum systems using FCIQMC and exact diagonalisation, implemented in Julia with a focus on Fock-state based representations. The package provides a unified interface via CommonSolve.jl, supports excited states, time series observables, and general operator expectations, and includes extensive data structures for Fock addresses, operators, and parallel vector handling. It introduces stochastic and semistochastic variants of the FCIQMC algorithm, discusses estimator techniques, and demonstrates performance on multi-core and MPI distributed platforms with detailed benchmarks. The work enables solving significantly larger problems than conventional dense or sparse matrix methods, by avoiding explicit matrix storage and leveraging matrix-free operator evaluations and replica-based observables. Overall, Rimu.jl advances high-performance simulation of quantum many-body systems, including ultra-cold atoms and lattice models, with practical tools for accuracy estimation, scalability, and extensibility.

Abstract

Rimu.jl is a Julia package for solving many-body quantum problems. The core of the package is a matrix-free implementation of Hamiltonians and other operators and compact representation of Fock states, which together allow for efficient methods suitable for high-performance computing. Rimu.jl includes a Julia implementation of the full configuration interaction quantum Monte Carlo (FCIQMC) algorithm which is a type of projector QMC algorithm for stochastically solving the time-independent Schrödinger equation. It also includes many well-known model Hamiltonians, and an interface for exact diagonalisation based on external eigenvalue solvers. Both the stochastic and exact diagonalisation methods are accessed with a CommonSolve.jl interface. We describe the FCIQMC algorithm and how to obtain estimators of observables as well as the key features of the implementation.

Figures (9)

• Figure 1: FCIQMC time series for the data obtained in the example. In panels (a) and (b), the shift and norm time series are shown. A close-up of the equilibrated region of the shift time series is shown in panel (c), which also shows the exact energy $E_0$ (dashed red line) and the estimate given by the mean of the shift $\bar{S}$ (orange line) with a blocked standard error $\sigma_{\bar{S}}$ (orange band).
• Figure 2: Estimates of the density-density correlator expectation values $\hat{G}^{(2)}(d)$ over a range of distances $d$. In panel (a), the orange bars show the estimates obtained from FCIQMC, where their widths are equal to 3 times the error bar. The exact values are plotted as dotted red lines. Panel (b) shows the error calculated from the differences between the estimates and the exact values.
• Figure 3: Stochastic vector compression. The plot shows the values $c_i$ at a given index $i$ of an uncompressed vector (yellow line) with a stochastically compressed vector (blue fill) superimposed on it. In the compressed vector, the coefficients with value $c_i \ge 1$ are unaffected by compression, while those below 1 are stochastically set to either 0 or 1 with probability proportional to $c_i$.
• Figure 4: The $M$-test performed on points of synthetic data. Panel a) shows how the standard error $\sigma_{\bar{X}}$ changes as the number of blocking steps $k$ increases. Note that past $k=8$, the values of $\sigma_{\bar{X}}$ are equivalent up to their error bars, hence the $M$-test identifies $k=8$ as the appropriate level of blocking. The bottom panel b) shows the raw data with the calculated mean and blocked standard error superimposed on it.
• Figure 5: A diagram showing a parallel (multi-threaded) operator-vector product $\check{\vb{T}}\vb{c}^{(n)}$. Segmentation of the working memory allows for fully parallel processing without race conditions as each core operates on its own memory segment. The operation is performed in two phases. First, the operator $\check{\vb{T}}$ is applied to each segment of $\vb{c}^{(n)}$ with the result stored column-wise in the working memory $W$. Then, the segments in $W$ are summed row-wise, (optionally) compressed, and stored to the result vector $\vb{c}^{(n+1)}$. The diagram shows a vector split into four segments, which would allow both phases to be performed in parallel on four CPU threads.