Quantum Package 2.0: An Open-Source Determinant-Driven Suite of Programs

Yann Garniron; Thomas Applencourt; Kevin Gasperich; Anouar; Benali; Anthony Fert\'e; Julien Paquier; Barth\'el\'emy Pradines and; Roland Assaraf; Peter Reinhardt; Julien Toulouse; Pierrette; Barbaresco; Nicolas Renon; Gr\'egoire David; Jean-Paul Malrieu and; Micka\"el V\'eril; Michel Caffarel; Pierre-Fran\c{c}ois Loos and; Emmanuel Giner; Anthony Scemama

arXiv:1902.08154·physics.comp-ph·June 5, 2020

Quantum Package 2.0: An Open-Source Determinant-Driven Suite of Programs

Yann Garniron, Thomas Applencourt, Kevin Gasperich, Anouar, Benali, Anthony Fert\'e, Julien Paquier, Barth\'el\'emy Pradines and, Roland Assaraf, Peter Reinhardt, Julien Toulouse, Pierrette, Barbaresco, Nicolas Renon, Gr\'egoire David, Jean-Paul Malrieu and, Micka\"el V\'eril

PDF

1 Repo

TL;DR

Quantum Package 2.0 is an open-source software suite for quantum chemistry that enables determinant-driven selected configuration interaction methods, supporting large-scale parallel computations and new methodological developments.

Contribution

It introduces a flexible determinant-driven framework with new perturbative and stochastic methods, enhancing efficiency and scalability for quantum chemistry calculations.

Findings

01

Supports runs on thousands of CPU cores with millions of determinants.

02

Achieves high parallel efficiency up to 12,288 cores.

03

Incorporates new renormalized PT2 correction and stochastic CIPSI selection.

Abstract

\textsc{Quantum Package} is an open-source programming environment for quantum chemistry specially designed for wave function methods. Its main goal is the development of determinant-driven selected configuration interaction (sCI) methods and multi-reference second-order perturbation theory (PT2). The determinant-driven framework allows the programmer to include any arbitrary set of determinants in the reference space, hence providing greater methodological freedoms. The sCI method implemented in \textsc{Quantum Package} is based on the CIPSI (Configuration Interaction using a Perturbative Selection made Iteratively) algorithm which complements the variational sCI energy with a PT2 correction. Additional external plugins have been recently added to perform calculations with multireference coupled cluster theory and range-separated density-functional theory. All the programs are…

Tables4

Table 1. Table 1: Zeroth-order energy E ( 0 ) superscript 𝐸 0 E^{(0)} , second-order perturbative correction E ( 2 ) superscript 𝐸 2 E^{(2)} and its renormalized version Z E ( 2 ) 𝑍 superscript 𝐸 2 ZE^{(2)} (in hartree) of CN3 for increasingly large wave functions. The excitation energy Δ E Δ 𝐸 \Delta E (in eV) is the energy difference between the ground state (GS) and the excited state (ES). The statistical error, corresponding to one standard deviation, is reported in parenthesis.

	$E^{(0)}$		$E^{(0)} + E^{(2)}$			$E^{(0)} + Z E^{(2)}$
$N_{det}$	GS (a.u.)	ES (a.u.)	GS (a.u.)	ES (a.u.)	$Δ E$ (eV)	GS (a.u.)	ES (a.u.)	$Δ E$ (eV)
$28$	$- 149.499 574$	$- 149.246 268$	$- 150.155 (1)$	$- 149.863 (1)$	$7.95 (5)$	$- 150.020 (1)$	$- 149.743 (1)$	$7.54 (5)$
$58$	$- 149.519 908$	$- 149.261 390$	$- 150.134 (1)$	$- 149.853 (1)$	$7.67 (5)$	$- 150.018 (1)$	$- 149.744 (1)$	$7.48 (5)$
$131$	$- 149.537 424$	$- 149.277 496$	$- 150.118 (1)$	$- 149.842 7 (9)$	$7.52 (4)$	$- 150.017 (1)$	$- 149.744 9 (9)$	$7.39 (4)$
$268$	$- 149.559 465$	$- 149.298 484$	$- 150.103 5 (9)$	$- 149.830 8 (9)$	$7.42 (4)$	$- 150.015 8 (9)$	$- 149.745 7 (9)$	$7.35 (4)$
$541$	$- 149.593 434$	$- 149.323 302$	$- 150.084 5 (8)$	$- 149.818 6 (8)$	$7.24 (4)$	$- 150.015 2 (8)$	$- 149.746 3 (8)$	$7.32 (4)$
$1 101$	$- 149.627 202$	$- 149.354 807$	$- 150.068 3 (8)$	$- 149.804 5 (8)$	$7.18 (3)$	$- 150.013 7 (8)$	$- 149.746 0 (8)$	$7.28 (3)$
$2 207$	$- 149.663 850$	$- 149.399 522$	$- 150.054 9 (7)$	$- 149.787 9 (7)$	$7.26 (3)$	$- 150.013 2 (7)$	$- 149.746 2 (7)$	$7.27 (3)$
$4 417$	$- 149.714 222$	$- 149.448 133$	$- 150.040 9 (6)$	$- 149.776 2 (6)$	$7.20 (3)$	$- 150.013 0 (6)$	$- 149.747 8 (6)$	$7.22 (3)$
$8 838$	$- 149.765 886$	$- 149.496 401$	$- 150.029 6 (5)$	$- 149.765 5 (5)$	$7.19 (2)$	$- 150.012 4 (5)$	$- 149.747 3 (5)$	$7.21 (2)$
$17 680$	$- 149.817 301$	$- 149.545 048$	$- 150.023 9 (4)$	$- 149.761 5 (4)$	$7.14 (2)$	$- 150.014 1 (4)$	$- 149.750 5 (4)$	$7.17 (2)$
$35 380$	$- 149.859 737$	$- 149.587 668$	$- 150.021 6 (3)$	$- 149.758 2 (3)$	$7.17 (1)$	$- 150.016 1 (3)$	$- 149.751 8 (3)$	$7.19 (1)$
$70 764$	$- 149.893 273$	$- 149.623 235$	$- 150.020 7 (2)$	$- 149.756 6 (3)$	$7.18 (1)$	$- 150.017 4 (2)$	$- 149.753 0 (3)$	$7.19 (1)$
$141 545$	$- 149.919 463$	$- 149.650 109$	$- 150.021 4 (2)$	$- 149.757 2 (2)$	$7.189 (8)$	$- 150.019 4 (2)$	$- 149.755 0 (2)$	$7.196 (8)$
$283 108$	$- 149.937 839$	$- 149.669 735$	$- 150.022 4 (2)$	$- 149.757 6 (2)$	$7.206 (7)$	$- 150.021 1 (2)$	$- 149.756 2 (2)$	$7.209 (7)$
$566 226$	$- 149.950 918$	$- 149.683 278$	$- 150.023 3 (1)$	$- 149.758 0 (1)$	$7.217 (6)$	$- 150.022 3 (1)$	$- 149.757 0 (1)$	$7.219 (6)$
$1 132 520$	$- 149.960 276$	$- 149.693 053$	$- 150.023 8 (1)$	$- 149.758 8 (1)$	$7.212 (5)$	$- 150.023 1 (1)$	$- 149.758 0 (1)$	$7.214 (5)$
$2 264 948$	$- 149.968 203$	$- 149.700 907$	$- 150.024 0 (1)$	$- 149.759 0 (1)$	$7.211 (4)$	$- 150.023 5 (1)$	$- 149.758 4 (1)$	$7.212 (4)$
$4 529 574$	$- 149.975 230$	$- 149.708 061$	$- 150.024 5 (1)$	$- 149.759 4 (1)$	$7.215 (4)$	$- 150.024 1 (1)$	$- 149.758 9 (1)$	$7.216 (4)$
$9 057 914$	$- 149.981 770$	$- 149.714 526$	$- 150.024 63 (9)$	$- 149.759 81 (8)$	$7.206 (3)$	$- 150.024 34 (9)$	$- 149.759 48 (8)$	$7.207 (3)$
$18 110 742$	$- 149.987 928$	$- 149.720 648$	$- 150.024 95 (7)$	$- 149.760 25 (8)$	$7.203 (3)$	$- 150.024 74 (7)$	$- 149.760 00 (8)$	$7.204 (3)$
$36 146 730$	$- 149.993 593$	$- 149.726 253$	$- 150.025 27 (6)$	$- 149.760 65 (7)$	$7.198 (3)$	$- 150.025 02 (6)$	$- 149.760 47 (7)$	$7.198 (3)$

Table 2. Table 2: Zeroth-order energy E ( 0 ) superscript 𝐸 0 E^{(0)} , second-order perturbative correction E ( 2 ) superscript 𝐸 2 E^{(2)} and its renormalized version Z E ( 2 ) 𝑍 superscript 𝐸 2 ZE^{(2)} (in hartree) as a function of the number of determinants N det subscript 𝑁 det N_{\text{det}} for the ground-state of the chromium dimer \ce Cr2 computed in the cc-pVQZ basis set. The statistical error, corresponding to one standard deviation, is reported in parenthesis.

$N_{det}$	$E^{(0)}$	$E^{(0)} + E^{(2)}$	$E^{(0)} + Z E^{(2)}$
$1 631$	$- 2086.742 321$	$- 2087.853 (3)$	$- 2087.679 (2)$
$3 312$	$- 2086.828 496$	$- 2087.821 (2)$	$- 2087.688 (1)$
$6 630$	$- 2086.920 161$	$- 2087.792 (1)$	$- 2087.694 (1)$
$13 261$	$- 2087.008 701$	$- 2087.764 (1)$	$- 2087.694 (1)$
$26 562$	$- 2087.091 669$	$- 2087.743 (1)$	$- 2087.692 (1)$
$53 129$	$- 2087.165 533$	$- 2087.725 (1)$	$- 2087.689 (1)$
$106 262$	$- 2087.234 564$	$- 2087.710 2 (9)$	$- 2087.685 0 (8)$
$212 571$	$- 2087.293 488$	$- 2087.703 0 (8)$	$- 2087.685 0 (7)$
$425 185$	$- 2087.343 762$	$- 2087.697 3 (7)$	$- 2087.684 4 (7)$
$850 375$	$- 2087.386 276$	$- 2087.697 8 (6)$	$- 2087.688 1 (6)$
$1 700 759$	$- 2087.422 707$	$- 2087.698 9 (6)$	$- 2087.691 6 (5)$
$3 401 504$	$- 2087.454 427$	$- 2087.700 7 (5)$	$- 2087.695 1 (5)$
$6 802 953$	$- 2087.482 238$	$- 2087.703 2 (4)$	$- 2087.698 8 (4)$
$13 605 580$	$- 2087.506 838$	$- 2087.705 6 (4)$	$- 2087.702 2 (4)$
$27 210 163$	$- 2087.528 987$	$- 2087.709 2 (4)$	$- 2087.706 4 (4)$
$54 415 174$	$- 2087.549 261$	$- 2087.711 6 (3)$	$- 2087.709 5 (3)$
Extrap.		$- 2087.734$	$- 2087.738$

Table 3. Table 3: Wall-clock time (in seconds) to perform a single Davidson iteration and a second-order correction E ( 2 ) superscript 𝐸 2 E^{(2)} calculation (which also includes the CIPSI selection) with an increasing number of 48-core compute nodes N nodes subscript 𝑁 nodes N_{\text{nodes}} . The statistical error obtained on E ( 2 ) superscript 𝐸 2 E^{(2)} , defining the stopping criterion, is 0.17 × 10 − 3 0.17 superscript 10 3 0.17\times 10^{-3} a.u.

$N_{nodes}$	Wall-clock time (in seconds)
	Davidson for	Davidson for	PT2/selection
	$N_{det} = 25 \times 10^{6}$	$N_{det} = 100 \times 10^{6}$	$N_{det} = 25 \times 10^{6}$
$1$	$3 340$	$65 915$	$406 840$
$32$	$142$	$2 168$	$12 711$
$48$	$109$	$1 497$	$8 515$
$64$	$93$	$1 181$	$6 421$
$96$	$93$	$834$	$4 323$
$128$	$93$	$674$	$3 287$
$192$	$96$	$522$	$2 435$
$256$	$96$	$519$	$1 996$

Table 4. Table 4: Time to access integrals (in nanoseconds/integral) with different access patterns. The time to generate random numbers (measured as 67 ns/integral) was not counted in the random access results.

Access	Array	Hash table
$i, j, k, l$	$9.72$	$125.79$
$i, j, l, k$	$9.72$	$120.64$
$i, k, j, l$	$10.29$	$144.65$
$l, k, j, i$	$88.62$	$125.79$
$l, k, i, j$	$88.62$	$120.64$
Random	$170.00$	$370.00$

