Highly Efficient and Scalable Compound Decomposition of Two

Subscriber access provided by University of Texas Libraries

Article

Highly efficient and scalable compound decomposition of two-electron integral tensor and its application in coupled cluster calculations Bo Peng, and Karol Kowalski J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00605 • Publication Date (Web): 23 Aug 2017 Downloaded from http://pubs.acs.org on August 25, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Highly efficient and scalable compound decomposition of two-electron integral tensor and its application in coupled cluster calculations Bo Peng∗ and Karol Kowalski∗ William R. Wiley Environmental Molecular Sciences Laboratory, Battelle, Pacific Northwest National Laboratory, K8-91, P.O. Box 999, Richland, WA 99352, USA E-mail: [email protected]; [email protected]

∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract The representation and storage of two-electron integral tensors are vital in largescale applications of accurate electronic structure methods. Low-rank representation and efficient storage strategy of integral tensors can significantly reduce the numerical overhead and consequently time-to-solution of these methods. In this paper, by combining pivoted incomplete Cholesky decomposition (CD) with a follow-up truncated singular vector decomposition (SVD), we develop a decomposition strategy to approximately represent the two-electron integral tensor in terms of low-rank vectors. A systematic benchmark test on a series of 1-D, 2-D, and 3-D carbon-hydrogen systems demonstrates high efficiency and scalability of the compound two-step decomposition of the two-electron integral tensor in our implementation. For the size of atomic basis set Nb ranging from ∼ 100 up to ∼ 2, 000, the observed numerical scaling of our implementation shows O(Nb2.5∼3 ) versus O(Nb3∼4 ) cost of performing single CD on the two-electron integral tensor in most of other implementations. More importantly, this decomposition strategy can significantly reduce the storage requirement of the atomic-orbital (AO) two-electron integral tensor from O(Nb4 ) to O(Nb2 log10 (Nb )) with moderate decomposition thresholds. The accuracy tests have been performed using ground- and excited-state formulations of coupled-cluster formalism employing single and double excitations (CCSD) on several benchmark systems including the C60 molecule described by nearly 1,400 basis functions. The results show that the decomposition thresholds can be generally set to 10−4 to 10−3 to give acceptable compromise between efficiency and accuracy.

Introduction Accurate electronic structure calculations often involve expensive tensor contractions. A typical example is the standard coupled-cluster (CC) theory, 1–8 of which the wave function of the system of interest is written as an exponential ansatz, eT |Φi, with T the cluster operator and |Φi the reference wave function. Higher excitations can be included in the T in an iterative 2


Page 2 of 48

Page 3 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


or perturbative manner to provide a systematically improved hierarchy of approximations. For this reason CC formalism has evolved into a method-of-choice for accurate predictions of geometrical structures, reaction dynamics, molecular properties, and excited-state processes (see Ref. 8–16 for recent reviews). However, the steep polynomial scaling and high storage requirements originating from complex contractions of high-dimensional tensors (e.g. two-electron integral tensor) are well-known bottlenecks that preclude CC methods from being applied to large systems in various areas of interest including molecular properties and excitation energies. 14–22 For example, when employing rudimentary CCSD approach (CC approach with singles and doubles) 5 the underlying tensor contractions can scale as O(Nb6 ) (Nb is the basis functions representing the size of a quantum chemical system). When performing excited-state calculations, such as equation-of-motion CCSD (EOM-CCSD) 23–25 and linear-response CCSD (LR-CCSD), 26,27 this polynomial scaling will be further increased by a significant prefactor. It has been shown that by properly selecting intermediate arrays and optimizing the O(Nb6 ) and O(Nb5 ) loops, the efficiency of a CCSD code can be increased by a factor of 5. 6 To a larger extent, high numerical costs associated with the polynomial scaling can be effectively addressed by the development of highly scalable implementations of CC methods, as evidenced by several recent benchmark calculations. 28–39 Growing interest in efficient utilization of peta- and soon-to-be exa-scale computational resources has stimulated an intensive development of various tensor libraries 40–53 that can be exploited in generating scalable CC codes for homogeneous as well as for many/multi-core computer systems. 54–61 Nevertheless, in all above mentioned examples of canonical CC implementations the storage requirement will quickly grow as a function of the system size to become a storage and communication bottleneck when going from mid- (102 ∼ 103 basis functions) to large-scale (103 ∼ 104 basis functions) CC calculations. Although it has been shown that by employing integral-direct algorithms the storage requirement can be greatly minimized, the integral-direct way might also bring frequent I/O operations and/or the necessity of recalculating “on-the-fly” atomic two-electron integrals, which would then increase the CPU

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

time and deteriorate the scaling with system size. Recent progress has been achieved in the development of low-order (or even linear) scaling techniques in the electronic structure calculations. Proposed techniques such as the fast multipole method (FMM), 62–64 fast quadratures, 65 and conjugate gradient density matrix search (CGDMS), 66 discontinuous Galerkin approach with adaptive local basis (DG-ALB) 67–70 enable linear-scaling density functional theory (DFT) calculations. In the CC regime, local approximations for the ground-state CC formulations significantly ameliorates these problems by primarily reducing the excitation manifold and in consequence the number of interactions involved in the calculations of ground-state CC energy. 71–85 The techniques such as local pair natural orbitals (LPNO) or domain-based local pair natural orbital (DLPNO) formulation 76,78,80,85 are capable of reproducing CC correlation energies obtained in canonical calculations. However, successful ground-state localization techniques may not be directly transferred to calculations of multiple excited states and properties (as shown in Ref. 22, 86, and 87, localization procedures may result in sizable errors in procedures of calculating dipole polarizabilities). From this perspective, CC formulations that first take advantage of controllable approximations to electron-electron interactions (and possible re-factorization of CC equations stemming from these approximations) prior to the redundant amplitudes elimination may be an alternative approach to properties and excited-state CC calculations. To reduce the high storage requirement of the two-electron integral tensor, various methods for low-rank representation of multi-dimensional tensors can be exploited (see Ref. 88 for a recent review). In electronic structure calculations, well-established density-fitting (DF) 89–97 and Cholesky decomposition (CD) 98–102 methods are often adapted to represent the four-index two-electron integral tensor in terms of products of three-index tensors. These methods have been applied to many molecular-orbital (MO) based methods to reduce the computational cost. 91,103–115 For example, DF is now the standard tool to reduce the computational cost of the second order Møller-Plesset perturbation calculations for large molecular systems and periodic systems. 91,92,116–125 In constrast to DF, the CD method provides a

4


Page 4 of 48

Page 5 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


means to control the accuracy of the decomposition to any arbitrary precision, and might be preferred in accurate electronic structure calculations. Since the length of each Cholesky vector is usually a quadratic function of Nb , even if the number of Cholesky vectors scales linearly as Nb , 126 the storage requirement of using Cholesky vectors (∼ O(Nb3 )) can still lead to a bottleneck in large-scale electronic structure calculations (e.g. for 10, 000 basis functions, the storage of Cholesky vectors could be as large as 15 ∼ 20 TB). In order to further reduce the storage requirement to avoid the bottleneck, Khoromskaia et al. 127 recently reported that by combining CD with the quantized tensor train (QTT) approximation of long Cholesky vectors the total storage requirement can be reduced by a factor of 10 for some moderate size systems. More recently, inspired by matrix bandwidth reduction algorithms in graph theory and their preliminary applications in quantum chemistry simulations, 128–132 we proposed to compute the optimal reordering of the original tensor prior to performing pivoted incomplete CD. 133 It turned out that the bandwidth reduced form of the original integral tensor is quite often a block-diagonal form, which can be easily implemented in a parallel manner to reduce the space and time required by the conventional incomplete CD procedure. In this paper, by combining pivoted incomplete CD with a follow-up truncated singular vector decomposition (SVD), we propose a compound decomposition strategy to approximately represent the two-electron integral tensor in terms of low-rank vectors. By carrying out systematic benchmark tests on a series of 1-D, 2-D, and 3-D carbon-hydrogen systems, we show that our implementation of this compound decomposition is highly efficient and scalable, and offers a significant reduction in the storage requirement. By further performing the ground- and excited-state CCSD calculations on these systems as well as a delocalized C60 molecule, we show that moderate truncation thresholds can be chosen for both CD and SVD procedures to give acceptable compromise between efficiency and accuracy.

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 6 of 48

Page 7 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The element of the Coulomb-type integral tensor Jµνλσ is defined as

Jµνλσ = (µν|λσ) =

Z Z

φµ (r1 )φν (r1 )

1 φλ (r2 )φσ (r2 ) dr1 dr2 , kr1 − r2 k

(1)

where φµ ’s (µ = 1, . . . , Nb ) are atomic orbitals (AO) or basis functions, subscripts µ, ν, λ, σ are AO indices, and r1 and r2 are electron coordinates. In practice, due to the positive semidefiniteness of J, its corresponding CD can be done through pivoting in an incomplete manner J ≈ LLT ,

(2)

where L contains the Cholesky vectors. The Frobenius norm of the residual associated with this approximation r = kJ − LLT k

(3)

can be controlled by picking a corresponding threshold θCD proportional to arbitrary precision ǫ for the largest diagonal element of the residual matrix. In such way, when the largest diagonal element in the residual matrix is smaller than θCD , the procedure will be terminated and r will be confined to be within a range proportional to ǫ. Obtained from the pivoted incomplete CD, L = {L1 . . . Lk . . . Lm } is a rank-m matrix with m the number of Cholesky vectors. For large systems, L is usually sparse, which is thus expected to be further compressed. To do that, we opt to do truncated singular vector decomposition (SVD) for each single Cholesky vector. Given an arbitrary Cholesky vector Lk = {Lk1 . . . Lkµν . . . LkN 2 } from L, its SVD gives b

Lk = Uk Λk (Vk )T ,

(4)

where Uk and Vk are matrices containing orthonormal vectors, and Λk is a diagonal matrix with Nb diagonal entries sorted in descending order, |σ1 | ≥ . . . ≥ |σn | ≥ . . . ≥ |σNb |, indicating the importance of the corresponding vectors in Uk and Vk . From this perspective, 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 48

Lk can be approximated by a rank-n tensor Lk(n) (n < Nb ), which is composed of n “most important” vectors from Uk and Vk (denoted as Uk(n) and Vk(n) ), and the first n singular values from Λk (denoted as Λk(n) ),

