Challenges
Datasets
Workspaces
Papers
Discussions
Leaderboard
Log inSign up
Challenges
Datasets
Workspaces
Papers
Discussions
Leaderboard
Blog
Job Board
Q3AS

© 2026 Aqora Quantum S.A.S.

TermsPrivacyLegal Notice

arXiv:2506.12501 · The Journal of Physical Chemistry A

Implementation of McMurchie–Davidson Algorithm for Gaussian AO Integrals Suited for SIMD Processors

Andrey Asadchev, Edward F. Valeev

Published Oct 13, 2025• 1 citation
arXivHTMLPDF

Abstract

High Resolution Image Download MS PowerPoint Slide We report an implementation of the McMurchie–Davidson evaluation scheme for 1- and 2-particle Gaussian AO integrals designed for processors with Single Instruction Multiple Data (SIMD) instruction sets. Like in our recent MD implementation for graphical processing units (GPUs) [ Asadchev, A.; Valeev, E. F. . J. Chem. Phys. 2024, 160, 244109 .], variable-sized batches of shellsets of integrals are evaluated at a time. By optimizing for the floating point instruction throughput rather than minimizing the number of operations, this approach achieves up to 50% of the theoretical hardware peak FP64 performance for many common SIMD-equipped platforms (AVX2, AVX512, NEON), which translates to speedups of up to 30 over the state-of-the-art one-shellset-at-a-time implementation of Obara–Saika-type schemes in Libint for a variety of primitive and contracted integrals. As with our previous work, we rely on the standard C++ programming language─such as the std::simd standard library feature to be included in the 2026 ISO C++ standard─without any explicit code generation to keep the code base small and portable. The implementation is part of the open source LibintX library freely available at https://github.com/ValeevGroup/libintx .

DOI
10.1021/acs.jpca.5c04136

Related papers

3-center and 4-center 2-particle Gaussian AO integrals on modern accelerated processors

Andrey Asadchev, Edward F. Valeev

We report an implementation of the McMurchie-Davidson (MD) algorithm for 3-center and 4-center 2-particle integrals over Gaussian atomic orbitals (AOs) with low and high angular momenta l and varying degrees of contraction for graphical processing units (GPUs). This work builds upon our recent implementation of a matrix form of the MD algorithm that is efficient for GPU evaluation of 4-center 2-particle integrals over Gaussian AOs of high angular momenta (l ≥ 4) [A. Asadchev and E. F. Valeev, J. Phys. Chem. A 127, 10889-10895 (2023)]. The use of unconventional data layouts and three variants of the MD algorithm allow for the evaluation of integrals with double precision and sustained performance between 25% and 70% of the theoretical hardware peak. Performance assessment includes integrals over AOs with l ≤ 6 (a higher l is supported). Preliminary implementation of the Hartree-Fock exchange operator is presented and assessed for computations with up to a quadruple-zeta basis and more than 20 000 AOs. The corresponding C++ code is part of the experimental open-source LibintX library available at https://github.com/ValeevGroup/libintx.

arXiv:2405.01834•The Journal of Chemical Physics•Jun 27, 2024•15 citations

Optimizing two-electron repulsion integral calculations with McMurchie–Davidson method on graphic processing unit

Yingqi Tian, Bingbing Suo, Yingjin Ma +1 more

In this article, several optimization methods of two-electron repulsion integral calculations on a graphic processing unit (GPU) are presented. These methods are based on the investigations of the method presented by McMurchie and Davidson (MD). A new Boys function evaluation method for the GPU calculation is introduced. The series summation, the error function, and the finite sum formula method are combined; thus, good performance on the GPU can be achieved. By taking some theoretical study of the McMurchie-Davidson recurrence relations, three major optimization approaches are derived from the deduction of the general term formula for the Hermite expansion coefficient. The three approaches contain a new form of the Hermite expansion coefficients with corresponding recurrence relations, which is more efficient for one-electron integrals and [ss|∗∗] or [∗∗|ss] type two-electron integrals. In addition, a simple yet efficient new recurrence formula for the coefficient evaluation is derived, which is more efficient both in float operations and memory operations than its original one. In average, the new recurrence relation can save 26% float operations and 37% memory operations. Finally, a common sub-expression elimination (CSE) method is implemented. This CSE method is directly generated from some equalities we discovered from the general term formula other than by computer algebra system software. This optimized method achieved up to 3.09 speedups compared to the original MD method on the GPU and up to 92.75 speedups compared to the GAMESS calculation on the central processing unit.

The Journal of Chemical Physics•Jul 20, 2021•21 citations

Efficient GPU Implementation of the McMurchie–Davidson Method for Shell-Based ERI Computations

Haruto Fujii, Yasuaki Ito, Nobuya Yokogawa +5 more