Equations60

E_{c} = E_{exact} - E_{HF},

E_{c} = E_{exact} - E_{HF},

ϕ_{p} (r) = μ \sum N_{orb} C_{μ p} χ_{μ} (r) .

ϕ_{p} (r) = μ \sum N_{orb} C_{μ p} χ_{μ} (r) .

∣ * ⟩ Ψ_{k} = I \sum N_{det} c_{I k} ∣ I ⟩,

∣ * ⟩ Ψ_{k} = I \sum N_{det} c_{I k} ∣ I ⟩,

\matrixelement * I \hat{H} I

\matrixelement * I \hat{H} I

\matrixelement * I \hat{H} I_{p}^{r}

\matrixelement * I \hat{H} I_{pq}^{r s}

(p ∣ \hat{h} ∣ q) = \int ϕ_{p} (r) \hat{h} (r) ϕ_{q} (r) d r

(p ∣ \hat{h} ∣ q) = \int ϕ_{p} (r) \hat{h} (r) ϕ_{q} (r) d r

(pq ∣ r s) = \iint ϕ_{p} (r_{1}) ϕ_{q} (r_{1}) r_{12}^{- 1} ϕ_{r} (r_{2}) ϕ_{s} (r_{2}) d r_{1} d r_{2}

(pq ∣ r s) = \iint ϕ_{p} (r_{1}) ϕ_{q} (r_{1}) r_{12}^{- 1} ϕ_{r} (r_{2}) ϕ_{s} (r_{2}) d r_{1} d r_{2}

∣ * ⟩ Ψ^{(0)} = I \sum c_{I} ∣ I ⟩

∣ * ⟩ Ψ^{(0)} = I \sum c_{I} ∣ I ⟩

e_{α}^{(2)} = \frac{\matrixelement * Ψ ^{(0)} H ^ α ^{2}}{E ^{(0)} - \matrixelement * α H ^ α},

e_{α}^{(2)} = \frac{\matrixelement * Ψ ^{(0)} H ^ α ^{2}}{E ^{(0)} - \matrixelement * α H ^ α},

E^{(0)} = \frac{\matrixelement * Ψ ^{(0)} H ^ Ψ ^{(0)}}{\innerproduct * Ψ ^{(0)} Ψ ^{(0)}} \geq E_{FCI}

E^{(0)} = \frac{\matrixelement * Ψ ^{(0)} H ^ Ψ ^{(0)}}{\innerproduct * Ψ ^{(0)} Ψ ^{(0)}} \geq E_{FCI}

E^{(2)} = α \sum e_{α}^{(2)} = α \sum \frac{\matrixelement * α H ^ Ψ ^{(0)} ^{2}}{E ^{(0)} - \matrixelement * α H ^ α}

E^{(2)} = α \sum e_{α}^{(2)} = α \sum \frac{\matrixelement * α H ^ Ψ ^{(0)} ^{2}}{E ^{(0)} - \matrixelement * α H ^ α}

e_{α}^{SHCI} = I max \quantity (\absolutevalue c_{I} \matrixelement * I \hat{H} α) .

e_{α}^{SHCI} = I max \quantity (\absolutevalue c_{I} \matrixelement * I \hat{H} α) .

E = E^{(0)} + E^{(2)},

E = E^{(0)} + E^{(2)},

Σ^{(2)} [E] = α \sum \frac{\matrixelement * α H ^ Ψ ^{(0)} ^{2}}{E - \matrixelement * α H ^ α} .

Σ^{(2)} [E] = α \sum \frac{\matrixelement * α H ^ Ψ ^{(0)} ^{2}}{E - \matrixelement * α H ^ α} .

E = E^{(0)} + Σ^{(2)} [E],

E = E^{(0)} + Σ^{(2)} [E],

Σ^{(2)} [E] \approx Σ^{(2)} [E^{(0)}] + (E - E^{(0)}) \partialderivative Σ^{(2)} [E] E_{E = E^{(0)}} .

Σ^{(2)} [E] \approx Σ^{(2)} [E^{(0)}] + (E - E^{(0)}) \partialderivative Σ^{(2)} [E] E_{E = E^{(0)}} .

E = E^{(0)} + Σ^{(2)} [E^{(0)}] + (E - E^{(0)}) \partialderivative Σ^{(2)} [E] E_{E = E^{(0)}} = E^{(0)} + Z E^{(2)},

E = E^{(0)} + Σ^{(2)} [E^{(0)}] + (E - E^{(0)}) \partialderivative Σ^{(2)} [E] E_{E = E^{(0)}} = E^{(0)} + Z E^{(2)},

Z = \quantity [1 - \partialderivative Σ^{(2)} [E] E_{E = E^{(0)}}]^{- 1},

Z = \quantity [1 - \partialderivative Σ^{(2)} [E] E_{E = E^{(0)}}]^{- 1},

\partialderivative Σ^{(2)} [E] E_{E = E^{(0)}} = - α \sum \frac{\matrixelement * α H ^ Ψ ^{(0)} ^{2}}{( E ^{(0)} - \matrixelement * α H ^ α ) ^{2}} < 0.

\partialderivative Σ^{(2)} [E] E_{E = E^{(0)}} = - α \sum \frac{\matrixelement * α H ^ Ψ ^{(0)} ^{2}}{( E ^{(0)} - \matrixelement * α H ^ α ) ^{2}} < 0.

∣ * ⟩ Ψ^{(0)} = I \in I_{n} \sum c_{I} ∣ I ⟩

∣ * ⟩ Ψ^{(0)} = I \in I_{n} \sum c_{I} ∣ I ⟩

E^{(2)} = α \sum c_{α} \matrixelement * Ψ^{(0)} \hat{H} α

E^{(2)} = α \sum c_{α} \matrixelement * Ψ^{(0)} \hat{H} α

E^{(2)} = I \sum α \in A_{I} \sum c_{α} \matrixelement * Ψ^{(0)} \hat{H} α = I \sum ε_{I},

E^{(2)} = I \sum α \in A_{I} \sum c_{α} \matrixelement * Ψ^{(0)} \hat{H} α = I \sum ε_{I},

c_{α} = \frac{\matrixelement * Ψ ^{(0)} H ^ α}{E ^{(0)} - \matrixelement * α H ^ α}

c_{α} = \frac{\matrixelement * Ψ ^{(0)} H ^ α}{E ^{(0)} - \matrixelement * α H ^ α}

E^{(2)} = J \in D \sum ε_{J} + K \in S \sum ε_{K},

E^{(2)} = J \in D \sum ε_{J} + K \in S \sum ε_{K},

\tilde{e}_{α}^{(2)} = k \sum N_{states} \frac{c _{α k}}{max _{I} c _{I k}^{2}} \matrixelement * Ψ_{k}^{(0)} \hat{H} α,

\tilde{e}_{α}^{(2)} = k \sum N_{states} \frac{c _{α k}}{max _{I} c _{I k}^{2}} \matrixelement * Ψ_{k}^{(0)} \hat{H} α,

c_{α k} = \frac{\matrixelement * Ψ _{k}^{(0)} H ^ α}{\matrixelement * Ψ _{k}^{(0)} H ^ Ψ _{k}^{(0)} - \matrixelement * α H ^ α} .

c_{α k} = \frac{\matrixelement * Ψ _{k}^{(0)} H ^ α}{\matrixelement * Ψ _{k}^{(0)} H ^ Ψ _{k}^{(0)} - \matrixelement * α H ^ α} .

\tilde{H} = H + γ \quantity (S^{2} - I ⟨ S^{2} ⟩_{target})^{2},

\tilde{H} = H + γ \quantity (S^{2} - I ⟨ S^{2} ⟩_{target})^{2},

P_{I} [i] = 1 \land k = 0 \sum i I [k],

P_{I} [i] = 1 \land k = 0 \sum i I [k],

{+ (- 1)^{P_{I} [h - 1] \oplus P_{I} [p - 1]}, - (- 1)^{P_{I} [h - 1] \oplus P_{I} [p - 1]}, if p > h, if h > p,

{+ (- 1)^{P_{I} [h - 1] \oplus P_{I} [p - 1]}, - (- 1)^{P_{I} [h - 1] \oplus P_{I} [p - 1]}, if p > h, if h > p,

∣ * ⟩ Ψ^{(0)} = I \sum N_{det}^{↑} J \sum N_{det}^{↓} C_{I J} ∣ I_{↑} J_{↓} ⟩,

∣ * ⟩ Ψ^{(0)} = I \sum N_{det}^{↑} J \sum N_{det}^{↓} C_{I J} ∣ I_{↑} J_{↓} ⟩,

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Code & Models

Repositories

QuantumPackage/qp2
noneOfficial

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Quantum Package 2.0: An Open-Source Determinant-Driven Suite of Programs

Yann Garniron

Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, France

Thomas Applencourt

Computational Science Division, Argonne National Laboratory, Argonne, IL 60439, United States

Kevin Gasperich

Computational Science Division, Argonne National Laboratory, Argonne, IL 60439, United States

Department of Chemistry, University of Pittsburgh, Pittsburgh, PA 15260, United States

Anouar Benali

Computational Science Division, Argonne National Laboratory, Argonne, IL 60439, United States

Anthony Ferté

Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France

Julien Paquier

Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France

Barthélémy Pradines

Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France

Institut des Sciences du Calcul et des Données, Sorbonne Université, F-75005 Paris, France

Roland Assaraf

Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France

Peter Reinhardt

Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France

Julien Toulouse

Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France

Pierrette Barbaresco

CALMIP, Université de Toulouse, CNRS, INPT, INSA, UPS, UMS 3667, France

Nicolas Renon

CALMIP, Université de Toulouse, CNRS, INPT, INSA, UPS, UMS 3667, France

Grégoire David

Aix-Marseille Univ, CNRS, ICR, Marseille, France

Jean-Paul Malrieu

Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, France

Mickaël Véril

Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, France

Michel Caffarel

Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, France

Pierre-François Loos

[email protected]

Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, France

Emmanuel Giner

[email protected]

Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France

Anthony Scemama

[email protected]

Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, France

Abstract

Quantum chemistry is a discipline which relies heavily on very expensive numerical computations. The scaling of correlated wave function methods lies, in their standard implementation, between $\order*{N^{5}}$ and $\order*{e^{N}}$ , where $N$ is proportional to the system size. Therefore, performing accurate calculations on chemically meaningful systems requires i) approximations that can lower the computational scaling, and ii) efficient implementations that take advantage of modern massively parallel architectures. Quantum Package is an open-source programming environment for quantum chemistry specially designed for wave function methods. Its main goal is the development of determinant-driven selected configuration interaction (sCI) methods and multi-reference second-order perturbation theory (PT2). The determinant-driven framework allows the programmer to include any arbitrary set of determinants in the reference space, hence providing greater methodological freedom. The sCI method implemented in Quantum Package is based on the CIPSI (Configuration Interaction using a Perturbative Selection made Iteratively) algorithm which complements the variational sCI energy with a PT2 correction. Additional external plugins have been recently added to perform calculations with multireference coupled cluster theory and range-separated density-functional theory. All the programs are developed with the IRPF90 code generator, which simplifies collaborative work and the development of new features. Quantum Package strives to allow easy implementation and experimentation of new methods, while making parallel computation as simple and efficient as possible on modern supercomputer architectures. Currently, the code enables, routinely, to realize runs on roughly 2 000 CPU cores, with tens of millions of determinants in the reference space. Moreover, we have been able to push up to 12 288 cores in order to test its parallel efficiency. In the present manuscript, we also introduce some key new developments: i) a renormalized second-order perturbative correction for efficient extrapolation to the full CI limit, and ii) a stochastic version of the CIPSI selection performed simultaneously to the PT2 calculation at no extra cost.

I Introduction

In 1965, Gordon Moore predicted that the number of transistors in an integrated circuit would double about every two years (the so-called Moore’s law). Moore (1965) Rapidly, this “law” was interpreted as an expected two-fold increase in performance every 18 months. This became an industrial goal. The development of today’s most popular electronic structure codes was initiated in the 1990’s (or even before). At that time, the increase of computational power from one supercomputer generation to the next was mostly driven by an increase of processors’ frequency. Indeed, the amount of random access memory was small, the time to access data from disk was slow, and the energy consumption of the most powerful computer was 236 kW, hence far from being an economical concern. top (2018) At the very beginning of the 21st century, having increased continuously, both the number of processors and their frequency raised the supercomputer power consumption by two orders of magnitude, inflating accordingly the electricity bill. The only way to slow down this frenetic growth of power consumption while keeping alive Moore’s dream was to freeze the processor’s frequency (between 1 and 4 GHz), and increase the number of CPU cores. The consequence of such a choice was that “free lunch” was over: the programmers now had to parallelize their programs to make them run faster. Sutter and Larus (2005) At the same time, computer scientists realized that the increase of performance in memory access was slower than the increase in computational power, Wulf and McKee (1995) and that the floating-point operation (or flop) count would soon stop being the bottleneck. From now on, data movement would be the main concern. This paradigm shift was named the memory wall. Moore’s law is definitely near the end of its life. Khan et al. (2018)