Lk ≈ Lk(n) = Uk(n) Λk(n) (Vk(n) )T .

(5)

A proper n needs to satisfy the following inequality q q 2 2 σ12 + . . . + σn2 ≥ σn+1 + . . . + σN , b

(6)

and a proper range of the SVD truncation threshold θSVD is then given by |σn+1 | < θSVD < |σn |. The right hand side of (6) defines the Frobenius norm of the error vector introduced by the truncated SVD of Lk ,

kLk − Lk(n) k =

q

2 2 σn+1 + . . . + σN . b

(7)

Since Lk is usually real symmetric (so is Lk(n) ), the SVD corresponds to eigen-decomposition, and Uk(n) and Vk(n) are identical. Then the storage requirement of the right hand side of Eq. (5) is only O(nNb + n), in contrast to O(Nb2 ) space required for storing Lk . Fig. 1 shows a simple example of this low-rank CD-SVD procedure working on the Coulomb-type two-electron integral tensor of a methane molecule in 6-31+G(d) basis, where both θCD and θSVD were set to 10−6 . After the compound decomposition, in order to obtain the integral tensor in molecular orbital (MO) space, the singular vectors (i.e. Uk(n) ) can be first contracted with the MO coefficient matrix C,

k(n) U pδ

=

Nb X

k(n)

Cµp Uµδ ,

µ=1

8


(8)

Page 9 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where Uk(n) denotes the transformed singular vector and p, as well as q, r, and s, are MO indices. Then, the MO integral tensor (pq|rs), when needed, can be approximated through

(pq|rs) ≈

m X

Lkpq Lkrs , with Lkpq =

k=1

n X

k(n)

k(n)

k(n)

U pδ Λδδ U qδ

.

(9)

δ=1

The generated MO integral tensor will be used in ground- and excited-state CC calculations. In this paper we have explored the accuracy of approximate CC formalisms utilizing the approximation of the two electron integral tensor obtained from low-rank SVD vectors in calculating ground-state energies, excitation energies, and properties (dipole polarizabilities). All tests were performed using singles and doubles approximations, i.e., CCSD, EOM-CCSD, and LR-CCSD formalisms. All CC formulations mentioned earlier exploit CCSD ground-state parametrization, where the electronic wave function |ΨCC i is given by the exponential ansatz |ΨCC i = eT1 +T2 |Φi ,

(10)

where |Φi is the so-called reference function (usually a Hartree-Fock Slater determinant). T1 and T2 represent singly and doubly excited cluster operators, respectively, which take the form

T1 =

X

tia a†a ai ,

(11)

i,a

T2 =

1 X ab † † t a a aj ai . 4 i,j,a,b ij a b

(12)

Here tai and tij ab are cluster amplitudes (indices i, j, . . . (a, b, . . .) stand for occupied (unoccupied) spinorbitals) and operators ap (a†p ) correspond to annihilation (creation) operators. Cluster operators T1 and T2 can iteratively be determined from the energy independent

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 48

coupled cluster equations hΦai | HN eT1 +T2 T1 +T2 hΦab ij | HN e

C

|Φi = 0 ,

(13)

C

|Φi = 0 ,

(14)

where subscript “C” designates a connected part of a given operator expression and HN represents electronic Hamiltonian H in a normal product form, i.e. HN = H − hΦ|H|Φi. Having obtained T1 and T2 operators, the CCSD correlation energy can be obtained from the expression ∆E = hΦ|(HN eT1 +T2 )C |Φi .

(15)

Once the ground-state CCSD problem is solved, the excited-state wave for n-th state in the EOM-CCSD parametrization takes the following form, |Ψn i = Rn eT1 +T2 |Φi,

(16)

where the excitation operator Rn in the EOM-CCSD approach is defined as

Rn = Rn,0 + Rn,1 + Rn,2

where Rn,0 is a scalar operator and Rn,1 =

X

ria (n)a†a ai and Rn,2 =

i,a

(17) 1 X ab r (n)a†a a†b aj ai . 4 i,j,a,b ij

The energy of the n-th excited state En can be obtained by solving eigenvalue problem ¯ n |Φi = En Rn |Φi, HR

(18)

¯ is given by H ¯ = e−T HeT . where similarity transformed Hamiltonian H The static and frequency dependent dipole polarizability can be calculated using LRCCSD formulation. Generally, employing the LR-CCSD formulation, the second-order prop-

10


Page 11 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


erties for a given frequency ω are given by i i hh i o n h 1 ¯ T (1) , T (1) |Φi , ¯ T (1) |Φi + hΦ|(1 + Λ) H, hhA; Biiω = Cˆ ±ω PÂ,B hΦ|(1 + Λ) A, A,ω B,ω B,−ω 2 (19) where Cˆ enforces time-reversal symmetry and Pˆ permutes operators A and B. Evaluating this quantity requires the evaluation of the Λ amplitudes (Lagrange multipliers) of gradient theory and the first-order response with respect to operators A and B at both +ω and −ω frequencies. For the calculation of dipole polarizabilities, A and B are replaced by dipole moment operators (denoted as d), and the dipole polarizability is given by

αIJ (ω) = −hhdI ; dJ iiω , I, J = x, y, z.

(20)

Results and discussion The low-rank compound decomposition, or CD-SVD, of the two-electron integral tensor in AO space has been implemented in the NWChem quantum chemistry platform 30 using the Global Array (GA) 134,135 shared-memory programming tools. The GA provides a convenient global-shared view for the two-electron integral tensor (a multi-dimensional array), and a complete functionalities for parallel code development.

Scalability test The original algorithm to carry out the incomplete pivoted CD was taken from Ref. 99, 126, and 136, which can be summarized in Algorithm 1.

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 48

Algorithm 1: Pivoted incomplete Cholesky algorithm for Coulomb-type integral tensor 1:

procedure IncompleteCholesky(m, L)

2:

m←0

3:

m ← Jλσλσ Dλσ

4: 5: 6: 7:

⊲ Get diagonals of the integral tensor

m m λσ ← λσ max Dλσ

m

that gives the maximal diagonal

m do while Dλσ m > θCD if m > 0 then

m m Rµνλσ m ← J µνλσ −

8:

end if

9:

m Lm µν ← Rµνλσ m /

q

Pm−1 i=0

m m Dλσ ← Dλσ − Lm Lm λσ λσ

11:

m←m+1 m m λσ ← λσ max Dλσ

13:

end while

14:

return m, L

15:

end procedure

i Liµν Lλσ m

m Dλσ m

10:

12:

⊲ Find λσ

⊲ Compute Cholesky vector ⊲ Update diagonals

⊲ Find λσ

m

that gives the maximal diagonal

As can be seen, this algorithm does not require the calculation and storage of the entire integral tensor. At each iteration, it only needs to calculate new O(Nb2 ) integrals, which can be easily parallelized since there is no communication between the computation of each integral. The generated Cholesky tensor, depending on its size, can be either on hard disk or in memory for the following truncated SVD decomposition. Overall, as pointed out in

12


Page 13 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


previous studies, 99,126,136 the full complexity of the incomplete CD algorithm is O(m2 Nb2 ), and the total storage requirement is O(mNb2 ). It should be pointed out that the storage of the Cholesky tensor is temporary, and distributed over available processors for the following truncated SVD. For each Choleksy vector, once the truncated SVD is done, the corresponding storage space will be released right away. Beside the parallelization of the computation of two-electron integrals at each iteration in the CD step, the algorithm can be further parallelized by explicitly controlling data locality and granularity. In our current implementation, this is achieved by thoroughly applying the GA functionalities in the following parts, (1) the GA distribution functionality is utilized to distribute all the vectors only according to the compound index (e.g. µν and λσ), and (2) the GA get-compute-put parallel model is used in step 7, 9,10, and the following truncated SVD decomposition. Note that in the get-compute-put model, part of the shared vector will be first copied to the local memory of each processor for parallel computation, and when finished the result will be copied back to the shared vector. A systematic scalability test has been done for a series of carbon-hydrogen systems, whose structures are shown in Fig. 2. As can be seen, starting from a simple methane molecule, three system-growth patterns were chosen to increase the system size based on the number of growing dimensions, namely, 1-D growing, 2-D growing, and 3-D growing. For the sake of accelerating the following electron correlation calculations, all structures are constrained to the corresponding highest abelian group, i.e. 1-D and 2-D structures are constrained to the C2h symmetry group, while 3-D structures are constrained to the C2v symmetry group. Only C-C single bond and C-H bond are included in all the structures shown in Fig. 2. Unless stated otherwise, the 6-31+G(d) basis set was used in most tests, which consists of 18 basis functions localized on each carbon atom, and 2 basis functions localized on each hydrogen atom. The basis set used throughout was not selected to provide highly accurate results but rather to allow us to consider a wide range of system size. For the 6-31+G(d) basis set, there are totally 1984, 1972, and 1806 basis functions generated for the largest 1-D, 2-D, and 3-D

13




Page 14 of 48

Page 15 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


structure, respectively. All the tests were performed on Constance supercomputer located in the Pacific Northwest National Laboratory’s Institutional Computing (PIC) facility. The Constance supercomputer comprises about 520 computer nodes with 56 Gb/s Fourteen Data Rate (FDR) Infiniband interconnect. Each node is a Dual Intel Haswell E5-2670 CPU (2.3 GHz) giving 24 cores, and each node is equipped with 64 GB 2133 MHz DDR4 memory. The summary of the scalability test can be seen in Fig. 3, which shows the change of the wall time of the CD-SVD (denoted as TWall ) for the studied systems as a function of the number of cores (denoted as Ncore ) and the number of basis functions (i.e. Nb ). Fig. 3a shows the accelaration curves of the CD-SVD for the three largest carbon-hydrogen systems when increasing Ncore from 288 to 384 and 576. The linear speedup curves are also given (dashed lines) for comparison. As can be seen, both C90 H182 (1-D) and C96 H122 (2-D) show superlinear speedup when comparing timings for Ncore = 288 and Ncore = 384. Here, we attribute the superlinear speedup to the efficient use of the large global memory available on the supercomputing cluster enabled by the implementation of shared-memory model provided by the GA tools. Similar phenomenon has also been reported previously in the direct MP2 calculation of morphine molecule using GA tools. 137 The efficient use of global memory has two merits, (a) the cost of memory accesses can be reduced by avoiding the need for virtual memory paging, and (b) the so-called “cache effect” becomes more significant, in particular for a sparser tensor. To understand (b), note that when computation is executed on a large amount of processors, in comparison to a dense tensor, more of significant data in a sparse tensor would be placed in fast memory (i.e. the cache), which makes TWall tend to decrease. If the reduction in TWall offsets increases in communication time and ideal wall time due to the use of additional processors, then the superlinear speedup will occur. Here, the two-electron integral tensors of C90 H182 (1-D) and C96 H122 (2-D) are obviously sparser than that of C91 H84 (3-D). Thus, when Ncore increases from 288 to 384, C90 H182 (1D) and C96 H122 (2-D) exhibit superlinear speedup while the speedup curve of C91 H84 (3-D) is almost linear. When further increasing Ncore to 576, the speedups of C90 H182 (1-D), C96 H122

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

