Quantum Monte Carlo (QMC) simulations for the recent studies on complex materials were confronted by new computational challenges. Traditional approach to accelerate the simulations by parallel Monte Carlo chains faces serious scalability problems since the speedup is reaching the limitation predicted by Amdahl's law. Fine-grained parallelization of matrix kernels is essential to achieve better performance. In this paper, we investigate the performance optimization techniques on GPU for the most time consuming computational kernel in the Determinant Quantum Monte Carlo (DQMC) simulation: multiplication of matrix exponentials. The matrix, derived from the kinetic Hamiltonian, is highly sparse, and its exponential is approximated by the block checkerboard method, which can represent a matrix exponential as a product of a sequence of sparse matrices. The matrix exponentials from 2D and 3D toruses are focused, and various optimization techniques, such as data streaming and concurrent kernels, are proposed. Experiments show that the proposed optimization techniques can improve the SpMM (Sparse Matrix Multiplication) function, modified from the CUDA SKD SpMV function, up to 16 times and 117 times for 2D and 3D problems respectively.
Quantum Monte Carlo (QMC) simulations for the recent studies on complex materials were confronted by new computational challenges. Traditional approach to accelerate the simulations by parallel Monte Carlo chains faces serious scalability problems since the speedup is reaching the limitation predicted by Amdahl's law. Fine-grained parallelization of matrix kernels is essential to achieve better performance. In this paper, we investigate the performance optimization techniques on GPU for the most time consuming computational kernel in the Determinant Quantum Monte Carlo (DQMC) simulation: multiplication of matrix exponentials. The matrix, derived from the kinetic Hamiltonian, is highly sparse, and its exponential is approximated by the block checkerboard method, which can represent a matrix exponential as a product of a sequence of sparse matrices. The matrix exponentials from 2D and 3D toruses are focused, and various optimization techniques, such as data streaming and concurrent kernels, are proposed. Experiments show that the proposed optimization techniques can improve the SpMM (Sparse Matrix Multiplication) function, modified from the CUDA SKD SpMV function, up to 16 times and 117 times for 2D and 3D problems respectively.
Determinant Quantum Monte Carlo (DQMC) simulation is one of few numerical methods that can explore the micro properties of fermions, which has many technically important applications in chemistry and material science. Conventionally, its parallelization relies on parallel Monte Carlo method, whose speedup is limited by the thermalization process and the underlying matrix computation. To achieve better performance, fine-grained parallelization on its numerical kernel is essential to utilize the massive parallel processing units, which are multicores and/or GPUs interconnected by high performance network. In this paper, we address the parallelization on one of the matrix kernel in the DQMC simulations: the multiplication of matrix exponentials. The matrix is derived from the kinetic Hamiltonian, which is highly sparse. We approximate its exponential by the checkerboard method, which decomposes the matrix exponential into a product of a sequence of block sparse matrices. We analyze the block sparse matrices of two common used lattice geometry: 2D torus and 3D cubic, and parallelize the computational kernel of multiplying them to a general matrix. The parallel algorithm is designed for multicore CPU and GPU. The results of experiments showed on a quad core processor, 3 times speedup can be observed in average, and on GPU, 145 times speedup is achievable.
QMCPACK is open source scientific software designed to perform Quantum Monte Carlo (QMC) simulation, a first-principles method for describing many-body systems. The evaluation of each Monte Carlo move requires finding the determinant of a dense matrix in the wave functions. This calculation forms a key computational kernel in QMCPACK. After each accepted event, the wave function matrix undergoes a rank-one update to represent a single particle move within the system. The matrix inverse is updated via the Sherman-Morrison formula; occasionally, the explicit inverse must be recomputed to maintain numerical stability.; We propose an alternate approach to this kernel that utilizes QR factorization to maintain stability without refactorization. In addition, algorithms based on a novel delayed update scheme are explored in this effort. This strategy enables probability evaluation for multiple successive Monte Carlo moves, with application of accepted moves to the wave function matrix delayed until an event is denied, or a predetermined limit on acceptances p is reached. Accepted events grouped in this manner are then applied to the matrix en bloc with enhanced arithmetic intensity and computational efficiency.
Quantum Monte Carlo (QMC) methods deliver highly accurate electronic structure calculations but are computationally intensive. The quantum Monte Carlo kernel library (QMCkl) provides a modular, portable collection of high-performance kernels implementing the core building blocks of QMC calculations. It offers a C-compatible application programming interface, supports the TREXIO standard for input, and covers essential QMC kernels including atomic and molecular orbitals, cusp corrections, the Jastrow factor, and the necessary derivatives also to perform variational and structural optimization. QMCkl separates algorithmic development from hardware-specific tuning by combining human-readable reference implementations with performance-optimized kernels that produce identical numerical results. The library enables consistent, efficient, and reproducible simulations across different QMC codes and architectures and achieves substantial speedups in the evaluation of the energy and its derivatives. Beyond QMC, QMCkl can accelerate deterministic quantum chemistry workflows and visualization tools, promoting cross-code interoperability and simplifying high-performance scientific software development.
We present jQMC, a Python-based computational package for {\it ab initio} Quantum Monte Carlo (QMC) simulations, designed for modern GPU-accelerated computing environments. jQMC implements two well-established QMC algorithms: Variational Monte Carlo (VMC) and the lattice-regularized variant of Diffusion Monte Carlo (LRDMC). The employed wave function is a Jastrow factor combined with the antisymmetrized geminal power with spin-singlet and spin-triplet pairings, which contains the single Slater determinant as its special lowest-rank case. The wave function can be initialized from external Hartree-Fock/Density Functional Theory calculations through the TREX-IO library (a common wave-function format across electronic-structure packages) and optimized by stochastic reconfiguration and linear-method energy minimization. One of the prominent features of jQMC is its use of JAX, which enables automatic differentiation for wave function optimization and atomic force calculations, and allows the main QMC algorithms to be Just-In-Time (JIT) compiled and portable across CPU and GPU. jQMC is vectorized over walkers at the top level of the QMC algorithms, providing efficient intra-GPU~(CPU) vectorization. The multi-GPU~(CPU) parallelization is also supported through MPI and JAX sharding. To assess the practical performance of this implementation, we benchmarked jQMC performance on NVIDIA GPUs (A100 and H100) and analyzed CUDA kernels. For the test cases analyzed here, with system sizes up to 160 electrons, the current version of jQMC is faster than TurboRVB, a Fortran90 code implementing the same algorithms and wave functions, once jQMC is run on GPUs. In terms of wall-time, the gain can reach an order of magnitude for VMC, while it is more moderate for LRDMC.
Determinant Quantum Monte Carlo (DQMC) provides numerically exact solutions for strongly correlated fermionic systems but faces significant computational challenges with increasing system size. While submatrix updates were originally developed for Hirsch-Fye QMC with onsite interactions at finite temperatures, their comprehensive application in DQMC has remained unexplored despite noted algorithmic similarities. We present the first comprehensive application of submatrix updates in DQMC, significantly extending beyond the original scope by enabling simulations with extended interactions and at zero temperature. Building upon conventional fast updates and delay updates, our generalized implementation achieves an order-of-magnitude improvement in computational efficiency, enabling simulations of the half-filled Hubbard model on lattices up to 8,000 sites - a scale previously challenging with standard DQMC implementations. This enhanced computational capability allows us to accurately determine the finite-temperature phase diagram of the 3D Hubbard model at half-filling. Our findings not only shed light on the phase transitions within these complex systems but also pave the way for more effective simulations of strongly correlated electrons, potentially guiding experimental efforts in cold atom simulations of the 3D Hubbard model.
Density matrix renormalization group (DMRG), as a numerical method of solving the ground state of one-dimensional strongly-correlated lattice model with very high accuracy, requires expensive computational and memory cost when applied to two- and quasi-two-dimensional problems. The number of DMRG kept states is generally very large to achieve a reliable accuracy for these applications, which results in numerous matrix and vector operations and unbearably consuming time in the absence of the proper parallelization. However, due to its sequential nature, the parallelization of DMRG algorithm is usually not straightforward. In this work, we propose a new hybrid parallelization strategy for the DMRG method. It takes advantage of the computing capability of both central processing unit (CPU) and graphics processing unit (GPU) of the computer. In order to achieve as many as DMRG kept states within a limited GPU memory, we adopt the four-block formulation of the Hamiltonian rather than the two-block formulation. The later consumes much more memories, which has been used in another pioneer work on the hybrid parallelization of the DMRG algorithm, and only a small number of DMRG kept states are available. Our parallel strategy focuses on the diagonalization of the Hamiltonian, which is the most time-consuming part of the whole DMRG procedure. A hybrid parallelization strategy of diagonalization method is implemented, in which the required data for diagonalization are distributed on both the host and GPU memory, and the data exchange between them is negligible in our data partitioning scheme. The matrix operations are also shared on both CPU and GPU when the Hamiltonian acts on a wave function, while the distribution of these operations is determined by a load balancing strategy. Taking fermionic Hubbard model for example, we examine the running performance of the hybrid parallelization strategy with different DMRG kept states and provide corresponding performance benchmark. On a 4-leg ladder, we employ the conserved quantities with U(1) symmetry of the model and a good-quantum number based task scheduling to further reduce the GPU memory cost. We manage to obtain a moderate speedup of the hybrid parallelization for a wide range of DMRG kept states. In our example, the ground state energy with high accuracy is obtained by the extrapolation of the results, with different numbers of states kept, and we show charge stripes which are usually experimentally observed in high-temperature superconductors. In this case, we keep 104 DMRG states and the GPU memory cost is less than 12 Gigabytes.