Quantum chemistry offers the formal machinery to derive molecular and physical properties arising from (sub)atomic interactions. However, as molecules of practical interest are largely polyatomic, contemporary approximation schemes such as the Hartree–Fock scheme are computationally expensive due to the large number of electron repulsion integrals (ERIs). Central to the Hartree–Fock method is the efficient computation of ERIs over Gaussian functions (GTO-ERIs). Here, the well-known McMurchie–Davidson method (MD) offers an elegant formalism by incrementally extending Hermite Gaussian functions and auxiliary tabulated functions. Although the MD method offers a high degree of versatility to acceleration schemes through Graphics Processing Units (GPUs), the current GPU implementations limit the practical use of supported values of the azimuthal quantum number. In this paper, we propose a generalized framework capable of computing GTO-ERIs for arbitrary azimuthal quantum numbers, provided that the intermediate terms of the MD method can be stored. Our approach benefits from extending the MD recurrence relations through shells, batches, and triple-buffering of the shared memory, and ordering similar ERIs, thus enabling the effective parallelization and use of GPU resources. Furthermore, our approach proposes four GPU implementation schemes considering the suitable mappings between Gaussian basis and CUDA blocks and threads. Our computational experiments involving the GTO-ERI computations of molecules of interest on an NVIDIA A100 Tensor Core GPU (NVIDIA, Santa Clara, CA, USA) have revealed the merits of the proposed acceleration schemes in terms of computation time, including up to a 72× improvement over our previous GPU implementation and up to a 4500× speedup compared to a naive CPU implementation, highlighting the effectiveness of our method in accelerating ERI computations for both monatomic and polyatomic molecules. Our work has the potential to explore new parallelization schemes of distinct and complex computation paths involved in ERI computation.

Applied Sciences•Feb 27, 2025•6 citations

The updates in Libcint 6: More integrals, API refinements, and SIMD optimization techniques

Qiming Sun

Libcint is a library designed for the evaluation of analytical integrals for Gaussian type orbitals. It prioritizes simplicity, ease of use, and efficiency for the development of quantum chemistry programs. In the release of version 6.0, Libcint supports the computation of integrals for various operators, such as overlap, Coulomb, Gaunt, Breit, attenuated Coulomb, Slater-type geminals, and Yukawa potential, as well as arbitrary orders of derivatives for these operators. To enhance the usability of the library, Libcint provides a uniform function signature for all integral functions. A code generator is included to automate the implementation of new integrals. To achieve better performance on modern central processing unit architectures, the library employs explicit single instruction multiple data parallelization in the code implementation.

The Journal of Chemical Physics•May 6, 2024•15 citations

GPU ‐Accelerated One‐Electron Integral Computation for Quantum Chemistry

Nobuya Yokogawa, Yasuaki Ito, Satoki Tsuji +5 more

ABSTRACT In Quantum chemical computation, numerical schemes such as the Hartree–Fock (HF) and density functional theory (DFT) are widely used to solve the Schrödinger equation numerically, to realize experiment‐free prediction and analysis of key molecular properties such as structure and energy. Computing one‐electron integrals, such as kinetic energy integrals and nuclear attraction integrals, is essential in both HF and DFT to characterize the molecular electronic states. However, as molecules of practical interest grow in size and angular momentum, computing one‐electron orbitals becomes computationally expensive in most cases. Although computing kinetic energy integrals on CPUs is straightforward, bottlenecks in CPU‐GPU data transfer have often been overlooked. In this study, we propose an efficient method to compute both the kinetic‐energy and nuclear‐attractive integrals on GPUs. First, we explicitly and symbolically expand recurrence relations based on the Obara–Saika and McMurchie–Davidson methods to eliminate redundant operations, thus improving computational efficiency. Second, we implemented a hybrid method that selects the best/fastest of both methods depending on the integration task. Third, we achieved further speedups by using CUDA streams to parallelize the execution of multiple kernels and efficiently utilize multiprocessor resources on the GPU. Computational experiments using NVIDIA A100 GPUs and Intel Xeon Gold 6338 CPU on relevant molecules of interest demonstrated the superiority of our one‐electron integral GPU implementations, achieving a speedup of 20.2 times over PySCF, and a speedup of 132.6 times over GPU4PySCF.

Concurrency and Computation Practice and Experience•Feb 20, 2026•0 citations

Data-Driven Optimization of Seminumerical Exchange for SIMD CPUs

Jesse DiCenso, Henryk Laqua, Martin Head-Gordon