!

"

Figure 3: Scalability test of the low-rank compound decomposition (CD-SVD) of the Coulomb-type two-electron integral tensor in AO space for the studied carbon-hydrogen systems. (a) TWall changes as a function of Ncore for C90 H182 (1-D structure, Nb = 1984, green), C96 H122 (2-D structure, Nb = 1972, blue), and C91 H84 (3-D structure, Nb = 1806, red). The dashed lines are ideal (or linear) speedup curves relative to TWall ’s for Ncore = 288. (b) log10 (TWall) changes as a function of log10 (Nb ). All 1-D, 2-D, and 3-D systems except the largest three are included in (b). k denotes the slope of a linear fitting curve y = k · x + b. For all the tests in (a) and (b), both θCD and θSVD were set to 10−4 .

16


Page 16 of 48

Page 17 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(2-D), and C91 H84 (3-D) with respect to the timing for Ncore = 288 are 2.14, 1.63, and 1.53, respectively. Fig. 3b reveals the actual scaling of our current implementation as a function of system size. As can be seen, for Ncore = 192, the slope k of fitting curve for log10 (TWall ) vs. log10 (Nb ) for all 1-D, 2-D, and 3-D structures (except the largest three) is about 2.6, i.e. TWall numerically scales as ∼ O(Nb2.6 ). This scaling is a bit lower than the previously reported O(Nb3 ) numerical scaling of incomplete CD on alkane chains, 99,138 and the O(m2 Nb2 ) scaling (proportional to O(Nb4 ) for m ∼ O(Nb )) reported in the numerical analysis of other implementations. 136,139 However, it should be noted that the observed scaling is empirical, and can be affected by many factors. As can be seen when Ncore is decrease to 96, k increases to ∼ 3, indicating that too many processors might not be optimal for smaller systems involved in the test. We also compared the wall time of our CD-SVD procedure plus the transformation step (denoted as TCD−SVD+Trans. ) with the wall time of conventional AO-to-MO transformation (denoted as TAO→MO ) for nine selected carbon-hydrogen systems. The results are shown in Fig. S1 in the supplementary information. Generally speaking, CD-SVD plus the transformation step (Eq. (9)) runs faster than the conventional AO-to-MO transformation for all the selected systems. Besides, TCD−SVD+Trans. is more sensitive to θCD than to θSVD which indicates that incomplete CD is the dominant step in our CD-SVD procedure. For a typical pick of θCD (i.e. 10−4 ), TCD−SVD+Trans. is 20 ∼ 40% of TAO→MO . Further tightening the threshold to 10−6 will increases the ratio of TCD−SVD+Trans. to TAO→MO , but more significantly for larger loose-packed 1-D and 2-D systems than for larger dense-packed 3-D systems.

Storage requirement The size of the storage space for the singular vectors obtained from the compound decomposition is related to the system-growth pattern (1-D, 2-D, or 3-D), θCD , θSVD , and Nb . Fig. 4a shows the relationship between the size of storage space and the growth pattern, θCD , 17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and θSVD for selected carbon-hydrogen systems. Given the same θCD and θSVD , and similar Nb , for the three smallest systems, C5 H12 (3-D, Nb = 114) requires larger storage space than C6 H12 (2-D, Nb = 132) and C6 H12 (1-D, Nb = 136). However, this difference becomes vague as the system size increases, e.g. for C16 H34 (1-D, Nb = 356), C16 H26 (2-D, Nb = 340), and C18 H24 (3-D, Nb = 372), the corresponding storage spaces required by the CD-SVD are more or less the same. For each selected structure (in particular structures with smaller size), the ratio of the required storage space in CD-SVD to the conventional storage space behaves like a monotonically decreasing function of the logarithmic value of θCD and/or θSVD . Similar behavior can also be observed when the size of the system is increasing according to a fixed growth pattern (1-D, 2-D, or 3-D). The largest ratio belongs to the smallest 3-D structure C5 H12 , and is only about 50%, which indicates at least half of the storage space can be saved for this structure in comparison to the conventional compact storage of the two-electron integral tensor, i.e. Nt (Nt + 1)/2 with Nt = Nb (Nb + 1)/2. Furthermore, as the system size increases according to the 3-D growth pattern, the ratio will quickly decrease to be < 10% when it reaches C18 H24 . Similar trends can also be found in 1-D and 2-D structures. For a moderate choice of 10−4 for both θCD and θSVD , the relationship between the space size and Nb by using the CD-SVD is more comprehensively exhibited in Fig. 4b,c, where all the studied systems are included. Fig. 4b shows a dramatic decrease in the storage requirement can be observed for the CD-SVD in comparison with the conventional compact storage. In particular, when Nb reaches over 1, 000, the storage space saving of the CD-SVD is over 99% with respect to the conventional compact storage of the full integral tensor. In comparison with the conventional single CD of the compact integral tensor, from which the length of each single Cholesky vector is Nb (Nb + 1)/2, as can be seen from Fig. 4c, CD-SVD starts to outperform when Nb is over ∼ 200 for 1-D structures, ∼ 300 for 2-D structures, and ∼ 350 for 3-D structures, respectively. When Nb is over 1, 000, the CD-SVD can save at least half of the space required by the conventional single CD, and the storage saving will keep increasing as the increasing of Nb (the biggest saving w.r.t. the single CD observed

18


Page 18 of 48

Page 19 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


from the current test almost reaches 80%). To figure out the cause that leads to the dramatic storage savings of CD-SVD, we further analyze the number of Cholesky vectors (denoted as m) and the number of singular vectors per Cholesky vector (denoted as nSVD ) as functions of Nb . As can be seen from Fig. 5a, consistent with other previous studies, 98–101,136,138 m is linearly dependent on Nb . Moreover, we found that the slope of the linear fitting curve between m and Nb shows a linear dependence on log10 (θCD ) (see the inset in Fig. 5a). Fig. 5b reveals the relationship between the average number of singular vectors per Cholesky vector, n ¯ SVD , and Nb for all the studied carbon-hydrogen systems. For each growth pattern, the growth of n ¯ SVD can be numerically well represented by a logarithmic growing curve of Nb . As can be seen from these curves, close-packed structures (e.g. 3-D structures) exhibit larger n ¯ SVD than loose-packed structures (e.g. 1-D structures), i.e. for a given Nb , n ¯ SVD,3−D > n ¯ SVD,2−D > n ¯ SVD,1−D . Despite of the difference, the logarithmic nature of n ¯ SVD as a function of Nb is barely affected by the growth pattern. Therefore, the actual total storage requirement of the CD-SVD, in general, scales approximately as O(m) × O(¯ nSVD ) × Nb ≈ O(Nb2 log10 (Nb )). For the purpose of comparison, the storage requirement of the conventional single CD of an [Nb (Nb + 1)/2] × [Nb (Nb + 1)/2] integral tensor (i.e. storage of m × Nb (Nb + 1)/2 double precision float numbers) is implied by a black dashed line, n ¯ SVD = (Nb + 1)/2, in Fig. 5b. Since each singular vector contains Nb double precision float numbers, the black dashed line then corresponds to the storage of Nb (Nb + 1)/2 double precision float numbers for each Cholesky vector, which is equivalent to the storage requirement using the conventional single CD. Several intersection points between the black dashed line and the colored fitting curves can be observed in the small Nb regime in Fig. 5b. Below these points (i.e. using smaller Nb ’s), the conventional single CD method exhibits slightly smaller storage requirement than the CD-SVD. However, for large-scale applications that usually require tens of thousands of basis functions for the description of the electronic wave function of the system, the O(Nb2 log10 (Nb )) storage required by the CD-SVD apparently outperforms that of the conventional single CD. For example,

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

!

"#############################################$

Figure 4: (a) Ratio of the storage space for singular vectors using the CD-SVD (denoted as SpaceCD−SVD ) to the storage space for the two-electron integral tensor in the compact form in the conventional method (denoted as Spaceconv. ) as a function of the logarithmic value of θCD and θSVD for the selected carbon-hydrogen systems. Totally nine systems were selected from three growth patterns, and 25 tests have been done for each system for θCD and θSVD ranging from 10−6 to 10−2 . (b,c) Storage space saving from the CD-SVD with respect to the conventional storage of the compact integral tensor (b) and to the conventional single CD (c) as a function of Nb for all the studied carbon-hydrogen systems. The storage saving SpaceCD−SVD is defined as 1 − Space . For all the tests in (b,c), both θCD and θSVD were set conv. (or SpaceCD ) −4 to 10 .

20


Page 20 of 48

Page 21 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for 10, 000 basis functions, which is typical for a large system consisting of 500 ∼ 1, 000 atoms, the required storage space by the conventional single CD is 15 ∼ 20 TB, while the storage requirement of the CD-SVD would be just 1 ∼ 5 TB, which may even outperform, for example, the standard ZLIB-type compression of Cholesky vectors that usually provides 2:1 to 5:1 compression (i.e. 50% ∼ 80% storage saving), and is tractable in most of the state-of-the-art supercomputing clusters. Note that the basis set employed in our tests, 6-31+G(d), contains diffuse functions. Usually, the diffuse functions are necessary to be included in the basis set for the purpose of properly describing excited states. However, this will apparently reduce the sparsity of the AO integral tensor in comparison with using non-diffuse basis sets, especially for very large molecules. Generally speaking, the CD of a more sparse AO integral tensor will produce more sparse Cholesky vectors. In particular, the CD of a sparse AO integral tensor with more sparsity in its diagonal will further reduce number of Cholesky vectors. For smaller number of more sparse Cholesky vectors, because n ¯ SVD will scale no higher than O(log10 (Nb )) (as numerically verified in Fig. 5b for 1-D, 2-D, and 3-D systems with Nb ranging from 100 ∼ 2, 000), the storage saving using CD-SVD with respect to the conventional compact storage of the entire AO integral tensor can therefore be expected to be larger.