The traditional sequential algorithms of quantum chemistry are currently being redesigned and replaced by parallel equivalents by multiple groups around the world. Booth et al. (2009); Booth and Alavi (2010); Cleland et al. (2010); Sharma et al. (2017); Garniron et al. (2017); Smith et al. (2017); Neuhauser et al. (2013); Willow and Hirata (2014); Willow et al. (2012); Johnson et al. (2017, 2016); Gruneis et al. (2017); Doran and Hirata (2016) This has obviously a significant influence on methodological developments. The most iconic example of this move towards parallel-friendly methods is the recently developed full configuration interaction quantum Monte Carlo (FCIQMC) method by Alavi and coworkers. Booth et al. (2009) FCIQMC can be interpreted as a Monte Carlo equivalent of older selected configuration interaction (sCI) algorithms Bender and Davidson (1969); Whitten and Hackmeyer (1969); Huron et al. (1973); Shih et al. (1978); Buenker et al. (1978); Evangelisti et al. (1983); Cimiraglia (1985); Cimiraglia and Persico (1987); Illas et al. (1988); Povill et al. (1992); Engels et al. (2001); Abrams and Sherrill (2005); Bunge and Carbó-Dorca (2006); Musch and Engels (2006); Bytautas and Ruedenberg (2009); Giner et al. (2013); Caffarel et al. (2014); Giner et al. (2015); Garniron et al. (2017); Caffarel et al. (2016, 2016); Holmes et al. (2016); Sharma et al. (2017); Holmes et al. (2017); Chien et al. (2018); Scemama et al. (2018, 2018); Loos et al. (2018); Garniron et al. (2018); Evangelista (2014); Schriber and Evangelista (2016, 2017); Liu and Hoffmann (2016); Per and Cleland (2017); Ohtsuka and ya Hasegawa (2017); Zimmerman (2017); Li et al. (2018); Loos et al. (2019) such as CIPSI (Configuration Interaction using a Perturbative Selection made Iteratively),Huron et al. (1973) that are iterative and thus a priori not well adapted to massively parallel architecture. As we shall see here, things turn out differently, and the focus of the present article is to show that sCI methods can be made efficient on modern massively parallel supercomputers.

Quantum Package QP (2019) is an open-source suite of wave function quantum chemistry methods mainly developed at the Laboratoire de Chimie et Physique Quantiques (LCPQ) in Toulouse (France), and the Laboratoire de Chimie Théorique (LCT) in Paris. Its source code is freely available on GitHub at the following address: https://github.com/QuantumPackage/qp2. Quantum Package strives to allow easy implementation and experimentation of new methods, while making parallel computation as simple and efficient as possible. Accordingly, the initial choice of Quantum Package was to go towards determinant-driven algorithms. Assuming a wave function expressed as a linear combination of determinants, a determinant-driven algorithm essentially implies that the outermost loop runs over determinants. On the other hand, more traditional integral-driven algorithms have their outermost loop running on the two-electron integrals appearing in the expression of the matrix elements in the determinant basis (see Sec. II.2). Determinant-driven algorithms allow more flexibility than their integral-driven counterparts, Povill and Rubio (1995) but they have been known for years to be less efficient than their integral-driven variant for solving electronic structure problems. In high-precision calculations, the number of determinants is larger than the number of integrals, justifying the integral-driven choice. However, today’s programming standards impose parallelism, and if determinant-driven calculations prove to be better adapted to parallelism, such methods could regain popularity. More conventional approaches have also been very successfully parallelized: CCSD(T), Olson et al. (2007); Kjaergaard et al. (2017) DMRG, Kantian et al. (2019) GW, Blase et al. (2018) QMC, Scemama et al. (2013, 2016); Kim et al. (2018) and many others.

Quantum Package was used in numerous applications, in particular to obtain reference ground-state energies Scemama et al. (2014); Caffarel et al. (2014); Giner et al. (2013, 2015); Caffarel et al. (2016, 2016) as well as excitation energies Giner et al. (2018); Loos et al. (2018, 2019) for atomic and molecular systems. For example, in Ref. Loos et al., 2018, Quantum Package has been used to compute more than hundred very accurate transition energies for states of various characters (valence, Rydberg, $n\rightarrow\pi^{*}$ , $\pi\rightarrow\pi^{*}$ , singlet, triplet, …) in 18 small molecules. The high quality and compactness of the CIPSI wave function was also used for quantum Monte Carlo calculations to characterize the ground state of the water and the \ceFeS molecules, Caffarel et al. (2016); Scemama et al. (2018) and obtained highly accurate excitation energies. Scemama et al. (2018); Dash et al. (2018); Pineda Flores and Neuscamman (2018) Of course, the technical considerations were not the main concern of the different articles that were produced. Because the present work focused on the actual implementation of the methods at least as much as on the theory behind them, this article is a perfect opportunity to discuss in depth their implementation.

This manuscript is organized as follows. In Sec. II, we briefly describe the main computational methods implemented in Quantum Package as well as newly developed methods and extrapolation techniques. Section III deals with their implementation. In particular, Sec. III.1 discusses the computation of the Hamiltonian matrix elements using determinant-driven algorithms, while Sec. III.3 focuses on the acceleration of the Davidson diagonalization, a pivotal point of sCI methods. In Sec. III.4, we focus on the determinant selection step used to build compact wave functions. In a nutshell, the principle is to incrementally build a reference wave function by scavenging its external space for determinants that interact with it. To make this step more affordable, we designed a new stochastic scheme which selects on the fly the more important determinants while the second-order perturbative (PT2) energy is computed using a hybrid stochastic-deterministic scheme. Garniron et al. (2017) Therefore, the selection part of this new stochastic CIPSI selection is virtually free as long as one is interested in the second-order perturbative correction, which is crucial in many cases in order to obtain near full configuration interaction (FCI) results. Section IV briefly explains how we produce spin-adapted wave functions, and Sec. V describes parallelism within Quantum Package. The efficiency of the present algorithms is demonstrated in Sec. VI.3 where illustrative calculations and parallel speedups are reported. Finally, Sec. VII discusses the development philosophy of Quantum Package as well as other relevant technical details. Unless otherwise stated, atomic units are used throughout.

II Methods

II.1 Generalities

The correlation energy is defined as Löwdin (1959)

[TABLE]

where $E_{\text{exact}}$ and $E_{\text{HF}}$ are, respectively, the exact (non-relativistic) energy and the Hartree-Fock (HF) energy in a complete (one-electron) basis set.

To include electron correlation effects, the wave function associated with the $k$ th electronic state, $\ket*{\Psi_{k}}$ , may be expanded in the set of all possible $N$ -electron Slater determinants, $\ket{I}$ , built by placing $N_{\uparrow}$ spin-up electrons in $N_{\text{orb}}$ orbitals and $N_{\downarrow}$ spin-down electrons in $N_{\text{orb}}$ orbitals (where $N=N_{\uparrow}+N_{\downarrow}$ ). These so-called molecular orbitals (MOs) are defined as linear combinations of atomic orbitals (AOs)

[TABLE]

Note that the MOs are assumed to be real valued in the context of this work. The eigenvectors of the Hamiltonian $\hat{H}$ are consequently expressed as linear combinations of Slater determinants, i.e.,

[TABLE]

where $N_{\text{det}}$ is the number of determinants. For sake of conciseness, we will restrict the discussion to the ground state (i.e. $k=0$ ) and drop the subscript $k$ accordingly. Solving the eigenvalue problem in this basis is referred to as FCI and yields, for a given basis set, the exact solution of the Schrödinger equation. Unfortunately, FCI is usually computationally intractable because of its exponential scaling with the size of the system.

II.2 Matrix elements of the Hamiltonian

In the $N$ -electron basis of Slater determinants, one expects the matrix elements of $\hat{H}$ to be integrals over $3N$ dimensions. However, given the two-electron nature of the Hamiltonian, and because the MOs are orthonormal, Slater determinants that differ by more than two spinorbitals yield a zero matrix element. The remaining elements can be expressed as sums of integrals over one- or two-electron coordinates, which can be computed at a reasonable cost. These simplifications are known as Slater-Condon’s rules, and reads

[TABLE]

where $\hat{h}$ is the one-electron part of the Hamiltonian (including kinetic energy and electron-nucleus attraction operators),

[TABLE]

are one-electron integrals, $i\in\ket{I}$ means that $\phi_{i}$ belongs to the Slater determinant $\ket{I}$ , $\ket*{I_{p}^{r}}$ and $\ket*{I_{pq}^{rs}}$ are determinants obtained from $\ket{I}$ by substituting orbitals $\phi_{p}$ by $\phi_{r}$ , and $\phi_{p}$ and $\phi_{q}$ by $\phi_{r}$ and $\phi_{s}$ , respectively,

[TABLE]

are two-electron electron repulsion integrals (ERIs), $r_{12}^{-1}=\absolutevalue{{\mathbf{r}}_{1}-{\mathbf{r}}_{2}}^{-1}$ is the Coulomb operator, and $(pq||rs)=(pq|rs)-(ps|rq)$ are the usual antisymmetrized two-electron integrals.

Within the HF method, Roothaan’s equations allow to solve the problem in the AO basis.Roothaan (1951) In this context, one needs to compute the $\order*{N_{\text{orb}}^{4}}$ two-electron integrals $(\mu\nu|\lambda\sigma)$ over the AO basis. Thanks to a large effort in algorithmic development and implementation,Obara and Saika (1986); Head-Gordon and Pople (1988); Ten-no (1993); Gill et al. (1989); Gill and Pople (1991); Valeev (2018); Barca and Loos (2017); Zhang (2018) these integrals can now be computed very fast on modern computers. However, with post-HF methods, the computation of the two-electron integrals is a potential bottleneck. Indeed, when computing matrix elements of the Hamiltonian in the basis of Slater determinants, ERIs over MOs are required. Using Eq. (2), the cost of computing a single integral $(pq|rs)$ scales as $\order*{N_{\text{orb}}^{4}}$ . A naive computation of all integrals in the MO basis would cost ${\cal O}(N_{\text{orb}}^{8})$ . Fortunately, computing all of them can be scaled down to ${\cal O}(N_{\text{orb}}^{5})$ by transforming the indices one by one.Wilson (1987) This step is known as the four-index integral transformation. In addition to being very costly, this step is hard to parallelize in a distributed way, because it requires multiple collective communications.Rajbhandari et al. (2017); Limaye and Gadre (1994); Fletcher et al. (1999); Covick and Sando (1990) However, techniques such as density fitting (also called the resolution of the identity), Whitten (1973); Eichkorn et al. (1995); Schmitz et al. (2017) low-rank approximations, Beebe and Linderberg (1977); Aquilante et al. (2007); Røeggen and Johansen (2008); Peng and Kowalski (2017) or the combination of both Pham and Gordon (2019) are now routinely employed to overcome the computational and storage bottlenecks.

II.3 Selected CI methods

The sCI methods rely on the same principle as the usual configuration interaction (CI) approaches, except that determinants are not chosen a priori based on occupation or excitation criteria, but selected among the entire set of determinants based on their estimated contribution to the FCI wave function. Indeed, it has been noticed long ago that, even inside a predefined subspace of determinants, only a small number of them significantly contributes.Bytautas and Ruedenberg (2009); Anderson et al. (2018) Therefore, an on-the-fly selection of determinants is a rather natural idea that has been proposed in the late 1960’s by Bender and Davidson Bender and Davidson (1969) as well as Whitten and Hackmeyer.Whitten and Hackmeyer (1969) sCI methods are still very much under active development. The main advantage of sCI methods is that no a priori assumption is made on the type of electronic correlation. Therefore, at the price of a brute force calculation, a sCI calculation is less biased by the user’s appreciation of the problem’s complexity.

The approach that we have implemented in Quantum Package is based on the CIPSI algorithm developed by Huron, Rancurel and Malrieu in 1973, Huron et al. (1973) that iteratively selects external determinants ${\ket{\alpha}}$ — determinants which are not present in the (reference or variational) zeroth-order wave function

[TABLE]

at a given iteration — using a perturbative criterion

[TABLE]

where

[TABLE]

is the zeroth-order (variational) energy, and $e_{\alpha}^{(2)}$ the (second-order) estimated gain in correlation energy that would be brought by the inclusion of ${\ket{\alpha}}$ . The second-order perturbative correction

[TABLE]

is an estimate of the total missing correlation energy, i.e., $E^{(2)}\approx E_{\text{FCI}}-E^{(0)}$ , for large enough expansions.