We present a seminumerical exchange algorithm (SNK) based on the sn-LinK method [Laqua et al., J. Chem. Theory Comput. 16, 1456 (2020)], optimized for modern multi-core central processing units with single instruction multiple data vector instructions and report its implementation into the Q-Chem quantum chemistry program. Vectorization and memory management strategies for the Boys function and 3-center-1-electron (3c1e) integral kernels are given. Several new considerations are made in the computation of the 3c1e integrals: the Head-Gordon Pople horizontal recurrence relations, same-center optimizations, and multishell optimizations, all implemented symbolically employing the SymPy symbolic algebra library and exported as optimized C++ code leveraging common-sub-expression elimination (CSE) optimization techniques. A general scheme for choosing optimal integral kernels is developed through rigorous timing benchmarks. The resulting implementation is tested against traditional 4c2e exchange for several systems and basis sets, showing 1-2 orders of magnitude speedups for systems of practically relevant sizes (100-500 atoms, 2000-10000 basis functions). For smaller molecules, comparisons are also made against the more efficient occ-RI-K method showing a break-even point around 50-100 atoms beyond which SNK dramatically outperforms due to its superior sub-quadratic (between linear and quadratic) time complexity as well as its effectively negligible memory demand.

physics.chem-phCross-listedarXiv:2610.21075•Oct 16, 2026
High Resolution Image Download MS PowerPoint Slide We report an implementation of the McMurchie–Davidson evaluation scheme for 1- and 2-particle Gaussian AO integrals designed for processors with Single Instruction Multiple Data (SIMD) instruction sets. Like in our recent MD implementation for graphical processing units (GPUs) [ Asadchev, A.; Valeev, E. F. . J. Chem. Phys. 2024, 160, 244109 .], variable-sized batches of shellsets of integrals are evaluated at a time. By optimizing for the floating point instruction throughput rather than minimizing the number of operations, this approach achieves up to 50% of the theoretical hardware peak FP64 performance for many common SIMD-equipped platforms (AVX2, AVX512, NEON), which translates to speedups of up to 30 over the state-of-the-art one-shellset-at-a-time implementation of Obara–Saika-type schemes in Libint for a variety of primitive and contracted integrals. As with our previous work, we rely on the standard C++ programming language─such as the std::simd standard library feature to be included in the 2026 ISO C++ standard─without any explicit code generation to keep the code base small and portable. The implementation is part of the open source LibintX library freely available at https://github.com/ValeevGroup/libintx .
We report an implementation of the McMurchie-Davidson (MD) algorithm for 3-center and 4-center 2-particle integrals over Gaussian atomic orbitals (AOs) with low and high angular momenta l and varying degrees of contraction for graphical processing units (GPUs). This work builds upon our recent implementation of a matrix form of the MD algorithm that is efficient for GPU evaluation of 4-center 2-particle integrals over Gaussian AOs of high angular momenta (l ≥ 4) [A. Asadchev and E. F. Valeev, J. Phys. Chem. A 127, 10889-10895 (2023)]. The use of unconventional data layouts and three variants of the MD algorithm allow for the evaluation of integrals with double precision and sustained performance between 25% and 70% of the theoretical hardware peak. Performance assessment includes integrals over AOs with l ≤ 6 (a higher l is supported). Preliminary implementation of the Hartree-Fock exchange operator is presented and assessed for computations with up to a quadruple-zeta basis and more than 20 000 AOs. The corresponding C++ code is part of the experimental open-source LibintX library available at https://github.com/ValeevGroup/libintx.
In this article, several optimization methods of two-electron repulsion integral calculations on a graphic processing unit (GPU) are presented. These methods are based on the investigations of the method presented by McMurchie and Davidson (MD). A new Boys function evaluation method for the GPU calculation is introduced. The series summation, the error function, and the finite sum formula method are combined; thus, good performance on the GPU can be achieved. By taking some theoretical study of the McMurchie-Davidson recurrence relations, three major optimization approaches are derived from the deduction of the general term formula for the Hermite expansion coefficient. The three approaches contain a new form of the Hermite expansion coefficients with corresponding recurrence relations, which is more efficient for one-electron integrals and [ss|∗∗] or [∗∗|ss] type two-electron integrals. In addition, a simple yet efficient new recurrence formula for the coefficient evaluation is derived, which is more efficient both in float operations and memory operations than its original one. In average, the new recurrence relation can save 26% float operations and 37% memory operations. Finally, a common sub-expression elimination (CSE) method is implemented. This CSE method is directly generated from some equalities we discovered from the general term formula other than by computer algebra system software. This optimized method achieved up to 3.09 speedups compared to the original MD method on the GPU and up to 92.75 speedups compared to the GAMESS calculation on the central processing unit.
Quantum chemistry offers the formal machinery to derive molecular and physical properties arising from (sub)atomic interactions. However, as molecules of practical interest are largely polyatomic, contemporary approximation schemes such as the Hartree–Fock scheme are computationally expensive due to the large number of electron repulsion integrals (ERIs). Central to the Hartree–Fock method is the efficient computation of ERIs over Gaussian functions (GTO-ERIs). Here, the well-known McMurchie–Davidson method (MD) offers an elegant formalism by incrementally extending Hermite Gaussian functions and auxiliary tabulated functions. Although the MD method offers a high degree of versatility to acceleration schemes through Graphics Processing Units (GPUs), the current GPU implementations limit the practical use of supported values of the azimuthal quantum number. In this paper, we propose a generalized framework capable of computing GTO-ERIs for arbitrary azimuthal quantum numbers, provided that the intermediate terms of the MD method can be stored. Our approach benefits from extending the MD recurrence relations through shells, batches, and triple-buffering of the shared memory, and ordering similar ERIs, thus enabling the effective parallelization and use of GPU resources. Furthermore, our approach proposes four GPU implementation schemes considering the suitable mappings between Gaussian basis and CUDA blocks and threads. Our computational experiments involving the GTO-ERI computations of molecules of interest on an NVIDIA A100 Tensor Core GPU (NVIDIA, Santa Clara, CA, USA) have revealed the merits of the proposed acceleration schemes in terms of computation time, including up to a 72× improvement over our previous GPU implementation and up to a 4500× speedup compared to a naive CPU implementation, highlighting the effectiveness of our method in accelerating ERI computations for both monatomic and polyatomic molecules. Our work has the potential to explore new parallelization schemes of distinct and complex computation paths involved in ERI computation.
Libcint is a library designed for the evaluation of analytical integrals for Gaussian type orbitals. It prioritizes simplicity, ease of use, and efficiency for the development of quantum chemistry programs. In the release of version 6.0, Libcint supports the computation of integrals for various operators, such as overlap, Coulomb, Gaunt, Breit, attenuated Coulomb, Slater-type geminals, and Yukawa potential, as well as arbitrary orders of derivatives for these operators. To enhance the usability of the library, Libcint provides a uniform function signature for all integral functions. A code generator is included to automate the implementation of new integrals. To achieve better performance on modern central processing unit architectures, the library employs explicit single instruction multiple data parallelization in the code implementation.
ABSTRACT In Quantum chemical computation, numerical schemes such as the Hartree–Fock (HF) and density functional theory (DFT) are widely used to solve the Schrödinger equation numerically, to realize experiment‐free prediction and analysis of key molecular properties such as structure and energy. Computing one‐electron integrals, such as kinetic energy integrals and nuclear attraction integrals, is essential in both HF and DFT to characterize the molecular electronic states. However, as molecules of practical interest grow in size and angular momentum, computing one‐electron orbitals becomes computationally expensive in most cases. Although computing kinetic energy integrals on CPUs is straightforward, bottlenecks in CPU‐GPU data transfer have often been overlooked. In this study, we propose an efficient method to compute both the kinetic‐energy and nuclear‐attractive integrals on GPUs. First, we explicitly and symbolically expand recurrence relations based on the Obara–Saika and McMurchie–Davidson methods to eliminate redundant operations, thus improving computational efficiency. Second, we implemented a hybrid method that selects the best/fastest of both methods depending on the integration task. Third, we achieved further speedups by using CUDA streams to parallelize the execution of multiple kernels and efficiently utilize multiprocessor resources on the GPU. Computational experiments using NVIDIA A100 GPUs and Intel Xeon Gold 6338 CPU on relevant molecules of interest demonstrated the superiority of our one‐electron integral GPU implementations, achieving a speedup of 20.2 times over PySCF, and a speedup of 132.6 times over GPU4PySCF.
We present a seminumerical exchange algorithm (SNK) based on the sn-LinK method [Laqua et al., J. Chem. Theory Comput. 16, 1456 (2020)], optimized for modern multi-core central processing units with single instruction multiple data vector instructions and report its implementation into the Q-Chem quantum chemistry program. Vectorization and memory management strategies for the Boys function and 3-center-1-electron (3c1e) integral kernels are given. Several new considerations are made in the computation of the 3c1e integrals: the Head-Gordon Pople horizontal recurrence relations, same-center optimizations, and multishell optimizations, all implemented symbolically employing the SymPy symbolic algebra library and exported as optimized C++ code leveraging common-sub-expression elimination (CSE) optimization techniques. A general scheme for choosing optimal integral kernels is developed through rigorous timing benchmarks. The resulting implementation is tested against traditional 4c2e exchange for several systems and basis sets, showing 1-2 orders of magnitude speedups for systems of practically relevant sizes (100-500 atoms, 2000-10000 basis functions). For smaller molecules, comparisons are also made against the more efficient occ-RI-K method showing a break-even point around 50-100 atoms beyond which SNK dramatically outperforms due to its superior sub-quadratic (between linear and quadratic) time complexity as well as its effectively negligible memory demand.