Accuracy test To test the accuracy, ground and excited state CCSD calculations have been carried out. Fig. 6 shows the computed ground-state CCSD correlation energy ∆E CCSD , the first exLR−CCSD LR−CCSD citation energy ω EOM−CCSD , and static polarizaibilities αxx (the plots for αyy LR−CCSD and αzz are given in Fig. S2 and S3 the supplementary information) of nine selected

carbon-hydrogen systems for different pairs of θCD and θSVD in comparison with the corresponding conventional CCSD results. For the ground state, as can be seen from Fig. 6a, the deviation of ∆E CCSD with respect to conventional results (|∆∆E CCSD |) varies with the change of θCD and θSVD (from 10−6 to 21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

!

"

Figure 5: (a) Linear relationship between the number of Cholesky vectors, m, and Nb for all the studied carbon-hydrogen systems. θCD is ranging from 10−6 to 10−2 . The inset reveals the relationship between the slope of linear fitting curve (y = k · x + b) and log10 (θCD ). The black dashed line shows the dimension of the Coulomb-type two-electron integral tensor, i.e. Nb (Nb + 1)/2. (b) Average number of singular vectors per Cholesky vector (¯ nSVD ) as a function of Nb for all the carbon-hydrogen systems. The fitting function is y = a · log10 (x + c) + b with a, b, and c fitting parameters. For all the tests in (b), both θCD and θSVD were set to 10−4 .

22


Page 22 of 48

Page 23 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


10−2 ) and the system size. For a single structure, no matter in 1-D, 2-D, or 3-D category, |∆∆E CCSD | is reduced with tighter θCD and θSVD . For fixed θCD and θSVD , as the system size grows, along with any pattern, |∆∆E CCSD | slightly grows. Generally, for the system studied here with Nb ranging from 100 to 500, if one wants |∆∆E CCSD | to be < 10−3 a.u. (or < 1 mhartree), the θCD and θSVD need to be < 10−3 . It is worthwhile to mention that a higher agreement may be needed in some other CCSD calculations (e.g. geometry optimization and calculations of energy gradient). In these scenarios, tighter θCD and θSVD (e.g. < 10−6 ) will help reduce |∆∆E CCSD | to be < 1 µHartree for most of the cases in Fig. 6a. To have an idea about how the truncation thresholds used in our CD-SVD procedure affects the energy profile when geometry changes, we performed a series of CCSD(T)/augcc-pVDZ calculations on two typical cases, the bond breaking process of HF molecule and the dissociation of a sandwich-type benzene dimer. The plots are shown in Fig. S4 in the supplementary information. Generally, |∆∆E CCSD(T) | can be well controlled by θCD and θSVD in a wide range of the energy profile. For the first excitation energy, as can be seen from Fig. 6b, the trend of |∆ω EOM−CCSD | with respect to the increase of the system size depends on the tightness of θCD and θSVD used in the compound decomposition. For example, for 1-D growth pattern, if θCD and θSVD are set to 10−3∼−2 (that might be too loose to be ever used in practice), |∆ω EOM−CCSD | will increase with the system size for the studied carbon-hydrogen systems. However, if θCD and θSVD are set to 10−4 , |∆ω EOM−CCSD | will first slightly increase (but still well below 0.1 eV) and then almost stay same for larger system. Further tightening θCD and θSVD (e.g. θCD < 10−4 and θSVD = 10−6 ) might lead to no observable change of |∆ω EOM−CCSD |. Similar complication can also be observed for 2-D and 3-D growth patterns. Fortunately, if both θCD and θSVD are set to be < 10−3 , |∆ω EOM−CCSD | can be well controlled to be well below 0.1 eV for the studied systems. LR−CCSD Finally, for the static polarizaibilities, as shown in Fig. 6c, the change of |∆αxx | LR−CCSD /αxx can be well controlled to be < 0.2%, if θCD and θSVD are < 10−3 (similar situations

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

LR−CCSD LR−CCSD for αyy and αzz , see Fig. S2 and S3). Overall, for the current studied systems,

a moderate choice of 10−4∼−3 can be picked for θCD and θSVD to properly reproduce the conventional CCSD results.

!"""""""""""""""""""""""""""""""""""""""""""

#"""""""""""""""""""""""""""""""""""""""""""""""""$

Figure 6: Deviations of correlation energies (a), first excitation energies (b), and static polarizabilities (c) calculated in the CCSD level for the selected carbon-hydrogen systems along with the change of θCD and θSVD . Calculations with different basis sets have also been performed. The comparison of LR−CCSD LR−CCSD |∆∆E CCSD |, |∆ω EOM−CCSD |, and |∆αxx |/αxx for three smallest structures, each

from one growth pattern, using 6-31+G(d), ZPolX (X = C and H), 141–143 and aug-cc-pVDZ basis sets, is shown in Fig. S5 in the supplementary information. Generally speaking, for larger basis sets, slightly tighter θCD and θSVD (i.e. at least < 10−4 ) might be needed to minimize |∆∆E CCSD | to be well below 1 mHartree. On the other hand, in the excited 24


Page 24 of 48

Page 25 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


state CCSD calculations, θCD and θSVD can be set to be a little loose (e.g. 10−3 to make LR−CCSD LR−CCSD |∆ω EOM−CCSD | < 0.1 eV and |∆αxx |/αxx,conv. < 1%). For fixed θCD and θSVD , LR−CCSD LR−CCSD |∆ω EOM−CCSD | and |∆αxx |/αxx,conv. tend to get smaller with larger basis sets.

A further accuracy test has been done for an isolated C60 molecule. Owing to the delocalized character of the electronic structure of C60 , the proper choice of the basis set and accurate inclusion of the correlation effects are important for obtaining reliable properties, in particular the excited-state properties. A recent intermolecular rigid scan on C60 dimer exhibits a strong dependence of DLPNO-CC and related methods on the basis set effects. 140 The results obtained with Dunning basis sets seemed to be more correct. Also, according to the McAlexander and Crawford, 22 as the dimensionality of the system increases, the corresponding localization errors and computational requirements of using paired natural orbitals (PNO), projected atomic orbitals (PAO), and orbital-specific virtual (OSV) in the computation of dipole polarizabilities and optical rotations, will increase. Here, we employ reduced-size polarized ZPolX basis set (X = C) together with the Dunning aug-cc-pVDZ basis sets for the CCSD calculation of C60 molecule. In particular, we focus on the performance of the CD-SVD on the accurate computation of excitation energy and dipole polarizabilities (static and frequency dependent). The ZPolC and aug-cc-pVDZ basis sets generated 1080 and 1380 basis functions for C60 molecule. As can be seen from Tab 1, the calculated excitation energy associated with 1 T1u ← 1 Ag transition using CD-SVD with θCD and θSVD set to 10−4 is almost same as the one calculated using canonical EOM-CCSD, and the deviation is < 0.002 eV. The excitation energy using aug-cc-pVDZ basis set is within the error bar of experimental results, 144–148 while the one obtained using ZPolC is ∼ 0.1 eV above the upper bound of the error bar. As to the static and dynamic dipole polarizability of the C60 molecule, in the canonical calculations, both ZPolC and aug-cc-pVDZ give reasonable α(ω)’s, which fall into the experimental error bar for ω = 0.0 a.u., and are slightly above the error bar for ω = 0.0428 a.u. (λ = 1064 nm). In particular, when using ZPolC for α(ω) with ω = 0.0428 a.u., the deviation between the computed value and the upper

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

˚3 in comparison with the deviation of ∼ 1.3 bound of the experimental error bar is ∼ 0.7 A ˚ A3 using aug-cc-pVDZ. This can be attributed to the electric field dependence of the ZPolC basis set, and is consistent with our previous study, 149 where we showed that the ZPolX is able to give the polarizabilities of benzene molecule of similar accuracy as those calculated with much larger aug-cc-pVTZ basis set. Furthermore, when compared with the LR-CC2 ˚3 , which indicates the data, 150 the LR-CCSD results apparently outperform by 10 ∼ 12 A importance of T2 -dependent terms in the calculation of α(ω) of a large delocalized system such as C60 . More importantly, for both static and dynamic dipole polarizability, no matter which basis set is used, the difference between the computed value using the CD-SVD and the one from canonical calculations are negligible, and the deviation is < 0.1 a.u. Remarkably, in Table 1, the number of generated Cholesky vectors, m, is smaller than the expected value, i.e. m is only about 2∼3 times Nb in comparison with ∼ 4Nb estimated from Fig. 5a given θCD = 10−4 . This further reduction of m can be attributed to the high symmetry and loose-packing pattern of the C60 molecule which helped construct a more sparse AO integral tensor. For the generated Cholesky vectors, n ¯ SVD still falls into the O(log10 (Nb )) regime as described in Fig. 5b and Tab 1. Overall, due to the small m and n ¯ SVD , the application of CD-SVD in this case requires ∼ 8 GB storage space with ZPolC basis set (∼ 20 GB with aug-cc-pVDZ basis set) for the CCSD calculation of C60 , which are much less than the ∼ 12 GB (∼ 33 GB) storage requirement in the single CD calculation and ∼ 1.2 TB (∼ 3.3 TB) in the canonical calculation.

Conclusion and Outlook In this paper, we propose a compound decomposition strategy to approximately represent the two-electron integral tensor in terms of low-rank vectors with the major aim to address the storage bottleneck in large-scale applications of accurate electronic structure methods. In contrast to the canonical CD of the integral tensor, the compound decomposition combines incomplete CD with a follow-up truncated SVD. A systematic benchmark test on a series 26


Page 26 of 48

Page 27 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Computed transition energy (1 T1u ← 1 Ag ) and dipole polarizabilities (at ω = 0.0 a.u. and 0.0428 a.u.) of C60 molecule. The geometry is adopted from Ref. 151, and is subject to D2h symmetry. All 1s electrons were frozen during calculations. Due to the size of the system, we converged all CCSD iterations to 10−4 , and set both θCD and θSVD to 10−4 . a ω1EOM−CCSD T1u ← 1 Ag CD-SVD Canonical 3.8683 eV 3.8661 eV