Let us emphasize that sCI methods can be applied to any determinant space. Although presented here for the FCI space, it can be trivially generalized to a complete active space (CAS), but also to standard CI spaces such as CIS, CISD or MR-CISD. The only required modification is to set to zero the contributions associated with the determinants which do not belong to the target space.

There is, however, a computational downside to sCI methods. In conventional CI methods, the rule by which determinants are selected is known a priori, and therefore, one can map a particular determinant to some row or column indices.Knowles and Handy (1984) As a consequence, it can be systematically determined to which matrix element of $\hat{H}$ a two-electron integral contributes. This allows for the implementation of so-called integral-driven methods that work essentially by iterating over integrals. On the contrary, in (most) sCI methods, the determinants are selected a posteriori, and an explicit list has to be maintained as there is no immediate way to know whether or not a determinant has been selected. Consequently, we must rely on the so-called determinant-driven approach in which iterations are performed over determinants rather than integrals. This can be a lot more expensive, since the number of determinants $N_{\text{det}}$ is typically much larger than the number of integrals. The number of determinants scales as $\order*{N_{\text{orb}}!}$ while the number of integrals scales (formally) as $\order*{N_{\text{orb}}^{4}}$ . What makes sCI calculations possible in practice is that sCI methods generate relatively compact wave functions, i.e. wave functions where $N_{\text{det}}$ is much smaller (by orders of magnitude) than the size of the FCI space. Furthermore, determinant-driven methods require an effective way to compare determinants in order to extract the corresponding excitation operators, and a way to rapidly fetch the associated integrals involved, as described in Sec. III.1.

Because of this high computational cost, approximations have been proposed.Evangelisti et al. (1983) Recently, the semi-stochastic heat-bath configuration interaction (SHCI) algorithm has taken further the idea of a more approximate but extremely cheap selection. Holmes et al. (2016); Sharma et al. (2017); Li et al. (2018) Compared to CIPSI, the selection criterion is simplified to

[TABLE]

This algorithmically allows for an extremely fast selection of doubly-excited determinants by an integral-driven approach. Nonetheless, the bottlenecks of the SHCI are the diagonalization step and the computation of $E^{(2)}$ , which remain determinant driven.

As mentioned above, FCIQMC is an alternative approach of stochastic nature recently developed in Alavi’s group, Booth et al. (2009); Booth and Alavi (2010); Cleland et al. (2010) where signed walkers spawn from one determinant to connected ones, with a probability that is a function of the associated matrix element. The average proportion of walkers on a determinant converges to its coefficient in the FCI wave function. A more “brute force” approach is the purely stochastic selection of Monte Carlo CI (MCCI), Greer (1995, 1998) where determinants are randomly added to the zeroth-order wave function. After diagonalization, the determinants of smaller coefficient are removed, and new random determinants are added.

A number of other variants have been developed but are not detailed here. Bender and Davidson (1969); Whitten and Hackmeyer (1969); Huron et al. (1973); Evangelisti et al. (1983); Cimiraglia (1985); Cimiraglia and Persico (1987); Illas et al. (1988); Povill et al. (1992); Abrams and Sherrill (2005); Bunge and Carbó-Dorca (2006); Bytautas and Ruedenberg (2009); Giner et al. (2013); Caffarel et al. (2014); Giner et al. (2015); Garniron et al. (2017); Caffarel et al. (2016, 2016); Holmes et al. (2016); Sharma et al. (2017); Holmes et al. (2017); Chien et al. (2018); Scemama et al. (2018, 2018); Loos et al. (2018); Garniron et al. (2018); Evangelista (2014); Schriber and Evangelista (2016, 2017); Liu and Hoffmann (2016); Per and Cleland (2017); Ohtsuka and ya Hasegawa (2017); Zimmerman (2017); Li et al. (2018); Ohtsuka and ya Hasegawa (2017); Coe (2018); Loos et al. (2019) Although these various approaches appear under diverse acronyms, most of them rely on the very same idea of selecting determinants iteratively according to their contribution to the wave function or energy.

II.4 Extrapolation techniques

II.4.1 Usual extrapolation procedure

In order to extrapolate the sCI results to the FCI limit, we have adopted the method recently proposed by Holmes, Umrigar and Sharma Holmes et al. (2017) in the context of the SHCI method. Holmes et al. (2016); Sharma et al. (2017); Holmes et al. (2017) It consists of extrapolating the sCI energy, $E^{(0)}$ , as a function of the second-order Epstein-Nesbet energy, $E^{(2)}$ , which is an estimate of the truncation error in the sCI algorithm, i.e $E^{(2)}\approx E_{\text{FCI}}-E^{(0)}$ . Huron et al. (1973) When $E^{(2)}=0$ , the FCI limit has effectively been reached. This extrapolation procedure has been shown to be robust, even for challenging chemical situations. Holmes et al. (2017); Sharma et al. (2017); Scemama et al. (2018, 2018); Chien et al. (2018); Garniron et al. (2018); Loos et al. (2018, 2019) Below, we propose an improved extrapolation scheme which renormalizes the second-order perturbative correction.

II.4.2 Renormalized PT2

At a given sCI iteration, the sCI+PT2 energy is given by

[TABLE]

where $E^{(0)}$ and $E^{(2)}$ are given by Eqs. (9) and (10), respectively. Let us introduce the following energy-dependent second-order self-energy

[TABLE]

Obviously, we have $\Sigma^{(2)}[E^{(0)}]=E^{(2)}$ . Now, let us consider the more general problem, which is somewhat related to Brillouin-Wigner perturbation theory, where we have

[TABLE]

and assume that $\Sigma^{(2)}[E]$ behaves linearly for $E\approx E^{(0)}$ , i.e.,

[TABLE]

This linear behavior is corroborated by the findings of Nitzsche and Davidson. Nitzsche and Davidson (1978) Substituting Eq. (15) into (14) yields

[TABLE]

where the renormalization factor is

[TABLE]

and

[TABLE]

Therefore, the renormalization factor fulfills the condition $0\leq Z\leq 1$ , and its actual computation does not involve any additional cost when computed alongside $E^{(2)}$ as they involve the same quantities. This renormalization procedure of the second-order correction, that we have named rPT2, bears obvious similarities with the computation of quasiparticle energies within the G0W0 method. Onida et al. (2002); Reining (2017); Loos et al. (2018); Veril et al. (2018) Practically, the effect of rPT2 is to damp the value of $E^{(2)}$ for small wave functions. Indeed, when $N_{\text{det}}$ is small, the sum $E^{(0)}+E^{(2)}$ usually overestimates (in magnitude) the FCI energy, yielding a pronounced non-linear behavior of the sCI+PT2 energy. Consequently, by computing instead the (renormalized) energy $E^{(0)}+Z\,E^{(2)}$ , one observes a much more linear behavior of the energy, hence an easier extrapolation to the FCI limit. Its practical usefulness is illustrated in Sec. VI.2.

III Implementation

In this section, we give an overview of the implementation of the various methods present in Quantum Package. The implementation of the crucial algorithms is explained in detail in the PhD thesis of Dr Y. Garniron Garniron (2019) as well as in the Appendix of the present manuscript.

III.1 Determinant-driven computation of the matrix elements

For performance sake, it is vital that some basic operations are done efficiently and, notably, the computation of the Hamiltonian matrix elements. This raises some questions about the data structures chosen to represent the two-electron integrals and determinants, as well as their consequences from an algorithmic point of view. This section is going to address these questions by going through the basic concepts of our determinant-driven approach.

III.1.1 Storage of the two-electron integrals

In Quantum Package, the two-electron integrals are kept in memory because they require a fast random access. Considering the large number of two-electron integrals, a hash table is the natural choice allowing the storage of only non-zero values with a data retrieval in near constant time. Maurer and Lewis (1975) However, standard hashing algorithms tend to shuffle data to limit the probability of collisions. Here, we favor data locality using the hash function given in Algorithm 1. This hash function i) returns the same value for all keys related by permutation symmetry, ii) keeps some locality in the storage of data, and iii) can be evaluated in 10 CPU cycles (estimated with MAQAODjoudi et al. (2005)) if the integer divisions by two are replaced by right bit shift instructions.

The hash table is such that each bucket can potentially store $2^{15}$ consecutive key-value pairs. The 15 least significant bits of the hash value are removed to give the bucket index [ $i_{\text{bucket}}=\lfloor\text{hash}(i,j,k,l)/2^{15}\rfloor$ ], and only those 15 bits need to be stored in the bucket for the key storage [ $\text{hash}(i,j,k,l)\mod 2^{16}$ ]. Hence, the key storage only requires two bytes per key, and they are sorted in increasing order, enabling a binary search within the bucket. The key search is always fast since the binary search is bounded by 15 misses and the maximum size of the key array is 64 kiB, the typical size of the L1 cache. The efficiency of the integral storage is illustrated in Appendix A.1.

III.2 Internal representation of determinants

Determinants can be conveniently written as a string of creation operators applied to the vacuum state ${\ket{}}$ , e.g., $a^{\dagger}_{i}a^{\dagger}_{j}a^{\dagger}_{k}{\ket{}}=\ket{I}$ . Because of the fermionic nature of electrons, a permutation of two contiguous creation operators results in a sign change $a^{\dagger}_{j}a^{\dagger}_{i}=-a^{\dagger}_{i}a^{\dagger}_{j}$ , which makes their ordering relevant, e.g., $a^{\dagger}_{j}a^{\dagger}_{i}a^{\dagger}_{k}{\ket{}}=-\ket{I}$ . A determinant can be broken down into two pieces of information: i) a set of creation operators corresponding to the set of occupied spinorbitals in the determinant, and ii) an ordering of the creation operators responsible for the sign of the determinant, known as phase factor. Once an ordering operator ${\hat{\mathcal{O}}}$ is chosen and applied to all determinants, the phase factor may simply be included in the CI coefficient.

The determinants are built using the following order: i) spin-up ( $\uparrow$ ) spinorbitals are placed before spin-down ( $\downarrow$ ) spinorbitals, as in the Waller-Hartree double determinant representationPauncz (1989) ${\hat{\mathcal{O}}}\ket{I}=\hat{I}{\ket{}}=\hat{I}_{\uparrow}\hat{I}_{\downarrow}{\ket{}}$ , and ii) within each operator $\hat{I}_{\uparrow}$ and $\hat{I}_{\downarrow}$ , the creation operators are sorted by increasing indices. For instance, let us consider the determinant $\ket{J}=a^{\dagger}_{j}a^{\dagger}_{k}a^{\dagger}_{\bar{i}}a^{\dagger}_{i}{\ket{}}$ built from the set of spinorbitals $\{i_{\uparrow},j_{\uparrow},k_{\uparrow},i_{\downarrow}\}$ with $i<j<k$ . If we happen to encounter such a determinant, our choice of representation imposes to consider its re-ordered expression ${\hat{\mathcal{O}}}\ket{J}=-a^{\dagger}_{i}a^{\dagger}_{j}a^{\dagger}_{k}a^{\dagger}_{\bar{i}}{\ket{}}=-\ket{J}$ , and the phase factor must be handled.

The indices of the creation operators (or equivalently the spinorbital occupations), are stored using the so-called bitstring encoding. A bitstring is an array of bits; typically, the 64-bit binary representation of an integer is a bitstring of size 64. Quite simply, the idea is to map each spinorbital to a single bit with value set to its occupation number. In other words, 0 and 1 are associated with the unoccupied and occupied states, respectively. Additional information about the internal representation of determinants can be found in Appendix A.2.

III.3

Davidson diagonalization

Finding the lowest root(s) of the Hamiltonian is a necessary step in CI methods. Standard diagonalization algorithms scale as $\order*{N_{\text{det}}^{3}}$ and $\order*{N_{\text{det}}^{2}}$ in terms of computation and storage, respectively. Hence, their cost is prohibitive as $N_{\text{det}}$ is usually, at least, of the order of few millions. Fortunately, not all the spectrum of $\hat{H}$ is required: only the first few lowest eigenstates are of interest. The Davidson diagonalization Davidson (1975); Liu (1978); Olsen et al. (1990); Gadea (1994); Crouzeix et al. (1994) is an iterative algorithm which aims at extracting the first $N_{\text{states}}$ lowest eigenstates of a large matrix. This algorithm reduces the cost of both the computation and storage to $\order*{N_{\text{states}}N_{\text{det}}^{2}}$ and $\order*{N_{\text{states}}N_{\text{det}}}$ , respectively. It is presented as Algorithm 2 and further details about the present Davidson algorithm implementation are gathered in Appendix A.3.

III.4

CIPSI selection and PT2 energy

III.4.1 The basic algorithm

