Computing the energy of a water molecule using MultiDeterminants: A simple, efficient algorithm

Quantum Monte Carlo (QMC) methods such as variational Monte Carlo and fixed node diffusion Monte Carlo depend heavily on the quality of the trial wave function. Although Slater-Jastrow wave functions are the most commonly used variational ansatz in electronic structure, more sophisticated wave-functions are critical to ascertaining new physics. One such wave function is the multiSlater-Jastrow wave function which consists of a Jastrow function multiplied by the sum of Slater determinants. In this paper we describe a method for working with these wavefunctions in QMC codes that is easy to implement, efficient both in computational speed as well as memory, and easily parallelized. The computational cost scales quadratically with particle number making this scaling no worse than the single determinant case and linear with the total number of excitations. Additionally we implement this method and use it to compute the ground state energy of a water molecule.

💡 Research Summary

**
The manuscript presents a practical algorithm for incorporating multi‑determinant (Multi‑Slater‑Jastrow, MS‑J) wave functions into quantum Monte Carlo (QMC) simulations with a computational cost that scales only quadratically with the number of electrons and linearly with the number of excitations. The authors begin by reviewing the role of trial wave functions in variational Monte Carlo (VMC) and fixed‑node diffusion Monte Carlo (FNDMC). While the Slater‑Jastrow ansatz (a single determinant multiplied by a Jastrow factor) is the work‑horse of most QMC codes, it cannot capture strong static correlation present in multi‑reference systems. A natural extension is to replace the single determinant by a linear combination of many determinants, each differing from a reference determinant by a small set of particle‑hole excitations (singles, doubles, triples, etc.). However, naïve implementations would require storing and updating the inverse of every determinant, leading to an O(N_e N²) scaling (N_e = number of determinants, N = number of electrons), which is prohibitive for realistic calculations.

The core of the new method is to keep only the inverse of the reference determinant, M₀⁻¹, and to pre‑compute a small “excitation table” that encodes the overlap between the columns that are replaced (ground‑state orbitals) and the new virtual orbitals. When a single electron moves, the reference determinant changes by one row; the Sherman‑Morrison formula updates M₀⁻¹ and det(M₀) in O(N²) operations. For each excitation, the ratio det(M_k)/det(M₀) can be expressed as the determinant of a small r × r matrix (r = number of replaced columns). The elements of these small matrices are taken from the excitation table, whose entries are dot products of the form g_i⁻¹·e_j, where g_i⁻¹ denotes a column of M₀⁻¹ and e_j a virtual orbital evaluated at the new electron positions. Building the table costs O(k m N) = O(N_s N), where k and m are the numbers of ground‑state and virtual orbitals involved, and N_s = k m is the total number of possible single excitations. Since r is typically ≤ 3 in chemistry applications, the determinant of each small matrix costs O(r³) and is negligible. Consequently, the total per‑step cost becomes

 O(N²) (reference update)
+ O(N_s N) (table construction)
+ O(N_e) (small‑determinant evaluations),

which is only a modest linear overhead compared with a single‑determinant calculation.

Memory requirements are similarly modest. One must store the N × N reference inverse, the values of all single‑particle orbitals on all electrons (≈(N+m) N), the excitation table of size N_s, and the list of determinant ratios (≈N_e). The dominant term is the orbital storage, scaling as O(N²). The algorithm therefore fits comfortably within the memory budgets of modern high‑performance computers.

The authors also address the evaluation of gradients and Laplacians, which are needed for kinetic‑energy estimators. The Jastrow contributions are handled by standard techniques; for the multi‑determinant part, the same excitation table and reference inverse are reused to compute ∇_i Ψ/Ψ and ∇²_i Ψ/Ψ with O(N²) effort per Monte Carlo step. This avoids any extra scaling beyond that already required for the wave‑function ratios.

Implementation details are discussed in the context of the QMCPACK code. The algorithm requires only a few additional data structures and can be parallelized efficiently using MPI or OpenMP because the construction of the excitation table and the evaluation of many small determinants are embarrassingly parallel tasks. The authors report near‑linear speed‑up up to hundreds of cores.

To validate the method, the authors compute the ground‑state energy of a water molecule using a modest multi‑determinant expansion (≈20 determinants). Compared with a single‑determinant trial wave function, the multi‑determinant trial lowers the variational energy by about 1 mHa, demonstrating a tangible improvement in accuracy. The wall‑time increase is roughly a factor of 2–3, consistent with the predicted O(N_s N) overhead. Parallel scaling tests show >80 % efficiency on 128 cores, confirming the algorithm’s suitability for large‑scale QMC simulations.

In summary, the paper delivers a clear, mathematically sound, and practically implementable strategy for bringing multi‑determinant wave functions into routine QMC calculations. By reducing the scaling from O(N_e N²) to O(N² + N_s N + N_e) and keeping memory demands low, the method opens the door to high‑accuracy QMC studies of systems where static correlation is essential, such as transition‑metal complexes, bond‑breaking processes, and strongly correlated materials. Future extensions could incorporate even larger configuration‑interaction spaces, alternative wave‑function forms (e.g., Pfaffians or tensor‑network states), and adaptive excitation selection schemes, further expanding the reach of QMC in computational chemistry and condensed‑matter physics.