α(ω)LR−CCSD ω CD-SVD ZPolC 1080 2456 367 0.0 a.u.b 555.29 a.u. (82.29 ˚ A3 ) 0.0428 a.u.c 564.85 a.u. (83.70 ˚ A3 ) b aug-cc-pVDZ 1380 4338 418 3.5139 eV 3.5130 eV 0.0 a.u. 559.37 a.u. (82.89 ˚ A3 ) c 0.0428 a.u. 569.08 a.u. (84.33 ˚ A3 ) a Experimental values are ranging from 3.04 eV to 3.78 eV. See Ref. 144–148. b Experimental value is 76.5 ± 8 ˚ A3 . See Ref. 152. c Experimental value is 79 ± 4 ˚ A3 . See Ref. 153. Basis set

Nb

m

n ¯ SVD

Canonical 554.71 a.u. (82.20 ˚ A3 ) 564.30 a.u. (83.62 ˚ A3 ) 559.44 a.u. (82.90 ˚ A3 ) 569.15 a.u. (84.34 ˚ A3 )

of 1-D, 2-D, and 3-D carbon-hydrogen systems has been performed. For Nb ranging from ∼ 100 up to ∼ 2, 000, the observed numerical scaling of our implementation shows O(Nb2.5∼3 ) versus O(Nb3∼4 ) cost of performing single CD on the AO integral tensor in most of other implementations. More importantly, this decomposition strategy can significantly reduce the storage requirement of atomic-orbital (AO) based two-electron integral tensor from O(Nb4 ) to O(Nb2 log10 (Nb )) without significant loss of accuracy. The accuracy test has been done by carrying out ground- and excited-state CCSD calculations, in which we found that the decomposition thresholds can be generally set to 10−4∼−3 to give a promising compromise between efficiency and accuracy. From the point of view of future implementations of CC methods, the CD-SVD procedure has far-going consequences stemming from the possibility of re-factoring CC equations (using the contracted singular vectors obtained from the CD-SVD of the two-electron integral tensor) and storing a large portion of data locally at the node or at the local disk drive. The latter factor can help design new algorithms to reduce time-to-solution of CC calculations based on the distribution of all cluster amplitudes across the network and performing contractions with the locally available data. For example, using the contracted singular

27




Page 28 of 48

Page 29 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the Environmental Molecular Sciences Laboratory (EMSL) at the Pacific Northwest National Laboratory (PNNL). EMSL is funded by the Office of Biological and Environmental Research in the U.S. Department of Energy. PNNL is operated for the U.S. Department of Energy by the Battelle Memorial Institute under Contract DE-AC06-76RLO-1830. B.P. acknowledges the Linus Pauling Postdoctoral Fellowship from PNNL. K.K. acknowledges support from the ECP-NWChemEX project.

Supporting Information Available Supporting information includes (1) comparison of the time spent in CD-SVD procedure plus transformation step with the time of conventional AO-to-MO transformation for selected carbon-hydrogen systems, (2) deviations of static polarizabilities (αyy and αzz tensor elements) along with the change of truncation thresholds for selected carbon-hydrogen systems, (3) deviation of CCSD(T) energy as function of geometry change and truncation thresholds in the HF molecule bond breaking and sandwich-type benzene dimer dissociation processes, (4) deviations of CCSD correlation energies, first excitation energies, and static polarizabilities (αxx tensor element) for selected carbon-hydrogen systems with different basis sets along with the change of truncation thresholds, and (5) deviations of the oscillator strength of the first excitation energy calculated in the CCSD level for the selected carbonhydrogen systems along with the change of truncation thresholds. This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Coester, F. Bound States of A Many-Particle System. Nucl. Phys. 1958, 7, 421–424. (2) Coester, F.; K¨ ummel, H. Short-Range Correlations in Nuclear Wave Functions. Nucl. Phys. 1960, 17, 477–485. ˇ ıˇzek, J. On the Correlation Problem in Atomic and Molecular Systems. Calcula(3) C´ 29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tion of Wavefunction Components in Ursell-Type Expansion Using Quantum-Field Theoretical Methods. J. Chem. Phys. 1966, 45, 4256–4266. ˇ ıˇzek, J.; Shavitt, I. Correlation Problems in Atomic and Molecular Sys(4) Paldus, J.; C´ tems. IV. Extended Coupled-Pair Many-Electron Theory and Its Application to the BH3 Molecule. Phys. Rev. A 1972, 5, 50–67. (5) Purvis, G.; Bartlett, R. A Full Coupled-Cluster Singles and Doubles Model: The Inclusion of Disconnected Triples. J. Chem. Phys. 1982, 76, 1910–1918. (6) Scuseria, G. E.; Janssen, C. L.; III, H. F. S. An Efficient Reformulation of the ClosedShell Coupled Cluster Single and Double Excitation (CCSD) Equations. J. Chem. Phys. 1988, 89, 7382–7387. (7) Paldus, J.; Li, X. Adv. Chem. Phys.; John Wiley & Sons, Inc., 2007; pp 1–175. (8) Bartlett, R. J.; Musial, M. Coupled-Cluster Theory in Quantum Chemistry. Rev. Mod. Phys. 2007, 79, 291–352. (9) Lyakh, D. I.; Musia˚ a, M.; Lotrich, V. F.; Bartlett, R. J. Multireference Nature of Chemistry: The Coupled-Cluster View. Chem. Rev. 2012, 112, 182–243. (10) Serrano-Andrés, L.; Merchán, M. Quantum chemistry of the Excited State: 2005 overview. J. Mol. Struct.: TheoCHEM 2005, 729, 99–108. (11) Cársky, P.; Paldus, J.; Pittner, J. Recent Progress in Coupled Cluster Methods: Theory and Applications; Challenges and Advances in Computational Chemistry and Physics; Springer Netherlands, 2010. (12) Krylov, A. I. Equation-of-Motion Coupled-Cluster Methods for Open-Shell and Electronically Excited Species: The Hitchhiker’s Guide to Fock Space. Annu. Rev. Phys. Chem. 2008, 59, 433–462.

30


Page 30 of 48

Page 31 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(13) Jagau, T.-C.; Bravaya, K. B.; Krylov, A. I. Extending Quantum Chemistry of Bound States to Electronic Resonances. Annu. Rev. Phys. Chem. 2017, 68, 525–553. (14) Crawford, T. D.; Tam, M. C.; Abrams, M. L. The Current State of Ab Initio Calculations of Optical Rotation and Electronic Circular Dichroism Spectra. J. Phys. Chem. A 2007, 111, 12057–12068. (15) Crawford, T. D. Ab Initio Calculation of Molecular Chiroptical Properties. Theo. Chem. Acc. 2006, 115, 227–245. (16) Pecul, M.; Ruud, K. In Response Theory and Molecular Properties (A Tribute to Jan Linderberg and Poul Jørgensen); Jensen, H., Ed.; Advances in Quantum Chemistry; Academic Press, 2005; Vol. 50; pp 185 – 212. (17) Kánnár, D.; Tajti, A.; Szalay, P. G. Accuracy of Coupled Cluster Excitation Energies in Diffuse Basis Sets. J. Chem. Theory Comput. 2017, 13, 202–209. (18) Sauer, S. P. A.; Schreiber, M.; Silva-Junior, M. R.; Thiel, W. Benchmarks for Electronically Excited States: A Comparison of Noniterative and Iterative Triples Corrections in Linear Response Coupled Cluster Methods: CCSDR(3) versus CC3. J. Chem. Theory Comput. 2009, 5, 555–564. (19) Kánnár, D.; Szalay, P. G. Benchmarking Coupled Cluster Methods on Valence Singlet Excited States. J. Chem. Theory Comput. 2014, 10, 3757–3765. (20) Watson, T. J.; Lotrich, V. F.; Szalay, P. G.; Perera, A.; Bartlett, R. J. Benchmarking for Perturbative Triple-Excitations in EE-EOM-CC Methods. J. Phys. Chem. A 2013, 117, 2569–2579. (21) Schreiber, M.; Silva-Junior, M. R.; Sauer, S. P. A.; Thiel, W. Benchmarks for Electronically Excited States: CASPT2, CC2, CCSD, and CC3. J. Chem. Phys. 2008, 128, 134110. 31



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(22) McAlexander, H. R.; Crawford, T. D. A Comparison of Three Approaches to the Reduced-Scaling Coupled Cluster Treatment of Non-Resonant Molecular Response Properties. J. Chem. Theory Comput. 2016, 12, 209–222. (23) Geertsen, J.; Rittby, M.; Bartlett, R. J. The Equation-of-Motion Coupled-Cluster Method: Excitation Energies of Be and CO. Chem. Phys. Lett. 1989, 164, 57–62. (24) Comeau, D. C.; Bartlett, R. J. The Equation-of-Motion Coupled-Cluster Method. Applications to Open- and Closed-Shell Reference States. Chem. Phys. Lett. 1993, 207, 414–423. (25) Stanton, J. F.; Bartlett, R. J. The Equation-of-Motion Coupled-Cluster Method. A Systematic Biorthogonal Approach to Molecular Excitation Energies, Transition Probabilities, and Excited State Properties. J. Chem. Phys. 1993, 98, 7029–7039. (26) Monkhorst, H. J. Calculation of Properties with the Coupled-Cluster Method. Int. J. Quantum Chem. 1977, 12, 421–432. (27) Koch, H.; Jørgensen, P. Coupled Cluster Response Functions. J. Chem. Phys. 1990, 93, 3333–3344. (28) Kowalski, K.; Krishnamoorthy, S.; Olson, R. M.; Tipparaju, V.; Aprá, E. Scalable Implementations of Accurate Excited-State Coupled Cluster Theories: Application of High-Level Methods to Porphyrin-Based Systems. 2011 International Conference for High Performance Computing, Networking, Storage and Analysis (SC). 2011; pp 1–10. (29) Kowalski, K.; Krishnamoorthy, S.; Villa, O.; Hammond, J. R.; Govind, N. ActiveSpace Completely-Renormalized Equation-of-Motion Coupled-Cluster Formalism: Excited-State Studies of Green Fluorescent Protein, Free-Base Porphyrin, and Oligoporphyrin Dimer. J. Chem. Phys. 2010, 132, 154103.

32


Page 32 of 48