The largest amount of work for this second version of Quantum Package has been devoted to the improvement of the CIPSI algorithm implementation. Giner (2014) As briefly described in Sec. II, this is an iterative selection algorithm, where determinants are added to the reference wave function according to a perturbative criterion.

The $n$ th CIPSI iteration can be described as follows:

The zeroth-order (reference or variational) wave function

[TABLE]

is defined over a set of determinants $\mathcal{I}_{n}$ — characterized as internal determinants — from which the lowest eigenvector of $\hat{H}$ are obtained. 2. 2.

For all external determinants ${\ket{\alpha}}\notin\mathcal{I}_{n}$ but connected to $\mathcal{I}_{n}$ , i.e., $\matrixelement*{\Psi^{(0)}}{\hat{H}}{\alpha}\neq 0$ , we compute the individual perturbative contribution $e_{\alpha}^{(2)}$ given by Eq. (8). This set of external determinants is labeled $\mathcal{A}_{n}$ . 3. 3.

Summing the contributions of all the external determinants $\alpha\in\mathcal{A}_{n}$ gives the second-order perturbative correction provided by Eq. (10) and the FCI energy can be estimated as $E_{\text{FCI}}\approx E^{(0)}+E^{(2)}$ . 4. 4.

We extract ${\ket{\alpha^{\star}}}\in\mathcal{A}^{\star}_{n}$ , the subset of determinants ${\ket{\alpha}}\in\mathcal{A}_{n}$ with the largest contributions $e_{\alpha}^{(2)}$ , and add them to the variational space $\mathcal{I}_{n+1}=\mathcal{I}_{n}\cup\mathcal{A}^{\star}_{n}$ . In practice, in the case of the single-state calculation, we aim at doubling the size of the reference wave function at each iteration. 5. 5.

Iterate until the desired convergence has been reached.

All the details of our current implementation are reported in Appendix A.4. In the remaining of this section, we only discuss the algorithm of our new stochastic CIPSI selection.

III.4.2 New stochastic selection

In the past, CIPSI calculations were only possible in practice thanks to approximations. The first approximation restricts the set $\mathcal{A}_{n}$ by defining a set of generators. Indeed, it is very unlikely that ${\ket{\alpha}}$ will be selected if it is not connected to any $\ket{I}$ with a large coefficient, so only the determinants with the largest coefficients are generators. A second approximation defines a set of selectors in order to reduce the cost of $e_{\alpha}^{(2)}$ by removing the determinants with the smallest coefficients in the expression of $\Psi^{(0)}$ in $E^{(2)}$ . This approximate scheme was introduced in the 80’s and is known as three-class CIPSI.Evangelisti et al. (1983) The downside of these approximations is that the calculation is biased and, consequently, does not strictly converge to the FCI limit. Moreover, similar to the initiator approximation in FCIQMC, Cleland et al. (2010) this scheme suffers from a size-consistency issue. Ten-no (2017) The stochastic selection that we describe in this section (asymptotically) cures this problem, as there is no threshold on the wave function: if the calculation is run long enough, the unbiased FCI solution is obtained.

Recently, some of us developed a hybrid deterministic/stochastic algorithm for the computation of $E^{(2)}$ .Garniron et al. (2017) The main idea is to rewrite the expression of

[TABLE]

into elementary contributions labeled by the determinants of the internal space:

[TABLE]

where

[TABLE]

is the corresponding coefficient estimated via first-order perturbation theory, and $\mathcal{A}_{I}$ is the subset of determinants ${\ket{\alpha}}$ connected to $\ket{I}$ by $\hat{H}$ such that ${\ket{\alpha}}\notin\cup_{K<I}\mathcal{A}_{K}$ . The sum is decomposed into a stochastic and a deterministic contribution

[TABLE]

where $\mathcal{D}$ and $\mathcal{S}$ are the sets of determinants included in the deterministic and stochastic components, respectively.

The $\ket{I}$ ’s are sorted by decreasing $c_{I}^{2}$ , and two processes are used simultaneously to compute the contributions $\varepsilon_{I}$ . The first process is stochastic and $\ket{I}$ is drawn according to $c_{I}^{2}$ . When a given $\varepsilon_{I}$ has been computed once, its contribution is stored such that if $\ket{I}$ is drawn again later the contribution does not need to be recomputed. The only update is to increment the number of times it has been drawn for the Monte Carlo statistics. In parallel, a deterministic process is run, forcing to compute the contribution $\varepsilon_{I}$ with the smallest index which has yet to be computed. The deterministic component is chosen as the first contiguous set of $\varepsilon_{I}$ . Hence, the computation of $E^{(2)}$ is unbiased, and the exact deterministic value can be obtained in a finite time if the calculation is run long enough. The stochastic part is only a convergence accelerator providing a reliable error bar. The computation of $E^{(2)}$ is run with a default stopping criterion set to $\absolutevalue*{\delta E^{(2)}/E^{(2)}}=0.002$ , where $\delta E^{(2)}$ is the statistical error associated with $E^{(2)}$ . We would like to stress that, thanks to the present semistochastic algorithm, the complete wave function is considered, and that no threshold is required. Consequently, size-consistency will be preserved if a size-consistent perturbation theory is applied.

While performing production runs, we have noticed that the computation of $E^{(2)}$ was faster than the CIPSI selection. Hence, we have slightly modified the routines computing $E^{(2)}$ such that the selection of determinants is performed alongside the computation of $E^{(2)}$ . This new on-the-fly CIPSI selection performed during the stochastic PT2 calculation completely removes the conventional (deterministic) selection step, and the determinants are selected with no additional cost. We have observed that, numerically, the curves of the variational energy as a function of $N_{\text{det}}$ obtained with either the deterministic or the stochastic selections are indistinguishable, so that the stochastic algorithm does not harm the selection’s quality.

For the selection of multiple states, one PT2 calculation is run for each state and, as proposed by Angeli et al., Angeli et al. (1997) the selection criterion is modified as

[TABLE]

with

[TABLE]

This choice gives a balanced selection between states of different multi-configurational nature.

IV spin-adapted wave functions

Determinant-based sCI algorithms generate wave functions expressed in a truncated space of determinants. Obviously, the selection presented in the previous section does not enforce that $\hat{H}$ commutes with $\hat{S}^{2}$ in the truncated space. Hence, the eigenstates of $\hat{H}$ are usually not eigenvectors of $\hat{S}^{2}$ , although the situation improves when the size of the internal space tends to be complete. A natural way to circumvent this problem is to work in the basis of configuration state functions (CSFs), but this representation makes the direct computation of the Hamiltonian less straightforward during the Davidson diagonalization.

Instead, we follow the same path as the MELD and SCIEL codes,Davidson et al. (1979); Kozlowski and Davidson (1994); Caballol et al. (1998) and identify all the spatial occupation patterns associated with the determinants.Applencourt et al. (2018) We then generate all associated spin-flipped configurations, and add to the internal space all the missing determinants. This procedure ensures that $\hat{H}$ commutes with $\hat{S}^{2}$ in the truncated space, and spin-adapted states are obtained by the diagonalization of $\hat{H}$ . In addition, we apply a penalty method in the diagonalization by modifying the Hamiltonian as Fales et al. (2017)

[TABLE]

where $\mathbf{I}$ is the identity matrix and $\gamma$ is a fixed parameter set to 0.1 by default. This improves the convergence to the desired spin state, but also separates degenerate states with different spins, a situation that can potentially occurs with Rydberg states. In the Davidson algorithm, this requires the additional computation of $\mathbf{S}^{2}\,\mathbf{U}$ , for which the cost is expected to be the same as the cost of $\mathbf{H}\,\mathbf{U}$ (see Algorithm 2). The cost of computing $\mathbf{H}\,\mathbf{U}$ and $\mathbf{S}^{2}\,\mathbf{U}$ is mostly due to the search of the connected pairs of determinants, namely the determinants $\bra{I}$ and $\ket{J}$ for which $\matrixelement*{I}{\hat{H}}{J}$ and $\matrixelement*{I}{\hat{S}^{2}}{J}$ are not zero due to Slater-Condon’s rules. We have modified the function computing $\mathbf{H}\,\mathbf{U}$ so that it also computes $\mathbf{S}^{2}\,\mathbf{U}$ at the same time. Hence, the search of connected pairs is done once for both operations and $\mathbf{S}^{2}\,\mathbf{U}$ is obtained with no extra computational cost.

Working with spin-adapted wave functions increases the size of the internal space by a factor usually between 2 and 3, but it is particularly important if one is willing to obtain excited states. Loos et al. (2018); Scemama et al. (2018, 2018); Loos et al. (2019) Therefore, the default in Quantum Package is to use spin-adapted wave functions.

V Parallelism

In Quantum Package, multiple parallelism layers are implemented: a fine-grained layer to benefit from shared memory, an intermediate layer to benefit from fast communication within a group of nodes, and a coarse-grained layer to interconnect multiple groups of nodes. Fine-grained parallelism is performed with OpenMPDagum and Menon (1998) in almost every single routine. Then, to go beyond a single compute node, Quantum Package does not use the usual single program/multiple data (SPMD) paradigm. A task-based parallelism framework is implemented with the ZeroMQ library.Hintjens (2013) The single-node instance runs a compute process as well as a task server process, while helper programs can be spawned asynchronously on different (heterogeneous) machines to run a distributed calculation. The helper programs can connect via ZeroMQ to the task server at any time, and contribute to a running calculation. As the ZeroMQ library does not take full advantage of the low latency hardware present in HPC facilities, the helper programs are parallelized also with the message passing interfaceForum (1994) (MPI) for fast communication among multiple client nodes, typically for fast broadcasting of large data structures.

Hence, we have 3 layers of parallelism in Quantum Package: OpenMP, MPI and ZeroMQ. This allows for an elastic management of resources: a running calculation taking too much time can be dynamically accelerated by plugging in more computing resources, by submitting more jobs in the queue or possibly in the cloud, i.e. outside of the HPC facility. This scheme has the advantage that it is not necessary to wait for all the nodes to be free to start a calculation, and hence minimizes the waiting time in the batch queue. It also gives the possibility to use altogether different helper programs. For instance, one could use a specific GPU-accelerated helper program on a GPU node while CPU-only helpers run on the CPU-only partition of the cluster. It is also possible to write a helper program that helps only one PT2/selection step and then exit, allowing to gather resources after the PT2/selection has started, and freeing them for the following diagonalization step.

The current limitation of Quantum Package is the memory of the single-node instance. We have not yet considered the possibility to add more compute nodes to increase the available memory, but this can be done by transforming the main program into an MPI program using scattered data structures.

We now describe how the Davidson and PT2/selection steps are parallelized.

V.0.1 Davidson diagonalization

In the direct Davidson diagonalization method, the computational bottleneck is the matrix product $\mathbf{W}=\mathbf{H}\,\mathbf{U}$ , and only this step needs to be distributed. The calculation is divided into independent tasks where each task builds a unique piece of $\mathbf{W}$ containing 40 000 consecutive determinants. Communicating the result of all the tasks scales as $\order*{N_{\text{det}}}$ , independently of the number of parallel processes. On the other hand, $\mathbf{U}$ needs to be broadcast efficiently at the beginning of the calculation to each slave process.

The computation of a task is parallelized with OpenMP, looping in a way that guarantees a safe write access to $\mathbf{W}$ , avoiding the need of a lock. When idle, a slave process requests a task to the ZeroMQ task server, computes the corresponding result and sends it to the collector thread of the master instance via ZeroMQ. As the OpenMP tasks are not guaranteed to be balanced, we have used a dynamic scheduling, with a chunk size of 64 elements. The reason for this chunk size is to force that multiple threads access to $\mathbf{W}$ at memory addresses far apart, avoiding the so-called false sharing performance degradation that occurs when multiple threads write simultaneously in the same cache line. Bolosky and Scott (1993) When the task is fully computed, the computed piece of $\mathbf{W}$ is sent back to the master process and a new task is requested, until the task queue is empty.

The $\mathbf{U}$ and $\mathbf{W}$ arrays are shared among threads, as well as all the large constant data needed for the calculation such as the ERIs. Sharing $\mathbf{U}$ also provides the benefit to reduce the amount of communication since $\mathbf{U}$ needs to be fetched only once for each node, independently of the number of cores. To make the broadcast of $\mathbf{U}$ efficient, the slave helper program is parallelized with MPI in a SPMD fashion, and each node runs a single MPI process. The $\mathbf{U}$ matrix is fetched from the ZeroMQ server by the process with rank zero, and then it is broadcast to the other slave processes within the same MPI job via MPI primitives. Then, each MPI process behaves independently and communicates via ZeroMQ with the task server, and with the master node which collects the results. A schematic view of the communication is presented in Fig 1.

V.0.2 CIPSI selection and PT2 energy