Page 33 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(30) Valiev, M.; Bylaska, E.; Govind, N.; Kowalski, K.; Straatsma, T.; Dam, H. V.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T. et al. NWChem: A Comprehensive and Scalable Open-Source Solution for Large Scale Molecular Simulations. Comput. Phys. Comm. 2010, 181, 1477–1489. (31) Fan, P.-D.; Valiev, M.; Kowalski, K. Large-Scale Parallel Calculations with Combined Coupled Cluster and Molecular Mechanics Formalism: Excitation Energies of ZincPorphyrin in Aqueous Solution. Chem. Phys. Lett. 2008, 458, 205–209. (32) Pavo˘sević, F.; Peng, C.; Pinski, P.; Riplinger, C.; Neese, F.; Valeev, E. F. Sparse Maps-A Systematic Infrastructure for Reduced Scaling Electronic Structure Methods. V. Linear Scaling Explicitly Correlated Coupled-Cluster Method with Pair Natural Orbitals. J. Chem. Phys. 2017, 146, 174108. (33) Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E. F.; Neese, F. Sparse Maps-A Systematic Infrastructure for reduced-scaling Electronic Structure Methods. II. Linear Scaling domain Based Pair Natural Orbital Coupled Cluster Theory. J. Chem. Phys. 2016, 144, 024109. (34) Flocke, N.; Bartlett, R. J. A Natural Linear Scaling Coupled-Cluster Method. J. Chem. Phys. 2004, 121, 10935–10944. (35) Hughes, T. F.; Flocke, N.; Bartlett, R. J. Natural Linear-Scaled Coupled-Cluster Theory with Local Transferable Triple Excitations: Applications to Peptides. J. Phys. Chem. A 2008, 112, 5994–6003. (36) Ibrahim, K. Z.; Epifanovsky, E.; Williams, S.; Krylov, A. I. Cross-Scale Efficient Tensor Contractions for Coupled Cluster Computations Through Multiple Programming Model Backends. J. Parallel Distrib. Comput. 2017, 106, 92–105. (37) Kaliman, I. A.; Krylov, A. I. New Algorithm for Tensor Contractions on Multi-Core CPUs, GPUs, and Accelerators Enables CCSD and EOM-CCSD Calculations with 33



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Over 1000 Basis Functions on A Single Compute Node. J. Comput. Chem. 2017, 38, 842–853. (38) Kjærgaard, T.; Baudin, P.; Bykov, D.; Kristensen, K.; Jørgensen, P. The DivideExpand-Consolidate Coupled Cluster Scheme. WIREs Comput. Mol. Sci 2017, 9, e1319. (39) Peng, C.; Calvin, J. A.; Pavo˘sević, F.; Zhang, J.; Valeev, E. F. Massively Parallel Implementation of Explicitly Correlated Coupled-Cluster Singles and Doubles Using Tiled Array Framework. J. Phys. Chem. A 2016, 120, 10231–10244. (40) Solomonik, E.; Matthews, D.; Hammond, J.; Stanton, J.; Demmel, J. A Massively Parallel Tensor Contraction Framework for Coupled-Cluster Computations. Tech. Rep. EECS Department, University of California, Berkeley 2014, (41) Baumgartner, G.; Auer, A.; Bernholdt, D. E.; Bibireata, A.; Choppella, V.; Cociorva, D.; Gao, X.; Harrison, R. J.; Hirata, S.; Krishnamoorthy, S. et al. Synthesis of High-Performance Parallel Programs for a Class of Ab Initio Quantum Chemistry Models. Proceedings of the IEEE 2005, 93, 276–292. (42) Deumens, E.; Lotrich, V. F.; Perera, A.; Ponton, M. J.; Sanders, B. A.; Bartlett, R. J. Software Design of ACES III with the Super Instruction Architecture. Wiley Interdisciplinary Reviews: Computational Molecular Science 2011, 1, 895–901. (43) Epifanovsky, E.; Wormit, M.; Ku´s, T.; Landau, A.; Zuev, D.; Khistyaev, K.; Manohar, P.; Kaliman, I.; Dreuw, A.; Krylov, A. I. New Implementation of High-Level Correlated Methods Using A General Block Tensor Library for High-Performance Electronic Structure Calculations. J. Comput. Chem. 2013, 34, 2293–2309. (44) Gao, X.; Krishnamoorthy, S.; Sahoo, S. K.; Lam, C.-C.; Baumgartner, G.; Ramanujam, J.; Sadayappan, P. In Languages and Compilers for Parallel Computing: 18th

34


Page 34 of 48

Page 35 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


International Workshop, LCPC 2005, Hawthorne, NY, USA, October 20-22, 2005, Revised Selected Papers; Ayguadé, E., Baumgartner, G., Ramanujam, J., Sadayappan, P., Eds.; Springer Berlin Heidelberg: Berlin, Heidelberg, 2006; pp 215–229. (45) Hanrath, M.; Engels-Putzka, A. An Efficient Matrix-Matrix Multiplication Based Antisymmetric Tensor Contraction Engine for General Order Coupled Cluster. J. Chem. Phys. 2010, 133, 064108. (46) Hirata, S. Tensor Contraction Engine: Abstraction and Automated Parallel Implementation of Configuration Interaction, Coupled Cluster, and Many Body Perturbation Theories. J. Phys. Chem. A 2003, 107, 9887–9897. (47) Lai, P.-W.; Stock, K.; Rajbhandari, S.; Krishnamoorthy, S.; Sadayappan, P. A Framework for Load Balancing of Tensor Contraction Expressions via Dynamic Task Partitioning. Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis. New York, NY, USA, 2013; pp 1–10. (48) Lotrich, V.; Flocke, N.; Ponton, M.; Sanders, B. A.; Deumens, E.; Bartlett, R. J.; Perera, A. An Infrastructure for Scalable and Portable Parallel Programs for Computational Chemistry. Proceedings of the 23rd International Conference on Supercomputing. New York, NY, USA, 2009; pp 523–524. (49) Rajbhandari, S.; Nikam, A.; Lai, P.-W.; Stock, K.; Krishnamoorthy, S.; Sadayappan, P. Framework for Distributed Contractions of Tensors with Symmetry. The Ohio State University, Tech. Rep. OSU-CISRC 2013, (50) Kállay, M.; Surján, P. R. Higher Excitations in Coupled-Cluster Theory. J. Chem. Phys. 2001, 115, 2945–2954. (51) Kats, D.; Manby, F. R. Sparse Tensor Framework for Implementation of General Local Correlation Methods. J. Chem. Phys. 2013, 138, 144101.

35



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(52) Parkhill, J. A.; Head-Gordon, M. A Sparse Framework for the Derivation and Implementation of Fermion Algebra. Mol. Phys. 2010, 108, 513–522. (53) Lee, J.; Small, D. W.; Epifanovsky, E.; Head-Gordon, M. Coupled-Cluster ValenceBond Singles and Doubles for Strongly Correlated Systems: Block-Tensor Based Implementation and Application to Oligoacenes. J. Chem. Theory Comput. 2017, 13, 602–615. (54) Intel Corporation, Intel®Xeon Phi™Coprocessor System Software Developers Guide. 2014, (55) Peng, B.; Govind, N.; Apr, E.; Klemm, M.; Hammond, J. R.; Kowalski, K. Coupled Cluster Studies of Ionization Potentials and Electron Affinities of Single-Walled Carbon Nanotubes. J. Phys. Chem. A 2017, 121, 1328–1335. (56) Aprà, E.; Kowalski, K. Implementation of High-Order Multireference Coupled-Cluster Methods on Intel Many Integrated Core Architecture. J. Chem. Theory Comput. 2016, 12, 1129–1138. (57) NVIDIA Corporation, NVIDIA®CUDA C Programming Guide. 2013, (58) Jindal, N.; Lotrich, V.; Deumens, E.; Sanders, B. A. Exploiting GPUs with the Super Instruction Architecture. Int. J. Parallel Prog. 2016, 44, 309–324. (59) DePrince, A. E.; Hammond, J. R. Coupled Cluster Theory on Graphics Processing Units I. The Coupled Cluster Doubles Method. J. Chem. Theory Comput. 2011, 7, 1287–1295. (60) Ma, W.; Krishnamoorthy, S.; Villa, O.; Kowalski, K. GPU-Based Implementations of the Noniterative Regularized-CCSD(T) Corrections: Applications to Strongly Correlated Systems. J. Chem. Theory Comput. 2011, 7, 1316–1327.

36


Page 36 of 48

Page 37 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(61) Bhaskaran-Nair, K.; Ma, W.; Krishnamoorthy, S.; Villa, O.; van Dam, H. J. J.; Aprià, E.; Kowalski, K. Noniterative Multireference Coupled Cluster Methods on Heterogeneous CPU-GPU Systems. J. Chem. Theory Comput. 2013, 9, 1949–1957. (62) Greengard, L.; Rokhlin, V. A Fast Algorithm for Particle Simulations. J. Comput. Phys. 1987, 73, 325–348. (63) Strain, M. C.; Scuseria, G. E.; Frisch, M. J. Achieving Linear Scaling for the Electronic Quantum Coulomb Problem. Science 1996, 271, 51–53. (64) White, C. A.; Johnson, B. G.; Gill, P. M.; Head-Gordon, M. Linear Scaling Density functional Calculations via the Continuous Fast Multipole Method. Chem. Phys. Lett. 1996, 253, 268–278. (65) Stratmann, R.; Scuseria, G. E.; Frisch, M. J. Achieving Linear Scaling in ExchangeCorrelation Density Functional Quadratures. Chem. Phys. Lett. 1996, 257, 213–223. (66) Millam, J. M.; Scuseria, G. E. Linear Scaling Conjugate Gradient Density Matrix Search as An Alternative to Diagonalization for First Principles Electronic Structure Calculations. J. Chem. Phys. 1997, 106, 5569–5577. (67) Arnold, D. N. An Interior Penalty Finite Element Method with Discontinuous Elements. SIAM J. Numer. Anal. 1982, 19, 742–760. (68) Lin, L.; Lu, J.; Ying, L.; E, W. Adaptive Local Basis Set for Kohn-Sham Density Functional Theory in A Discontinuous Galerkin Framework I: Total Energy Calculation. J. Comput. Phys. 2012, 231, 2140–2154. (69) Hu, W.; Lin, L.; Yang, C. DGDFT: A Massively Parallel Method for Large Scale Density Functional Theory Calculations. J. Chem. Phys. 2015, 143, 124110. (70) Banerjee, A. S.; Lin, L.; Hu, W.; Yang, C.; Pask, J. E. Chebyshev Polynomial Filtered

37



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Subspace Iteration in the Discontinuous Galerkin Method for Large-Scale Electronic Structure Calculations. J. Chem. Phys. 2016, 145, 154101. (71) Hampel, C.; Werner, H.-J. Local Treatment of Electron Correlation in Coupled Cluster Theory. J. Chem. Phys. 1996, 104, 6286–6297. (72) Sch¨ utz, M.; Werner, H.-J. Low-Order Scaling Local Electron Correlation Methods. IV. Linear Scaling Local Coupled-Cluster (LCCSD). J. Chem. Phys. 2001, 114, 661–681. (73) Sch¨ utz, M. Low-Order Scaling Local Electron Correlation Methods. V. Connected Triples Beyond (T): Linear Scaling Local CCSDT-1b. J. Chem. Phys. 2002, 116, 8772–8785. (74) Kinoshita, T.; Hino, O.; Bartlett, R. J. Singular Value Decomposition Approach for the Approximate Coupled-Cluster Method. J. Chem. Phys. 2003, 119, 7756–7762. (75) Fedorov, D. G.; Kitaura, K. Coupled-Cluster Theory Based upon the Fragment Molecular-Orbital Method. J. Chem. Phys. 2005, 123, 134103. (76) Neese, F.; Wennmohs, F.; Hansen, A. Efficient and Accurate Local Approximations to Coupled-Electron Pair Approaches: An Attempt to Revive the Pair Natural Orbital Method. J. Chem. Phys. 2009, 130, 114108. (77) Li, W.; Piecuch, P.; Gour, J. R.; Li, S. Local Correlation Calculations Using Standard and Renormalized Coupled-Cluster Approaches. J. Chem. Phys. 2009, 131, 114109. (78) Neese, F.; Hansen, A.; Liakos, D. G. Efficient and Accurate Approximations to the Local Coupled Cluster Singles Doubles Method Using a Truncated Pair Natural Orbital Basis. J. Chem. Phys. 2009, 131, 064103. (79) Yang, J.; Chan, G. K.-L.; Manby, F. R.; Schtz, M.; Werner, H.-J. The Orbital-SpecificVirtual Local Coupled Cluster Singles and Doubles Method. J. Chem. Phys. 2012, 136, 144105. 38


Page 38 of 48

Page 39 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(80) Riplinger, C.; Neese, F. An Efficient and Near Linear Scaling Pair Natural Orbital Based Local Coupled Cluster Method. J. Chem. Phys. 2013, 138, 034106. (81) Krause, C.; Werner, H.-J. Comparison of Explicitly Correlated Local Coupled-Cluster Methods with various choices of Virtual Orbitals. Phys. Chem. Chem. Phys. 2012, 14, 7591–7604. (82) Liakos, D. G.; Neese, F. Improved Correlation Energy Extrapolation Schemes Based on Local Pair Natural Orbital Methods. J. Phys. Chem. A 2012, 116, 4801–4816. (83) Li, W.; Piecuch, P. Improved Design of Orbital Domains within the Clusterin-Molecule Local Correlation Framework: Single-Environment Cluster-in-Molecule Ansatz and Its Application to Local Coupled-Cluster Approach with Singles and Doubles. J. Phys. Chem. A 2010, 114, 8644–8657. (84) Li, W.; Piecuch, P. Multilevel Extension of the Cluster-in-Molecule Local Correlation Methodology: Merging Coupled-Cluster and Møller-Plesset Perturbation Theories. J. Phys. Chem. A 2010, 114, 6721–6727. (85) Riplinger, C.; Sandhoefer, B.; Hansen, A.; Neese, F. Natural Triple Excitations in Local Coupled Cluster Calculations with Pair Natural Orbitals. J. Chem. Phys. 2013, 139, 134101. (86) Friedrich, J.; McAlexander, H. R.; Kumar, A.; Crawford, T. D. Incremental Evaluation of Coupled Cluster Dipole Polarizabilities. Phys. Chem. Chem. Phys. 2015, 17, 14284– 14296. (87) Kumar, A.; Crawford, T. D. Frozen Virtual Natural Orbitals for Coupled-Cluster Linear-Response Theory. J. Phys. Chem. A 2017, 121, 708–716. (88) Grasedyck, L.; Kressner, D.; Tobler, C. A literature survey of low-rank tensor Approximation techniques. GAMM-Mitteilungen 2013, 36, 53–78. 39



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(89) Weigend, F. A Fully Direct RI-HF Algorithm: Implementation, Optimised Auxiliary Basis Sets, Demonstration of Accuracy and Efficiency. Phys. Chem. Chem. Phys. 2002, 4, 4285–4291. (90) Vahtras, O.; Almlöf, J.; Feyereisen, M. Integral Approximations for LCAO-SCF Calculations. Chem. Phys. Lett. 1993, 213, 514–518. (91) Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Use of Approximate Integrals in Ab Initio Theory. An Application in MP2 Energy Calculations. Chem. Phys. Lett. 1993, 208, 359–363. (92) Werner, H.-J.; Manby, F. R.; Knowles, P. J. Fast Linear Scaling Second-Order MøllerPlesset Perturbation Theory (MP2) Using Local and Density Fitting Approximations. J. Chem. Phys. 2003, 118, 8149–8160. (93) Sodt, A.; Subotnik, J. E.; Head-Gordon, M. Linear Scaling Density Fitting. J. Chem. Phys. 2006, 125, 194109. (94) Rendell, A. P.; Lee, T. J. Coupled-Cluster Theory Employing Approximate Integrals: An Approach to Avoid the Input/Output and Storage Bottlenecks. J. Chem. Phys. 1994, 101, 400–408. (95) Whitten, J. L. Coulombic Potential Energy Integrals and Approximations. J. Chem. Phys. 1973, 58, 4496–4501. (96) Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. On Some Approximations in Applications of Xα Theory. J. Chem. Phys. 1979, 71, 3396–3402. (97) Boström, J.; Aquilante, F.; Pedersen, T. B.; Lindh, R. Ab Initio Density Fitting: Accuracy Assessment of Auxiliary Basis Sets from Cholesky Decompositions. J. Chem. Theory Comput. 2009, 5, 1545–1553.

40


Page 40 of 48

Page 41 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(98) Røeggen, I.; Wisloff-Nilssen, E. On the Beebe-Linderberg Two-Electron Integral Approximation. Chem. Phys. Lett. 1986, 132, 154–160. (99) Koch, H.; Sánchez de Merás, A.; Pedersen, T. B. Reduced Scaling in Electronic Structure Calculations Using Cholesky Decompositions. J. Chem. Phys. 2003, 118, 9481– 9484. (100) Beebe, N. H. F.; Linderberg, J. Simplifications in the Generation and Transformation of Two-Electron Integrals in Molecular Calculations. Int. J. Quantum Chem. 1977, 12, 683–705. (101) Aquilante, F.; Pedersen, T. B.; Lindh, R. Low-Cost Evaluation of the Exchange Fock Matrix from Cholesky and Density Fitting Representations of the Electron Repulsion Integrals. J. Chem. Phys. 2007, 126, 194106. (102) Røeggen, I.; Johansen, T. Cholesky Decomposition of the Two-Electron Integral Matrix in Electronic Structure Calculations. J. Chem. Phys. 2008, 128, 194107. (103) Neese, F. An Improvement of the Resolution of the Identity Approximation for the Formation of the Coulomb Matrix. Journal of Computational Chemistry 2003, 24, 1740–1747. (104) Aquilante, F.; Boman, L.; Boström, J.; Koch, H.; Lindh, R.; de Merás, A. S.; Pedersen, T. B. In Linear-Scaling Techniques in Computational Chemistry and Physics: Methods and Applications; Zalesny, R., Papadopoulos, M. G., Mezey, P. G., Leszczynski, J., Eds.; Springer Netherlands: Dordrecht, 2011; pp 301–343. (105) Aquilante, F.; Pedersen, T. B.; Lindh, R.; Roos, B. O.; de Meràs, A. S.; Koch, H. Accurate Ab Initio Density Fitting for Multiconfigurational Self-Consistent Field Methods. J. Chem. Phys. 2008, 129, 024113.

41



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(106) Aquilante, F.; Pedersen, T. B. Quartic Scaling Evaluation of Canonical Scaled Opposite Spin Second-Order Møller-Plesset Correlation Energy Using Cholesky Decompositions. Chemical Physics Letters 2007, 449, 354–357. (107) Okiyama, Y.; Nakano, T.; Yamashita, K.; Mochizuki, Y.; Taguchi, N.; Tanaka, S. Acceleration of Fragment Molecular Orbital Calculations with Cholesky Decomposition Approach. Chemical Physics Letters 2010, 490, 84–89. (108) Aquilante, F.; Malmqvist, P.; Pedersen, T. B.; Ghosh, A.; Roos, B. O. Cholesky Decomposition-Based Multiconfiguration Second-Order Perturbation Theory (CDCASPT2): Application to the Spin-State Energetics of CoIII(diiminato)(NPh). J. Chem. Theory Comput. 2008, 4, 694–702. (109) Boström, J.; Delcey, M. G.; Aquilante, F.; Serrano-Andrès, L.; Pedersen, T. B.; Lindh, R. Calibration of Cholesky Auxiliary Basis Sets for Multiconfigurational Perturbation Theory Calculations of Excitation Energies. J. Chem. Theory Comput. 2010, 6, 747–754. (110) Hättig, C.; Weigend, F. CC2 Excitation Energy Calculations on Large Molecules Using the Resolution of the Identity Approximation. J. Chem. Phys. 2000, 113, 5154–5161. (111) Köhn, A.; Hättig, C. Analytic Gradients for Excited States in the Coupled-Cluster Model CC2 Employing the Resolution-of-the-Identity Approximation. J. Chem. Phys. 2003, 119, 5021–5036. (112) Hättig, C. Geometry Optimizations with the Coupled-Cluster Model CC2 Using the Resolution-of-the-Identity Approximation. J. Chem. Phys. 2003, 118, 7751–7761. (113) Weigend, F.; Kattannek, M.; Ahlrichs, R. Approximated Electron Repulsion Integrals: Cholesky Decomposition versus Resolution of the Identity Methods. J. Chem. Phys. 2009, 130, 164106.