In the computation of $E^{(2)}$ and the CIPSI selection, each task corresponds to the computation of one $\varepsilon_{J}$ or $\varepsilon_{K}$ in Eq. (23), together with the selection of the associated external determinants. To establish the list of tasks, the Monte Carlo sampling is pre-computed on the master node. We associate to each task the number of drawn Monte Carlo samples such that running averages can be computed when the results of the tasks have been received by the collector thread. When the convergence criterion is reached, the task queue is emptied and the collector waits for all the running tasks to terminate.

As opposed to the Davidson implementation where each task is parallelized with OpenMP, here each OpenMP thread handles independently a task computed on a single core. Hence, there are multiple ZeroMQ clients per node, typically one per core, requesting tasks to the task server and sending the results back to the collector thread (see Fig. 2). Here, all the OpenMP threads are completely independent during the whole selection, and this explains the pleasing scaling properties of our implementation, as shown in Sec. VI.3. As in the Davidson distributed scheme, when the helper programs are run with MPI all the common data are communicated once from the ZeroMQ server to the rank-zero MPI process. Then, the data is broadcast to all the other processes with MPI primitives (there is one MPI process per node).

VI Results

VI.1 Capabilities of Quantum Package

Before illustrating the new features of Quantum Package in the next subsection. We propose to give an overview of what can be achieved (in terms of system and basis set sizes) with the current implementation of Quantum Package. To do so we propose to review some of our very recent studies.

In Ref. Loos et al., 2018, we studied 18 small molecules (water, hydrogen sulfide, ammonia, hydrogen chloride, dinitrogen, carbon monoxide, acetylene, ethylene, formaldehyde, methanimine, thioformaldehyde, acetaldehyde, cyclopropene, diazomethane, formamide, ketene, nitrosomethane, and the smallest streptocyanine) with sizes ranging from 1 to 3 non-hydrogen atoms. For such systems, using sCI expansions of several million determinants, we were able to compute more than hundred highly accurate vertical excitation energies with typically augmented triple- $\zeta$ basis sets. It allowed us to benchmark a series of 12 state-of-the-art excited-state wave function methods accounting for double and triple excitations.

Even more recently, Loos et al. (2019) we provided accurate reference excitation energies for transitions involving a substantial amount of double excitation using a series of increasingly large diffuse-containing atomic basis sets. Our set gathered 20 vertical transitions from 14 small- and medium-size molecules (acrolein, benzene, beryllium atom, butadiene, carbon dimer and trimer, ethylene, formaldehyde, glyoxal, hexatriene, nitrosomethane, nitroxyl, pyrazine, and tetrazine). For the smallest molecules, we were able to obtain well converged excitation energies with augmented quadruple- $\zeta$ basis set while only augmented double- $\zeta$ bases were manageable for the largest systems (such as acrolein, butadiene, hexatriene and benzene). Note that the largest sCI expansion considered in this study had more than 200 million determinants.

In Ref. Giner et al., 2018, Giner et al. studied even larger systems containing transition metals: \ce[CuCl4]^2-, \ce[Cu(NH3)4]^2+ and \ce[Cu(H2O)4]^2+. They were able, using large sCI expansions, to understand the physical phenomena that determine the relative energies of three of the lowest electronic states of each of these square-planar copper complexes.

VI.2 Extrapolation

To illustrate the extrapolation procedure described in Sec. II.4, we consider a cyanine dye Le Guennic and Jacquemin (2015) \ceH2N-CH=NH2+ (labeled as CN3 in the remaining) in both its ground state and first excited state. Send et al. (2011); Boulanger et al. (2014); Garniron et al. (2018) The geometry is the equilibrium geometry of the ground state optimized at the PBE0/cc-pVQZ level. Boulanger et al. (2014) The ground state is a closed shell, well described by a single reference, while the excited state is singly excited and requires, at least, two determinants to be properly modeled. The calculations were performed in the aug-cc-pVDZ basis set with state-averaged natural orbitals obtained from an initial CIPSI calculation. All the electrons were correlated, so the FCI space which is explored corresponds to a CAS(24,114) space. The reference excitation energy, obtained at the CC3/ANO-L-VQZP level is 7.18 eV Send et al. (2011) (see also Ref. Garniron et al., 2018). Note that this particular transition is fairly insensitive to the basis set as long as at least one set of diffuse functions is included. For example, we have obtained 7.14 and 7.13 eV at the CC3/aug-cc-pVDZ and CC3/aug-cc-pVTZ levels, respectively. Loos et al. (2018)

In Fig. 3, we plot the energy convergence of the ground state (GS) and the excited state (ES) as a function of the number of determinants $N_{\text{det}}$ , with and without the second-order perturbative contribution. From the data gathered in Table 1, one can see that, although $E^{(2)}$ is still large (roughly 0.02 a.u.), the sCI+PT2 and sCI+rPT2 excitation energies converge to a value of 7.20 eV compatible with the reference energy obtained in a larger basis set. We have also plotted the sCI+rPT2 energy given by $E^{(0)}+ZE^{(2)}$ (see Sec. II.4.2) and we clearly see that this quantity converges much faster than the usual sCI+PT2 energy. Even for very small reference wave function, the energy gap between GS and ES is qualitatively correct. The graph of Fig. 4, which shows the zeroth-order energy $E^{(0)}$ as a function of the second-order energy $E^{(2)}$ (dotted lines) or its renormalization variant $Z\,E^{(2)}$ (solid lines), also indicates that it is practically much easier to extrapolate to the FCI limit using the rPT2 correction.

As a second test case for rPT2, we consider the widely-studied example of the chromium dimer (\ceCr2) in its ${}^{1}\Sigma_{g}^{+}$ ground state. Scuseria and Schaefer III (1990); Roos and Andersson (1995); Brynda et al. (2009); Coe et al. (2014); Purwanto et al. (2015); Sokolov and Chan (2016); Sokolov et al. (2017); Tsuchimochi and Ten-no (2017); Li Manni et al. (2013); Vancoillie et al. (2016); Holmes et al. (2016); Guo et al. (2016); Sharma et al. (2017); Garniron et al. (2017) This system is notoriously challenging as it combines dynamic and static correlation effects hence requiring multi-configurational methods and large basis sets in order to have a balanced treatment of these two effects. Consequently, we compute its ground-state energy in the cc-pVQZ basis set with an internuclear distance $R_{\ce{Cr-Cr}}=1.68~{}\AA$ close to its experimental equilibrium geometry. Our full-valence calculation corresponds to an active space CAS(28,198) and the computational protocol is similar to the previous example. The second-order corrected value $E^{(0)}+E^{(2)}$ as well as its renormalized version $E^{(0)}+ZE^{(2)}$ as a function of the number of determinants in the reference wave function are reported in Table 2 and depicted in Fig. 5. Here also, we observe that rPT2 is clearly a superior extrapolation framework compared to the standard PT2 version as it yields a much straighter extrapolation curve, even in the case of a strongly correlated system such as \ceCr2. The renormalization factor $Z$ [see Eq. (17)] mitigates strongly the overestimation of the FCI energy for small wave functions by damping the second-order energy $E^{(2)}$ . Linear extrapolations of the PT2 and rPT2 energies based on the two largest wave functions yields extrapolated FCI energies of -2087.734 and -2087.738, respectively (see also Table 2). The difference between these two extrapolated FCI energies provides a qualitative idea of the extrapolation accuracy.

VI.3 Speedup

In this Section, we discuss the parallel efficiency of the algorithms implemented in Quantum Package. The system we chose for these numerical experiments is the benzene molecule \ceC6H6 for which we have performed sCI calculations with the 6-31G* basis set. The frozen-core approximation has been applied and the FCI space that we explore is a CAS(30,90). The measurements were made on GENCI’s Irene supercomputer. Each Irene’s node is a dual-socket Intel(R) Xeon(R) Platinum 8168 [email protected] with 192GiB of RAM, with a total of 48 physical CPU cores. Parallel speedup curves are made up to 12 288 cores (i.e. 256 nodes) for i) a single iteration of the Davidson diagonalization, and ii) the hybrid semistochastic computation of $E^{(2)}$ (which includes the CIPSI selection). The speedup reference corresponds to the single node calculation (48 cores).

First, we measure the time required to perform a single Davidson iteration as a function of the number of CPU cores for the two largest wave functions ( $N_{\text{det}}=25\times 10^{6}$ and $100\times 10^{6}$ ). The timings are reported in Table 3 while the parallel speedup curve is represented in Fig. 6. The parallel efficiency increases together with $N_{\text{det}}$ , as shown in Fig. 6. For the largest wave function, a parallel efficiency of 66% is obtained on 192 nodes (i.e. 9216 cores). We note that the speedup reaches a plateau at 3 072 cores (64 nodes) for $N_{\text{det}}=25\times 10^{6}$ . For this wave function, there are 625 tasks computing each 40 000 rows of $\mathbf{W}$ . When the number of nodes reaches 64, the number of tasks is too small for the load to be balanced between the nodes, and the computational time is limited by the time taken to compute the longest task. The same situation arises for $N_{\text{det}}=100\times 10^{6}$ with 9 408 cores (192 nodes), with 2 500 tasks to compute.

Second, we analyze the parallel efficiency of the calculation of $E^{(2)}$ for the sCI wave function with $N_{\text{det}}=25\times 10^{6}$ . The stopping criterion during the calculation of $E^{(2)}$ is given by a relative statistical error below $2\times 10^{-3}$ of the current $E^{(2)}$ value. The speedups are plotted in Fig. 6 (see also Table 3). For 192 nodes, one obtains a parallel efficiency of 89%. The present parallel efficiency is not as good as the one presented in the original paper.Garniron et al. (2017) The reason behind this is a faster computation of $e_{\alpha}^{(2)}$ , which reduces the parallel efficiency by increasing the ratio communication/computation.

VII Developing in Quantum Package

VII.1 The Quantum Package philosophy

Quantum Package is a standalone easy-to-use library for developers. The main goals of Quantum Package are to i) facilitate the development of new quantum chemistry methods, ii) minimize the dependency on external programs/libraries, and iii) encourage the collaborative and educative work through human readable programs. Therefore, from the developer point of view, Quantum Package can be seen as a standalone library containing all important quantities needed to perform quantum chemistry calculations, both involving wave function theory, through the determinant driven algorithms, and DFT methods, thanks to the presence of a quadrature grid for numerical integrations and basic functionals. These appealing features are made more concrete thanks to the organization of Quantum Package in terms of core modules and plugins (see Sec. VII.3) together with its programming language (see Sec. VII.2), which naturally creates a very modular environment for the programmer.

Although Quantum Package is able to perform all the required steps from the calculation of the one- and two-electron integrals to the computation of the sCI energy, interfacing Quantum Package, at any stage, with other programs is relatively simple. For example, canonical or CASSCF molecular orbitals can be imported from GAMESS, Schmidt et al. (1993) while atomic and/or molecular integrals can be read from text files like fcidump. Thanks to this flexibility, some of us are currently developing plugins for performing sCI calculations for periodic systems.

VII.2 The IRPF90 code generator

It is not a secret that large scientific codes written in Fortran (or in similar languages) are difficult to maintain. The program’s complexity originates from the inter-dependencies between the various entities of the code. As the variables are more and more coupled, the programs become more and more difficult to maintain and to debug. To keep a program under control, the programmer has to be aware of all the consequences of any source code modification within all possible execution paths. When the code is large and written by multiple developers, this becomes almost impossible. However, a computer can easily handle such a complexity by taking care of all the dependencies between the variables, in a way similar to how GNU Make handles the dependencies between source files.

IRPF90 is a Fortran code generator. Scemama (2009) Schematically, the programmer only writes computation kernels, and IRPF90 generates the glue code linking all these kernels together to produce the expected result, handling all relationships between variables. To illustrate in a few words how IRPF90 works, let us consider the simple example which consists of calculating the total energy of a molecular system as the sum of the nuclear repulsion and the electronic energy $E_{\text{tot}}=E_{\text{nuc}}+E_{\text{ele}}$ . The electronic energy is the sum of the kinetic and potential energies, i.e., $E_{\text{ele}}=E_{\text{kin}}+E_{\text{pot}}$ .

The production tree associated with the computation of the total energy is shown in Fig. 7. Within the IRPF90 framework, the programmer writes a provider for each entity, i.e., a node of the production tree. The provider is a subroutine whose only goal is to compute the value associated with the entity, assuming the values of the entities on which it depends are computed and valid. Hence, when an entity is used somewhere in the program (in a subroutine, a function or a provider), a call to its provider is inserted in the code before it is used such that the corresponding value is guaranteed to be valid.

Quantum Package is a library of providers designed to make the development of new wave function theory and DFT methods simple. Only a few programs using these providers are part of the core modules of Quantum Package, such as the sCI module using the CIPSI algorithm or the module containing the semi-stochastic implementation of the second-order perturbative correction. The main goal of Quantum Package is to be used as a library of providers, and programmers are encouraged to develop their own modules using Quantum Package.