42


Page 42 of 48

Page 43 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(114) Chwee, T. S.; Carter, E. A. Density Fitting of Two-Electron Integrals in Local Multireference Single and Double Excitation Configuration Interaction Calculations. Mol. Phys. 2010, 108, 2519–2526. (115) Chwee, T. S.; Carter, E. A. Cholesky Decomposition within Local Multireference Singles and Doubles Configuration Interaction. J. Chem. Phys. 2010, 132, 074104. (116) Werner, H.-J.; Manby, F. R. Explicitly Correlated Second-Order Perturbation Theory Using Density Fitting and Local Approximations. J. Chem. Phys. 2006, 124, 054114. (117) Manby, F. R. Density Fitting in Second-Order linear-r12 Møller-Plesset Perturbation Theory. J. Chem. Phys. 2003, 119, 4607–4613. (118) Katouda, M.; Nakajima, T. MPI/OpenMP Hybrid Parallel Algorithm of Resolution of Identity Second-Order Møller-Plesset Perturbation Calculation for Massively Parallel Multicore Supercomputers. J. Chem. Theory Comput. 2013, 9, 5373–5380. (119) Maschio, L. Local MP2 with Density Fitting for Periodic Systems: A Parallel Implementation. J. Chem. Theory Comput. 2011, 7, 2818–2830. (120) Nakajima, T.; Hirao, K. An Approximate Second-Order Møller-Plesset Perturbation Approach for Large Molecular Calculations. Chem. Phys. Lett. 2006, 427, 225–229. (121) Kendall, R. A.; Fr¨ uchtl, H. A. The Impact of the Resolution of the Identity Approximate Integral Method on Modern Ab Initio Algorithm Development. Theo. Chem. Acc. 1997, 97, 158–163. (122) Bernholdt, D. E.; Harrison, R. J. Large-Scale Correlated Electronic Structure Calculations: the RI-MP2 Method on Parallel Computers. Chem. Phys. Lett. 1996, 250, 477–484.

43



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(123) Maschio, L.; Usvyat, D.; Manby, F. R.; Casassa, S.; Pisani, C.; Sch¨ utz, M. Fast localMP2 Method with Density-Fitting for Crystals. I. Theory and Algorithms. Phys. Rev. B 2007, 76, 075101. (124) Katouda, M.; Nagase, S. Efficient Parallel Algorithm of Second-Order Møller-Plesset Perturbation Theory with Resolution-of-Identity Approximation (RI-MP2). Int. J. Quantum Chem. 2009, 109, 2121–2130. (125) Cremer, D. Møller-Plesset Perturbation Theory: From Small Molecule Methods to Methods for Thousands of Atoms. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 509–530. (126) Pedersen, T. B.; Aquilante, F.; Lindh, R. Density Fitting with Auxiliary Basis Sets from Cholesky Decompositions. Theor. Chem. Acc. 2009, 124, 1–10. (127) Khoromskaia, V.; Khoromskij, B. N.; Schneider, R. Tensor-Structured Factorized Calculation of Two-Electron Integrals in a General Basis. SIAM J. Sci. Comput. 2013, 35, A987–A1010. (128) Kussmann, J.; Ochsenfeld, C. Linear-scaling Method for Calculating nuclear magnetic Resonance chemical shifts Using gauge-including Atomic Orbitals within Hartree-Fock and Density-functional Theory. J. Chem. Phys. 2007, 127, 054103. (129) Maurer, S. A.; Clin, L.; Ochsenfeld, C. Cholesky-decomposed Density MP2 with Density Fitting: Accurate MP2 and Double-hybrid DFT energies for Large Systems. J. Chem. Phys. 2014, 140, 224112. (130) Moritz, G.; Hess, B. A.; Reiher, M. Convergence behavior of the Density-matrix renormalization group Algorithm for Optimized Orbital Orderings. J. Chem. Phys. 2005, 122, 024107.

44


Page 44 of 48

Page 45 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(131) Chan, G. K.-L.; Head-Gordon, M. Highly Correlated Calculations with a polynomial Cost Algorithm: A Study of the Density Matrix renormalization group. J. Chem. Phys. 2002, 116, 4462–4476. (132) Brandhorst, K.; Head-Gordon, M. Fast Sparse Cholesky Decomposition and Inversion Using Nested Dissection Matrix Reordering. J. Chem. Theory Comput. 2011, 7, 351– 368. (133) Peng, B.; Kowalski, K. Low-Rank Factorization of Electron Integral Tensors and Its Application in Electronic Structure Theory. Chem. Phys. Lett. 2017, 672, 47–53. (134) Nieplocha, J.; Harrison, R. J.; Littlefield, R. J. Global Arrays: a Portable “SharedMemory” Programming Model for Distributed Memory Computers. Proceedings of Supercomputing ’94. 1994; pp 340–349. (135) Nieplocha, J.; Harrison, R. J.; Littlefield, R. J. Global Arrays: A Nonuniform Memory Access Programming Model for High-Performance Computers. J. Supercomput. 1996, 10, 169–189. (136) Epifanovsky, E.; Zuev, D.; Feng, X.; Khistyaev, K.; Shao, Y.; Krylov, A. I. General Implementation of the Resolution-of-the-Identity and Cholesky Representations of Electron Repulsion Integrals within Coupled-Cluster and Equation-of-Motion Methods: Theory and Benchmarks. J. Chem. Phys. 2013, 139, 134105. (137) Wong, A. T.; Harrison, R. J.; Rendell, A. P. Parallel Direct Four-Index Transformations. Theoretica chimica acta 1996, 93, 317–331. (138) Krisiloff, D. B.; Krauter, C. M.; Ricci, F. J.; Carter, E. A. Density Fitting and Cholesky Decomposition of the Two-Electron Integrals in Local Multireference Configuration Interaction Theory. J. Chem. Theory Comput. 2015, 11, 5242–5251.

45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(139) Harbrecht, H.; Peters, M.; Schneider, R. On the Low-Rank Approximation by the Pivoted Cholesky Decomposition. Appl. Numer. Math. 2012, 62, 428–440. (140) Sharapa, D. I.; Margraf, J. T.; Hesselmann, A.; Clark, T. Accurate Intermolecular Potential for the C60 Dimer: The Performance of Different Levels of Quantum Theory. J. Chem. Theory Comput. 2017, 13, 274–285. (141) Baranowska, A.; Siedlecka, M.; Sadlej, A. J. Reduced-Size Polarized Basis Sets for Calculations of Molecular Electric Properties. IV. First-Row Transition Metals. Theor. Chem. Acc. 2007, 118, 959–972. (142) Benkova, Z.; Sadlej, A. J.; Oakes, R. E.; Bell, S. E. J. Reduced-Size Polarized Basis Sets for Calculations of Molecular Electric Properties. I. The Basis Set Generation. J. Comput. Chem. 2005, 26, 145–153. (143) Benkova, Z.; Sadlej, A. J.; Oakes, R. E.; Bell, S. E. J. Reduced-Size Polarized Basis Sets for Calculations of Molecular Electric Properties. III. Second-Row Atoms. Theor. Chem. Acc. 2005, 113, 238–247. (144) Krätschmer, W.; Fostiropoulos, K.; Huffman, D. R. The Infrared and Ultraviolet Absorption Spectra of Laboratory-Produced Carbon Dust: Evidence for the Presence of the C60 Molecule. Chem. Phys. Lett. 1990, 170, 167–170. (145) Hare, J.; Kroto, H.; Taylor, R. Preparation and UV / Visible Spectra of Fullerenes C60 and C70 . Chem. Phys. Lett. 1991, 177, 394–398. (146) Lucas, A.; Gensterblum, G.; Pireaux, J. J.; Thiry, P. A.; Caudano, R.; Vigneron, J. P.; Lambin, P.; Krätschmer, W. Elementary Excitations of C60 from the Far Infrared to the Far Vacuum Ultraviolet Studied by High-Resolution Electron-Energy-Loss Spectroscopy. Phys. Rev. B 1992, 45, 13694–13702.

46


Page 46 of 48

Page 47 of 48

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(147) Leach, S.; Vervloet, M.; Desprès, A.; Bréheret, E.; Hare, J. P.; Dennis, T. J.; Kroto, H. W.; Taylor, R.; Walton, D. R. Electronic Spectra and Transitions of the Fullerene C60 . Chem. Phys. 1992, 160, 451–466. (148) Orlandi, G.; Negri, F. Electronic States and Transitions in C60 and C70 Fullerenes. Photochem. Photobiol. Sci. 2002, 1, 289–308. (149) Kowalski, K.; Hammond, J. R.; de Jong, W. A.; Sadlej, A. J. Coupled Cluster Calculations for Static and Dynamic Polarizabilities of C60 . J. Chem. Phys. 2008, 129, 226101. (150) Pedersen, T. B.; de Mers, A. M. J. S.; Koch, H. Polarizability and Optical Rotation Calculated from the Approximate Coupled Cluster Singles and Doubles CC2 Linear Response Theory Using Cholesky Decompositions. J. Chem. Phys. 2004, 120, 8887– 8897. (151) Ruud, K.; Jonsson, D.; Taylor, P. R. The Dispersion of the Polarizability of C60 : A Confirmation of Recent Experimental Results through Theoretical Calculations. J. Chem. Phys. 2001, 114, 4331–4332. (152) Antoine, R.; Dugourd, P.; Rayane, D.; Benichou, E.; Broyer, M.; Chandezon, F.; Guet, C. Direct Measurement of the Electric Polarizability of Isolated C60 Molecules. J. Chem. Phys. 1999, 110, 9771–9772. (153) Ballard, A.; Bonin, K.; Louderback, J. Absolute Measurement of the Optical Polarizability of C60 . J. Chem. Phys. 2000, 113, 5732–5735.

47



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(pq|rs) ≈

m " k=1

Lkpq Lkrs , with Lkpq =

n "

k(n)

k(n)

k(n)

U pδ Λδδ U qδ

δ=1

!"#$%&'()'*+,$'-'." /0,.1'2(3."'&$%4

om O(Nb4 ) t 50,.1'2(67(8'9"#$:

d O(Nb3 ) n 50,.1'2(670!;7(8'9"#$:

o O(Nb2 log10 (Nb ))

48


Page 48 of 48

Highly Efficient and Scalable Compound Decomposition of Two

Recommend Documents