VII.3 The plugin system

External programmers should not add their contributions by modifying directly Quantum Package’s core, but by creating their own modules in independent repositories hosted and distributed by themselves. This model gives more freedom to the developers to distribute modules as we do not enforce them to follow any rule. The developers are entirely responsible for their own plugins. This model has the advantage to redirect immediately the users to the right developer for questions, installation problems, bug reports, etc.

Quantum Package integrates commands to download external repositories and integrate all the plugins of these repositories into the current installation of Quantum Package. External plugins appear exactly as if they were part of Quantum Package, and if a plugin is useful for many users, it can be easily integrated in Quantum Package’s core after all the coding and documentation standards are respected.

Multiple external plugins were developed by the authors. For instance, one can find a multi-reference coupled cluster program,Giner et al. (2016); Garniron et al. (2017) interfaces with the quantum Monte Carlo programs QMC=Chem,Scemama et al. (2013) QMCPackKim et al. (2018) and CHAMP,Umrigar et al. (2018) an implementation of the shifted-Bk method,Garniron et al. (2018) a program combining CIPSI with RSDFT,Giner et al. (2018) a four-component relativistic RSDFT code, Paquier and Toulouse (2018) and many others.

In particular, Quantum Package also contains the basic tools to use and develop range-separated density-functional theory (RSDFT, see, e.g., Refs. Savin, 1996; Toulouse et al., 2004) which allows to perform multi-configurational density-functional theory (DFT) calculations within a rigorous mathematical framework. In the core modules of Quantum Package, single-determinant approximations of RSDFT are available, which fall into the so-called range-separated hybrid Gerber and Ángyán (2005); Ángyán et al. (2005) (RSH) approximation. These approaches correct for the wrong long-range behavior of the usual hybrid approximations thanks to the inclusion of the long-range part of the HF exchange. Quantum Package contains all necessary integrals to perform RSDFT calculations, including the long-range interaction integrals and Hartree-exchange-correlation energies and potentials derived from the short-range version of the local-density approximation (LDA) Paziani et al. (2006) and a short-range generalized-gradient approximation (GGA) based on the Perdew-Burke-Ernzerhof (PBE) functional. Goll et al. (2006) All numerical integrals are performed using the standard Becke quadrature grid Becke (1988) associated with the improved radial grids of Mura et *al.*Mura and Knowles (1996) With these tools, more evolved schemes based on RSDFT have been developed, such as an energy correlation functional with multideterminantal reference depending on the on-top pair density Ferté et al. (2019) or a basis set correction. Giner et al. (2018) The corresponding source code can be found as external plugins (see, for example, https://gitlab.com/eginer/qp_plugins_eginer).

VIII Conclusion

Significant improvements were brought to the second version of Quantum Package. Some were single-core optimizations, and others focused on the algorithm adaptation to large-scale parallelism (load balancing in particular). Currently, the code has a parallel efficiency that enables routinely to realize runs on roughly 2 000 CPU cores, with tens of millions of determinants in the reference space. Moreover, we have been able to push up to 12 288 cores (256 nodes) on GENCI’s supercomputer Irene. Such a gain in efficiency has and will lead to many more challenging chemical applications. Caffarel et al. (2014); Giner et al. (2013, 2015); Caffarel et al. (2016, 2016); Loos et al. (2018, 2019); Scemama et al. (2018); Dash et al. (2018); Pineda Flores and Neuscamman (2018)

The Davidson diagonalization, which is at the center of sCI and FCI methods, suffers from the impossibility to fully store the Hamiltonian in the memory of a single node. The solution we adopted was to resort to direct methods, i.e., recomputing on the fly the matrix elements at each iteration. While an extremely fast method was already available to detect zero matrix elements,Scemama and Giner (2013) the former implementation still had to search over the $\order*{N_{\text{det}}^{2}}$ matrix elements for interacting determinant pairs. Now, determinants are split in disjoint sets entirely disconnected from each other. Thus, only a small fraction of the matrix elements need to be explored, and an algorithm with $\order*{N_{\text{det}}^{3/2}}$ scaling was proposed. While the parallelization of this method was somewhat challenging due to the extremely unbalanced nature of the elementary tasks, a distributed implementation was realized with satisfying parallel speedups (typically 35 for 50 nodes) with respect to the 48-core single-node reference.

Significant improvements were also realized in the computation of the second-order perturbative correction, $E^{(2)}$ . A natural idea was to take into account the tremendous number of tiny contributions via a stochastic Monte Carlo approach. $E^{(2)}$ being itself an approximate quantity used for estimating the FCI energy, its exact value is indeed not required, as long as the value is unbiased and the statistical error is kept under control. Our scheme allows to compute $E^{(2)}$ with a small error bar for a few percent of the cost of the fully deterministic computation.

Similarly, the CIPSI selection is now performed stochastically alongside the PT2 calculation. Therefore, the selection part of the new stochastic CIPSI selection is virtually free as long as one is interested in the second-order perturbative correction.

Finally, efforts have been made to make this software as developer friendly as possible thanks to a very modular architecture that allows any developer to create his/her own module and to directly benefit from all pre-existing work.

License

Quantum Package is licensed under GNU Affero General Public License (AGPLv3).

Acknowledgements.

The authors would like to thank the Centre National de la Recherche Scientifique (CNRS) for funding and Cyrus Umrigar for carefully reading the manuscript. Funding from Projet International de Coopération Scientifique (PICS08310) is also acknowledged. This work was performed using HPC resources from CALMIP (Toulouse) under allocation 2019-18005 and from GENCI-TGCC (Grant 2018-A0040801738). A.B. was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, as part of the Computational Materials Sciences Program and Center for Predictive Simulation of Functional Materials. K.G. acknowledges support from grant number CHE1762337 from the U.S. National Science Foundation.

Appendix A Implementation details

A.1 Efficiency of integral storage

The efficiency of the storage as a hash table was measured on a dual socket Intel Xeon E5-2680 [email protected] processor, taking the water molecule with the cc-pVQZ basis set (115 MOs). The time to access all the integrals was measured by looping over the entire set of ERIs using different loop orderings. The results are given in Table 4, the reference being the storage as a plain four-dimensional array.

In the array storage, the value of 170 ns/integral in the random access case is typical of the latency to fetch a value in the RAM modules, telling that the requested integral is almost never present in any level of cache. When the data is accessed with a stride of one ( $i,j,l,k$ storage), the hierarchical architecture of the cache levels accelerates the access by a factor of 18, down to 9.71 ns/integral, corresponding mostly to the overhead of the function call, the retrieval of the data being negligible.

With the hash table, the random access is only 2.18 times slower than the random access in the array. Indeed, two random accesses are required: one for the first element of the key bucket to do the search, and one for the value of the integral. The remaining time corresponds to the binary search. The results show that data locality is exploited: when the access is done with a regular access pattern, the data is fetched roughly 3 times faster than using a random access, giving a latency below the latency of a random access in the array.

A CIPSI calculation was run once with the array storage, and once with the hash table storage. With the hash storage, the total wall clock time was increased only by a factor of two. To accelerate the access to the most frequently used integrals and reduce this overhead, we have implemented a software cache. All the integrals involving the 128 MOs closest to the Fermi level are copied in a dense array of 1284 elements (2 GiB), and benefit from the fastest possible access.

A.2 Internal representation of determinants

Determinants can be conveniently written as a string of creation operators applied to the vacuum state ${\ket{}}$ , e.g., $a^{\dagger}_{i}a^{\dagger}_{j}a^{\dagger}_{k}{\ket{}}=\ket{I}$ . Because of the fermionic nature of electrons, a permutation of two contiguous creation operators results in a sign change $a^{\dagger}_{j}a^{\dagger}_{i}=-a^{\dagger}_{i}a^{\dagger}_{j}$ , which makes their ordering relevant, e.g., $a^{\dagger}_{j}a^{\dagger}_{i}a^{\dagger}_{k}{\ket{}}=-\ket{I}$ . A determinant can be broken down into two pieces of information: i) a set of creation operators corresponding to the set of occupied spinorbitals in the determinant, and ii) an ordering of the creation operators responsible for the sign of the determinant, known as phase factor. Once an ordering operator ${\hat{\mathcal{O}}}$ is chosen and applied to all determinants, the phase factor may simply be included in the CI coefficient.

The determinants are built using the following order: i) spin-up ( $\uparrow$ ) spinorbitals are placed before spin-down ( $\downarrow$ ) spinorbitals, as in the Waller-Hartree double determinant representationPauncz (1989) ${\hat{\mathcal{O}}}\ket{I}=\hat{I}{\ket{}}=\hat{I}_{\uparrow}\hat{I}_{\downarrow}{\ket{}}$ , and ii) within each operator $\hat{I}_{\uparrow}$ and $\hat{I}_{\downarrow}$ , the creation operators are sorted by increasing indices. For instance, let us consider the determinant $\ket{J}=a^{\dagger}_{j}a^{\dagger}_{k}a^{\dagger}_{\bar{i}}a^{\dagger}_{i}{\ket{}}$ built from the set of spinorbitals $\{i_{\uparrow},j_{\uparrow},k_{\uparrow},i_{\downarrow}\}$ with $i<j<k$ . If we happen to encounter such a determinant, our choice of representation imposes to consider its re-ordered expression ${\hat{\mathcal{O}}}\ket{J}=-a^{\dagger}_{i}a^{\dagger}_{j}a^{\dagger}_{k}a^{\dagger}_{\bar{i}}{\ket{}}=-\ket{J}$ , and the phase factor must be handled.

The indices of the creation operators (or equivalently the spinorbital occupations), are stored using the so-called bitstring encoding. A bitstring is an array of bits; typically, the 64-bit binary representation of an integer is a bitstring of size 64. Quite simply, the idea is to map each spinorbital to a single bit with value set to its occupation number. In other words, 0 and 1 are associated with the unoccupied and occupied states, respectively.

For simplicity and performance considerations, the occupations of the spin-up and spin-down spinorbitals are stored in different bitstrings, rather than interleaved or otherwise merged in the same one. This allows to straightforwardly map orbital index $p$ to bit index $p-1$ (orbitals are usually indexed from 1, while bits are indexed from 0). This makes the representation of a determinant a tuple of two bitstrings, associated with respectively spin-up and spin-down orbitals. A similar parity representation of the fermionic operators is commonly used in quantum computing. Seeley et al. (2012)

The storage required for a single determinant is, in principle, one bit per spinorbital, or $2\times N_{\text{orb}}$ bits. However, because CPUs are designed to handle efficiently 64-bit integers, each spin part is stored as an array of 64-bit integers, the unused space being padded with zeros. The actual storage needed for a determinant is $2\times 64\times N_{\text{int}}$ bits, where $N_{\text{int}}=\left\lfloor(N_{\text{orb}}-1)/64\right\rfloor+1$ is the number of 64-bit integers needed to store one spin part.

Taking advantage of low-level hardware instructions, Scemama and Giner (2013) we are able, given two arbitrary determinants $\ket{I}$ and $\ket{J}$ , to find with a minimal cost the excitation operator $\hat{T}$ such that $\ket{J}=\hat{T}\ket{I}$ . This is a necessary step to obtain the $(i,j,k,l)$ indices of the two-electron integral(s) involved in the Hamiltonian matrix element between $\ket{I}$ and $\ket{J}$ . Then, fetching the values of the integrals can be done quickly using the hash table presented in Sec. III.1.

Because the data structure used to store determinants implies an ordering of the MOs, we also need to compute a phase factor. Here, we propose an algorithm to perform efficiently the computation of the phase factor. For a determinant $\ket{I}$ that is going to be used repeatedly for phase calculations, we introduce a phase mask represented as a bitstring:

[TABLE]

where $\wedge$ denotes the and bitwise operation, and $I[k]$ is the $k$ th bit of bitstring $I$ , corresponding to the $(k+1)$ th spinorbital of determinant $\ket{I}$ (remember that the orbital indices start at 1 and the bit indices start at 0). In other words, the $i$ th bit of the phase mask is set to 1 if the number of electrons occupying the $i+1$ lowest spinorbitals is odd, and 0 otherwise. When an electron of determinant $\ket{I}$ is excited from orbital $h$ to $p$ , the associated phase factor is

[TABLE]

where $\oplus$ denotes the exclusive or (xor) operation. So if the phase mask is available, the computation of the phase factor only takes a few CPU cycles. Another important aspect is to create efficiently the phase masks. We propose Algorithm 3, which computes it in a logarithmic time for groups of 64 MOs, taking advantage of the associativity of the exclusive or operator. Indeed, the “for” loop executes 6 cycles to update the mask for $2^{6}=64$ MOs.

A.3 Davidson diagonalization

Within Quantum Package, the Davidson diagonalization algorithm is implemented in its multi-state version. Algorithmically, the expensive part of the Davidson diagonalization is the computation of the matrix product $\mathbf{H}\,\mathbf{U}$ . As mentioned above (see Sec. II), two determinants $\ket{I}$ and $\ket{J}$ are connected via $\mathbf{H}$ (i.e. $\matrixelement*{I}{\hat{H}}{J}\neq 0$ ) only if they differ by no more than two spinorbitals. Therefore, the number of non-zero elements per row in $\mathbf{H}$ is equal to the number of single and double excitation operators, namely $\order*{N_{\uparrow}^{2}(N_{\text{orb}}-N_{\uparrow})^{2}}$ . As $\mathbf{H}$ is symmetric, the number of non-zero elements per column is identical. This makes $\mathbf{H}$ very sparse. However, for large basis sets, the whole matrix may still not fit in a single node memory, as the number of non-zero entries to be stored is of the order of $N_{\text{det}}N_{\uparrow}^{2}(N_{\text{orb}}-N_{\uparrow})^{2}$ . One possibility would be to distribute the storage of $\mathbf{H}$ among multiple compute nodes, and use a distributed library such as PBLASChoi et al. (1995) to perform the matrix-vector operations. Another approach is to use a so-called direct algorithm, where the matrix elements are computed on the fly, and this is the approach we have chosen in Quantum Package. This effectively means iterating over all pairs of determinants $\ket{I}$ and $\ket{J}$ , checking whether $\ket{I}$ and $\ket{J}$ are connected by $\mathbf{H}$ and if they are, accessing the corresponding integral(s) and computing the phase factor. Even though it is possible to compute the excitation degree between two determinants very efficiently,Scemama and Giner (2013) the number of such computations scales as $N_{\text{det}}^{2}$ , which becomes rapidly prohibitively high. To get an efficient determinant-driven implementation it is mandatory to filter out all pairs of determinants that are not connected by $\mathbf{H}$ , and iterate only over connected pairs. To reach this goal, we have implemented an algorithm similar to the Direct Selected Configuration Interaction Using Strings (DISCIUS) algorithm.Povill and Rubio (1995)

The determinants of the internal space are re-ordered in linear time as explained in Ref. Scemama et al., 2016, such that the wave function can be expressed as

[TABLE]

where we take advantage of the Waller-Hartree double determinant representation. Pauncz (1989)

Moving along a row or a column of $\mathbf{C}$ keeps the spin-up or spin-down determinants fixed, respectively. For a given determinant, finding the entire list of same-spin single and double excitations can be performed in $\order*{{N_{\text{det}}^{\uparrow}}}=\order*{{N_{\text{det}}^{\downarrow}}}=\order*{\sqrt{N_{\text{det}}}}$ , while finding the opposite-spin double excitations is done via a two-step procedure. First, we look for all the spin-up single excitations. Then, starting from this list of spin-up single excitations, we search for the spin-down single excitation such that the resulting opposite-spin doubly-excited determinant belongs to $\Psi^{(0)}$ . Hence, the formal scaling is reduced to $\order*{N_{\text{det}}^{3/2}}$ . It could be further reduced to $\order*{N_{\text{det}}}$ at the cost of storing the list of all singly- and doubly-excited determinants for each spin-up and spin-down determinant, but we preferred not to follow this path in order to reduce the memory footprint as much as possible.

A.4 CIPSI selection and PT2 energy

There are multiple ways to compute the $e_{\alpha}^{(2)}$ ’s. One way is to loop over pairs of internal determinants $\ket{I}$ and $\ket{J}$ , generate the list of external determinants $\quantity{{\ket{\alpha}}}$ connecting $\ket{I}$ and $\ket{J}$ and increment the corresponding values $e_{\alpha}^{(2)}$ stored in a hash table. Using a hash table to store in memory a list of ${\ket{\alpha}}$ ’s without duplicates and their contributions $e_{\alpha}^{(2)}$ is obviously not a reasonable choice since the total number of ${\ket{\alpha}}$ ’s scales as $\order*{N_{\text{det}}N_{\uparrow}^{2}\quantity(N_{\text{orb}}-N_{\uparrow})^{2}}$ . To keep the memory growth in check, we must design a function that can build a stream of unique external determinants, compute their contribution $e_{\alpha}^{(2)}$ and retain in memory only the few most significant pairs $({\ket{\alpha}},e_{\alpha}^{(2)})$ .

In Quantum Package, we build the stream of unique external determinants as follows. We loop over the list of internal determinants (the generators) sorted by decreasing $c_{I}^{2}$ . For each generator $\ket{I}$ , we generate all the singly- and doubly-excited determinants $\quantity{{\ket{\alpha}}}$ , removing from this set the internal determinants and the determinants connected to any other generator $\ket{J}$ such that $J<I$ . This guarantees that the ${\ket{\alpha}}$ ’s are considered only once, without any additional memory requirement.

For each generator $\ket{I}$ , before generating its set of ${\ket{\alpha}}$ ’s, we pre-compute the diagonal of the Fock matrix associated with $\ket{I}$ . This enables to compute the diagonal elements $\matrixelement*{\alpha}{\hat{H}}{\alpha}$ involved in Eq. (8) for a few flops.Cimiraglia (1996) The computation of $\matrixelement*{\Psi^{(0)}}{\hat{H}}{\alpha}=\sum_{J}c_{J}\matrixelement*{J}{\hat{H}}{\alpha}$ is more challenging than the diagonal term since, at first sight, it appears to involve the $N_{\text{det}}$ internal determinants. Fortunately, most of the terms amongst this sum vanish due to Slater-Condon’s rules. Indeed, we know that the terms where $\ket{J}$ is more than doubly excited with respect to ${\ket{\alpha}}$ vanish, and these correspond to the determinants $\ket{J}$ which are more than quadruply excited with respect to $\ket{I}$ . Cimiraglia (1996) To compute efficiently $\matrixelement*{\Psi^{(0)}}{\hat{H}}{\alpha}$ , for each generator $\ket{I}$ , we create a filtered wave function $\ket*{\Psi^{(0)}_{I}}$ by projecting $\ket*{\Psi^{(0)}}$ on a subset $\mathcal{J}_{I}$ of internal determinants $\{\ket{J}\}$ where $\matrixelement*{J}{\hat{H}}{\alpha}$ is possibly non-zero. This yields $\matrixelement*{\Psi^{(0)}}{\hat{H}}{\alpha}=\matrixelement*{\Psi^{(0)}_{I}}{\hat{H}}{\alpha}$ , where $\Psi^{(0)}_{I}$ is a much smaller determinant expansion than $\Psi^{(0)}$ . In addition, as we have defined the ${\ket{\alpha}}$ ’s in such a way that they do not interact with $\ket{J}$ when $J<I$ , all these $\ket{J}$ ’s can also be excluded from $\mathcal{J}_{I}$ . This pruning process yielding to $\ket*{\Psi^{(0)}_{I}}$ will be referred to as the coarse-grained filtering. A fine-grained filtering of $\ket*{\Psi^{(0)}_{I}}$ is performed in a second stage to reduce even more the number of determinants, as we shall explain later.

To make the coarse-grained filtering efficient, we first filter out the determinants that are more than quadruply excited in the spin-up and spin-down sectors separately. Using the representation shown in Eq. (29), this filtering does not need to run through all the internal determinants and scales as $\order*{N_{\text{det}}^{\uparrow}}=\order*{\sqrt{N_{\text{det}}}}$ . It is important to notice that, at this stage, the size of $\mathcal{J}_{I}$ is bounded by the number of possible quadruple excitations in both spin sectors, and does not scale any more as $\order*{N_{\text{det}}}$ . Next, we remove the determinants that are i) quadruply excited in one spin sector and excited in the other spin sector, ii) triply excited in one spin sector and more than singly excited in the other spin sector, and iii) all the determinants that are doubly excited in one spin sector and more than doubly excited in the other spin sector.

The external determinant contributions are computed in batches. A batch $I_{pq}$ is defined by a doubly-ionized generator $\ket{I_{pq}}=a_{p}a_{q}\ket{I}$ . When a batch is created, the fine-grained filtering step is applied to $\mathcal{J}_{I}$ to produce $\mathcal{J}_{I_{pq}}$ and $\Psi^{(0)}_{I_{pq}}$ , such that $\matrixelement*{\Psi^{(0)}_{I_{pq}}}{\hat{H}}{\alpha}=\matrixelement*{\Psi^{(0)}_{I}}{\hat{H}}{\alpha}$ .

Each external determinant produced in the batch $I_{pq}$ is characterized by two indices $r$ and $s$ with ${\hat{\mathcal{O}}}a^{\dagger}_{r}a^{\dagger}_{s}a_{p}a_{q}\ket*{I}=\ket*{I^{rs}_{pq}}$ . The contribution associated with each determinant of a given batch will be computed incrementally in a two-dimensional array $A(r,s)$ as follows. A first loop is performed over all the determinants $\ket{J}$ belonging to the filtered internal space $\mathcal{J}_{I_{pq}}$ . Comparing $\ket{J}$ to $\ket{I_{pq}}$ allows to quickly identify if $\ket{J}$ will be present in the list of external determinants, and consequently tag the corresponding cell $A(r,s)$ as banned. Banned cells will not be considered for the computation of $e_{\alpha}^{(2)}$ nor the determinant selection, as they correspond to determinants already belonging to the internal space. A second loop over all the $\ket{J}\in\mathcal{J}_{I_{pq}}$ is then performed. During this loop, all the $(r,s)$ pairs where $\ket*{I^{rs}_{pq}}$ is connected to $\ket{J}$ are generated, and the corresponding cells $A(r,s)$ are incremented with $c_{J}\matrixelement*{J}{\hat{H}}{I^{rs}_{pq}}$ . After this second loop, $A(r,s)=\matrixelement*{\Psi}{\hat{H}}{I^{rs}_{pq}}$ and all the contributions $e_{\alpha}^{(2)}$ of the batch can be obtained using $A(r,s)$ . The running value of $E^{(2)}$ is then incremented, and the $N_{\text{det}}$ most significant determinants are kept in an array sorted by decreasing $\absolutevalue*{e_{\alpha}^{(2)}}$ .

Figure 8 shows the number of determinants retained in $\Psi^{(0)}_{I}$ or $\Psi^{(0)}_{I_{pq}}$ after filtering out disconnected determinants of the ground state of the CN3 molecule with 935 522 determinants (see Sec. A). This example shows that, starting from $\Psi^{(0)}$ , the coarse-grained process which consists of removing the determinants more than quadruply excited with respect to the generator $\ket{I}$ produces wave functions $\Psi^{(0)}_{I}$ with a typical size of 120 000 determinants, a reduction by a factor 8. Then, starting from $\Psi^{(0)}_{I}$ , the fine-grained filtering, specific to the batch generating $\Psi^{(0)}_{I_{pq}}$ , reduces even more the number of determinants (by a factor 3), down to a typical size of 40 000 determinants, which represents only $4\%$ of the total wave function $\Psi^{(0)}$ .

Bibliography155

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Moore (1965) Moore, G. Cramming More Components onto Integrated Circuits. Electronics 1965 , 38 , 114–117.
2top (2018) Lists | | | TOP 500 Supercomputer Sites. 2018; https://www.top 500.org/lists/1993/11 , [Online; accessed 9. Oct. 2018].
3Sutter and Larus (2005) Sutter, H.; Larus, J. Software and the concurrency revolution. Queue 2005 , 3 , 54.
4Wulf and Mc Kee (1995) Wulf, Wm. A.; Mc Kee, S. A. Hitting the memory wall: implications of the obvious. SIGARCH Comput. Archit. News 1995 , 23 , 20–24.
5Khan et al. (2018) Khan, H. N.; Hounshell, D. A.; Fuchs, E. R. H. Science and research policy at the end of Moore’s law. Nature Electronics 2018 , 1 , 14–21.
6Booth et al. (2009) Booth, G. H.; Thom, A. J. W.; Alavi, A. Fermion Monte Carlo without fixed nodes: A game of life, death, and annihilation in Slater determinant space. J. Chem. Phys. 2009 , 131 , 054106.
7Booth and Alavi (2010) Booth, G. H.; Alavi, A. Approaching chemical accuracy using full configuration-interaction quantum Monte Carlo: A study of ionization potentials. J. Chem. Phys. 2010 , 132 , 174104.
8Cleland et al. (2010) Cleland, D.; Booth, G. H.; Alavi, A. Communications: Survival of the fittest: Accelerating convergence in full configuration-interaction quantum Monte Carlo. J. Chem. Phys. 2010 , 132 , 041103.