State-Averaged Pair Natural Orbitals for Excited States: A Route

Subscriber access provided by University of Sunderland

Quantum Electronic Structure

State-Averaged Pair Natural Orbitals for Excited States: A Route Towards Efficient Equation of Motion Coupled-Cluster Chong Peng, Marjory C Clement, and Edward F. Valeev J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00171 • Publication Date (Web): 25 Sep 2018 Downloaded from http://pubs.acs.org on October 1, 2018

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 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 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.

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 36 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

State-Averaged Pair Natural Orbitals for Excited States: A Route Towards Efficient Equation of Motion Coupled-Cluster Chong Peng, Marjory C. Clement, and Edward F. Valeev∗ Department of Chemistry, Virginia Tech, Blacksburg, VA 24061, USA E-mail: [email protected] Abstract A reduced-complexity variant of equation-of-motion coupled-cluster singles and doubles (EOM-CCSD) method is formulated in terms of state-averaged excited state pair natural orbitals (PNO) designed to describe manifolds of excited states. Stateaveraged excited state PNOs for the target manifold are determined by averaging CIS(D) pair densities over the model manifold. The performance of PNO-EOM-CCSD approach has been tested with the help of a distributed-memory parallel canonical EOM-CCSD implementation within the Massively Parallel Quantum Chemistry program that allows treatment of systems with 50+ atoms using realistic basis sets with 1000+ functions. The use of state-averaged PNOs offers several potential advantages relative to the recently proposed state-specific PNOs: our approach is robust with respect to root flipping and state degeneracies, it is more economical when computing large manifolds of states, and it simplifies evaluation of transition-specific observables such as dipole moments. With the PNO truncation threshold of 10−7 , the errors in excitation energies are on average below 0.02 eV for the first six singlet states of 28 organic molecules included in the standard test set of Thiel and co-workers (J. Chem. Phys. 2008, 128, 134110) with 50-70 state-averaged PNOs per pair.

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

1

Introduction

Accurate description of electronic spectra of medium (< 100 atoms) and large (> 100 atoms) molecular systems has always been a challenge for quantum chemistry. The time-dependent density functional theory (TDDFT) is the most popular method for the analysis of excited states due to its computational efficiency, capable of treatment of systems with hundreds and thousands of atoms. Although TDDFT provides medium accuracy for one-electron excitations, the accuracy of TDDFT can be limited for certain types of excited states (e.g. Rydberg or charge-transfer) 1 and in general its accuracy depends strongly on the density functional. 2 In contrast to TDDFT, multiconfiguration/multireference (MR) wave function models, such as MR perturbation theory methods (e.g. complete active space perturbation theory (CASPT2) 3 and n-electron valence states perturbation theory (NEVPT2) 4 ) and MR configuration interaction 5,6 can recover both static and dynamic electron correlation, can treat multiple electronic states on equal footing, and attain high accuracy, albeit for rather small systems. 7 Among the challenges of the MR approaches is the need to select the active space, and the exponential growth of complexity with the size of active space. Although the latter can be avoided for certain types of systems by numerical approximations such as density matrix renormalization group 8,9 and other tensor network approaches, the MR methods are generally more difficult to use for nonspecialists. Accurate treatment of dynamical electron correlation in the context of MR methodologies is an ongoing direction of research. In this work we focus on the treatment of excited states by the coupled-cluster method. The highly robust coupled-cluster hierarchy provides unparalleled accuracy for the ground states by systematically including two-, three- and higher-body correlation effects from a single determinant reference. The CC ansatz can be extended to excited states through the use of the linear-response (LR) theory, 10 the symmetry-adapted cluster configuration interaction (SAC-CI) method, 11,12 or the equation of motion coupled-cluster (EOM-CC) method. 13,14 However, the high-order scaling of the coupled-cluster methods limits its application to small molecules. Even with truncation to singles and doubles excitations, the excited state CCSD 2


Page 2 of 36

Page 3 of 36 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


methods still have polynomial scaling with large factor O (N 6 ), and are constrained to systems containing only 20-30 atoms without access to campus-level or national computing resources. Recently, the development of reduced scaling variants of the coupled-cluster methods has been reinvigorated by Neese’s introduction 15 of Pair Natural Orbitals (PNOs) in the context of local correlation formalisms of CC initiated by Pulay 16 and pursued by Werner 17,18 and others. 19,20 PNOs were originally proposed in the 1960s under the name Pseudo-natural Orbitals. 21–24 Truncation of PNOs significantly reduces the number of unoccupied orbitals while only introducing small errors in correlation energies in post-Hartree-Fock calculations. However, the demanding computational cost of the pair-specific integral transformation to the PNO space, which scales O (N 7 ) if there is no truncation of the PNO space, prevents the development of PNO-based electronic structure theories. In 2009, Neese et al. “revived” the PNOs by using them in the framework of local correlation theories (hence, the LPNO moniker) for the CEPA 15 and CCSD 25 methods, making use of density fitting approximations to accelerate the integral transformation process. It made large-scale coupled-cluster possible for systems with up to 100 atoms using a small workstation. The LPNO approach was improved by imposing block sparsity into cluster operator amplitudes via domains of projected atomic orbitals (PAOs). 26,27 Compared to the pioneering linear-scaling PAO-based coupled-cluster methods developed by Werner and co-workers, 17,28 and other linear-scaling coupled-cluster approaches, 29 the use of PAOs and PNOs allows to greatly reduce the prefactor, while with linear-scaling local density fitting 30,31 the PAO-based PNO methods maintain formal linear scaling. Crucial to such development is the introduction of explicit correlation to reduce the basis set error, with pioneering ideas coming from Werner, Tew, and others 32–39 These developments are pursued in parallel by several groups, with polynomial scaling explicitly correlated PNO-CCSD(T) code demonstrated by Hättig, Tew, and co-workers, 40–42 linear scaling PNO-CCSD(T)-F12 demonstrated by Neese, Valeev, and co-workers, 36,37 and a scalable implementation of a linear-scaling PNO-CCSD(T)-F12 demonstrated by Werner

3



and co-workers. 38,39,43,44 Open-shell PNO coupled-cluster has also been demonstrated. 45 The key ideas of modern reduced-scaling coupled-cluster methods apply not only to the ground states but also to the excited states. Two competing visions have been proposed for how to formulate reduced-scaling excited state coupled-cluster methods. Hättig and Helmich have explored excited-state coupled-cluster methods by introducing PNO-EOMCC2 with state-specific PNOs; 46 the method had formal O (N 4 ) cost complexity but the dominant contributions were O (N 3 ). This effort was followed by the PNO based CIS(D) 47 and ADC(2)-x; 48 as we were revising this manuscript, Hättig and Frank also presented PNO-based EOM-CCSD for excitation energies. 49 The common theme of these methods is the use of state-specific PNOs to compress the cluster operator (computed in the ground state CC equation) and the excited state wave operators for each state, with excited states computed one at a time. Thus the total number of PNOs grows linearly with the number of excited states; this makes the cost of the integral transformation to grow linearly with the number of states. Alternative direction was explored by Dutta et al who presented a PNO-based coupled-cluster method for excited states utilizing the similarity-transformed EOM (STEOM) CCSD framework. 50,51 In their approach the use of PNOs is limited to the ground state only, with DLPNO-CCSD amplitudes subsequently transformed to the canonical basis and used to evaluate the bt-PNO-STEOM-CCSD energies of manifolds of states, at a O (N 6 ) complexity but reducible O (N 5 ) with additional improvements. Unlike the methods of Hättig et al. which utilize sparsity in ground and excited state wave functions, bt-PNO-STEOM-CCSD does not use sparsity for the excited state wave functions, thus limiting the size of the system that can be considered. We should also note that the use of local correlation ideas (namely, PAO domains) without the PNO-style compression, have been explored in the context of coupled-cluster methods for excitation energies, such as local EOM-CC2 52–54 and local EOM-CCSD. 18,19 Lastly, the PNO framework has been using with multireference perturbation theories, e.g. DLPNO-NEVPT2 55 and LPNO-CASPT2, 56 to treat excited states; in both cases state-specific PNO formalism was used.

4


Page 4 of 36

Page 5 of 36 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


In this work, we present a PNO-based approach suitable for robust treatment of manifolds of excited states with the EOM-CCSD methods. The key idea is to use state-averaged PNOs similar to those used in the ground-state PNO coupled-cluster methods through stateaveraged guess pair densities averaged over the target excited state manifold. To quickly explore the performance of our approach we simulated it using a parallel canonical (i.e. naive complexity) EOM-CCSD implementation. An efficient parallel canonical EOM-CCSD was implemented in the open-source Massively Parallel Quantum Chemistry (MPQC) package, 57 based on the ground-state CCSD implementation described previously. 58 The new implementation exhibits good strong-scaling parallel performance and allows the calculation of excitation energy for systems with more than 50 atoms and more than 1000 basis functions; this is crucial to the exploration of the state-averaged PNO ansatz for systems of realistic size. In Section 2, the theory and implementation of state-averaged PNOs are discussed. Section 3 describes the computational details as well as the computing resources used. Section 4 reports on the performance of the state-averaged PNOs for the prediction of EOM-CCSD excitation energies.

2

Methods

The coupled-cluster ground-state wave function, ˆ

Ψ(0) ≡ eT |0i ,

(1)

where the |0i stands for the zeroth order reference wave function (usually a Hartree-Fock determinant), is determined by projection of the Schrödinger equation against excited de-

5



Page 6 of 36

terminants

¯ |0i , 0 = hai | H

(2)

¯ 0 = hab ij | H |0i .

(3)

¯ ≡ e−Tˆ HeTˆ the usual similarity-transformed Hamiltonian. Within the equation with H of motion coupled-cluster method 13,14 kth excited-state wave function is obtained in a CI fashion, by the action of a linear excitation operator acting on the ground-state CC wave function: ˆ (k) eTˆ |0i , Ψ(k) ≡ R

(4)

ˆ (k) and the corresponding energies E(k) are obtained by diagonalizing the similarity-transformed R Hamiltonian: ¯R ˆ (k) |0i = E(k) R ˆ (k) |0i . H

(5)

In practice the ground and excited states are represented in terms of single and double excitations only: CCSD Tˆ ≡ Tˆ1 + Tˆ2 ,

(6)

Tˆ1 ≡ Tai Eia ,

(7)

1 ij ab Eij . Tˆ2 ≡ Tab 2

(8)

ˆ (k) CCSD ˆ 1(k) + R ˆ 2(k) , R ≡ δk0 + R

(9)

i ˆ 1(k) ≡ Ba(k) R Eia ,

(10)

ˆ 2(k) ≡ 1 B ij E ab , R 2 ab(k) ij

(11)

6


Page 7 of 36 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 Eia and Eijab are the spin-free one- and two-particle replacers: Eia ≡

X

a†aτ aiτ ,

(12)

τ =α,β

Eijab ≡

X

a†aτ a†bσ ajσ aiτ .

(13)

τ,σ=α,β

The Einstein summation convention utilized throughout this paper. The storage and operation costs of the CCSD and EOM-CCSD methods are O (N 4 ) and O (N 6 ), respectively. ij ij Two-body amplitude tensors, Tab and Bab , are efficiently rank-compressed by transform-

ing each ij block into the ij-specific subspace. For the ground-state amplitudes the optimal pair-specific subspaces are robustly approximated by the truncated singular subspace of the corresponding ground-state pair densities, Dij (0) , computed from guess amplitudes: Dij (0) = where Tij

ab

2 ˜ ij† + Tij† T ˜ ij ), (Tij T 1 + δij

(14)

ij ˜ ij ≡ 2Tij − Tij† . (Semi)canonical MP1 amplitudes, ≡ Tab and T

(T (1) )ij ab ≡

Gij ab , j i fi + fj − faa − fbb

(15)

are typically used as the guess 25 (In Eq. (15) f are the matrix elements of the Fock operator, and G stands for the two-electron integral).

Gqs pr ≡ hpr|qsi

(16)

PNOs are the basis for the singular subspace of Dij (0) ; they are obtained by solving the eigensystem: ij ij ij Dij (0) U(0) = U(0) n(0) .

(17)

where nij (0) are the PNO occupation numbers. PNOs with occupation numbers less than user-

7



Page 8 of 36

provided threshold TCutPNO are omitted, hence the number of PNOs per pair is independent of the system size (i.e. O (1)). One-body amplitudes, Tai and Bai , are compressed similarly to the two-body counterparts by transforming into the basis of orbital-specific virtuals (OSVs). 59 OSVs are traditionally defined to be identical to the PNOs of diagonal pairs but truncated according to a different threshold, TCutOSV . As pointed out by Hättig and Helmich, 47 the optimal singular subspaces for the ground and excited state amplitudes differ; as a result, the PNOs and OSVs must be constructed separately for the ground and excited states. Hättig and Helmich proposed the use of state-specific PNOs, where the PNOs for each state are constructed using CIS(D) doubles amplitudes with respect to that state: 46

ij Bab(k)

=

ij Kab(k)

ω(k) + fii + fjj − faa − fbb

,

ij j ij ij i ic l l Kab(k) = Bc(k) Gcj ab + Bc(k) Gab − Ba(k) Glb − Bb(k) Gal ,

(18)

(19)

where Bia(k) and Bij ab(k) are the CIS singles amplitudes and CIS(D) doubles amplitudes for excited state k, and ω(k) is the CIS excitation energy. The state-specific PNOs for excited states can be obtained from the state-specific pair density using the CIS(D) doubles amplitudes similar to the approach used in ground state:

Dij (k) =

2 ij† ˜ ij ˜ ij† (Bij (k) B(k) + B(k) B(k) ), 1 + δij

(20)

Such definition of excited state PNOs yields good accuracy in the context of PNO-EOM-CC2 method. 46 However, there are several factors that prompted us to look beyond the statespecific PNOs. First, and foremost, the cost of PNO construction and integral transformation grow linearly with the number of states. This is particularly notable since the cost of PNObased methods is often dominated by the cost of the integral transformation, even when domain approximations are employed. 30,31 Second, state-specific PNOs can make it more

8


Page 9 of 36 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


difficult to deal with degenerate state manifolds (which ideally need to be expressed in the basis invariant to mixings of orbitals from degenerate manifolds). Lastly, the use of statespecific PNOs increases the complexity of formalism and implementation. Thus we decided to investigate PNO-EOM-CCSD that uses one set of PNOs for all excited states, in particular, we propose to use the state-averaged PNOs. The state-averaged PNOs are defined as the eigenvectors of averaged pair densities over an N -state manifold: N 1 X ij D(k) . D = N k ij

(21)

(State-averaged OSVs will be defined in this work the PNOs of the diagonal pairs, in complete analogy with the construction of the ground-state OSVs). Note that motivated by the aforementioned potential shortcomings of state-averaged PNOs Helmich and Hättig briefly discussed the possibility of using “merged PNOs” as a single basis for describing multiple excited states in Ref. 46 but deemed the idea not likely to be practical and did not pursue the idea further. We briefly tested the “merged” PNO approach (in which truncated PNOs obtained for each state separately are concatenated and orthogonalized) against the state-averaged approach and found that the former leads to a much more rapid grows of the average number of “merged” PNOs with the number of states, e.g. for water dimer in aug-cc-pVDZ basis and TCutPNO = 10−8 the average number of merged PNOs per pair was 68, compared to 44 state-averaged PNOs. This can be easily understood: merging multiple sets of state-specific PNOs does not distinguish between similar PNOs, whereas in computing state-averaged PNOs similar state-specific PNOs will be mixed and represented by a smaller set of state-averaged PNOs. We expect this to be a typical behavior for the merged PNOs and thus expect the merged PNO approach to be significantly more costly to the state-averaged approach that we proposed here. Although the work is underway in our group to develop a production implementation of reduced-scaling CC, here our goal is more modest: we aim to evaluate the proposed state-

9



Page 10 of 36

averaged PNO formulation in the context of EOM-CCSD. Hence we initially implemented a simulation for PNO-EOM-CCSD based on a newly-developed massively parallel canonical (i.e., O (N 6 )) EOM-CCSD program in the MPQC code. Note that simulation has been used previously for initial evaluation of locally-correlated PAO-based EOM-CCSD by Russ and Crawford 19,20 and by Korona and Werner. 18 Werner et al. and Crawford et al. also used simulation to compare PAO-, OSV-, and PNO-based formulations of CCSD. 60,61 It should also be noted that Hättig and Helmich have demonstrated production-level O (N 4 ) PNOEOM-CC2 methods, 46 and a cubic scaling PNO-EOM-CCSD implementation at the time of writing this manuscript. 49 The ground-state PNO-CCSD simulation was implemented with modification to the Jacobi update in the following manner: 1. After the CCSD amplitude residuals R1 and R2 are computed, the unoccupied space of R1 is transformed into a semi-canonical OSV basis and R2 is transformed into a semi-canonical PNO basis: ¯ i = U a Ri , R ai ai a

(22)

ij b a ¯ ij R aij bij = Uaij Rab Ubij ,

(23)

where ai and aij are unoccupied orbitals in the truncated OSV and PNO basis, respectively. 2. The residuals are updated through a Jacobi update in the OSV and PNO space:

¯i = ∆ ai

¯ ij ∆ aij bij =

¯i R ai , fii − f¯aaii

¯ ij R aij bij , a b j i fi + fj − f¯aijij − f¯bijij

(24)

(25)

3. The updated residuals are extrapolated with DIIS and back-transformed into the

10


Page 11 of 36 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


canonical basis: ¯ ia , ∆ia = Uaai ∆ i

(26) b

ij aij ¯ ij ∆ij ab = Ua ∆aij bij Ub .

(27)

4. The new CCSD amplitudes are formed as an update to the current amplitudes:

Tin+1 = Tin + ∆i ,

(28)

ij ij Tij n+1 = Tn + ∆ ,

(29)

where n stands for the number of current iteration. 5. The CCSD residuals are recomputed using the new amplitudes, and this process is repeated from step 1 until convergences is reached. Similarly, the state-averaged PNO simulation in EOM-CCSD can be done with modification in the Davidson solver: 1. The unoccupied space of residuals produced by the Davidson algorithm is transformed into the OSV and PNO bases:

a i ¯i R ai (k) = Uai Ra(k) ,

(30)

ij a b ¯ ij R aij bij (k) = Uaij Rab(k) Ubij .

(31)

2. A preconditioner is applied to the residuals in the OSV and PNO spaces:

¯i B ai (k) =

¯ ij B aij bij (k)

=

¯i R ai (k)

, ω(k) + fii − f¯aaii ¯ ij R aij bij (k)

a b ω(k) + fii + fjj − f¯aijij − f¯bijij

11


(32)

,

(33)


Page 12 of 36

where ω(k) is the eigenvalue of state k. 3. The unoccupied space of updated trial vectors is projected back into the canonical space: i ¯i , Ba(k) = Uaai B ai (k)

(34)

bij ij ¯ ij Bab(k) = Uaaij B aij bij (k) Ub .

(35)

4. The new trial vectors are added to the next iteration of the Davidson algorithm to update the subspace, and the process is continued from step 1 until it reached convergence.

3

Computational Details

The canonical EOM-CCSD code was implemented and tested in the developmental version of the MPQC program. 57 All computations in Section 4.1 were performed on a commodity cluster at Virginia Tech, each node of which has two Intel Xeon E5-2670 CPUs (332 GFLOPS) and 64 GB of RAM. MPQC was compiled using GCC 5.3.0 with Intel MPI 5.0 and the serial Intel MKL version 11.2.3. All computations launched 1 MPI process per node with 16 threads per MPI process, with the orbital block size set to 20. The state-averaged PNO-EOM-CCSD simulation code was also implemented in MPQC. For simplicity we truncate PNOs and OSVs using single parameter, i.e. TCutOSV =TCutPNO /10 for both ground and excited states. Neither domain nor weak pair approximations were utilized. The occupied MOs were localized in all calculations via the Foster-Boys algorithm. 62,63 The density-fitting (resolution-of-identify) approximation 64,65 and frozen core approximation were used for all the calculations performed in this work. We have used the cc-pVTZ 66 and aug-cc-pVD/TZ 67 atomic orbital basis sets in our calculations, with the corresponding auxiliary basis sets cc-pVTZ-RI and aug-cc-pVD/ TZ-RI 68 for density-fitting. In Section 4.1, the geometries of the methylated uracil dimer with water were obtained from Ref. 69 and 12


Page 13 of 36 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 structure of 11-cis-retinal protonated Schiff base was obtained from Ref. 50. In Section 4.2, the structures of benzonitrile and phenolate form of the anionic chromophore of the photoactive yellow protein (PYPb) were obtained from Ref. 7. In Section 4.3, a total number of 10 excited-states were computed for the benchmark dataset of 28 organic molecules by Thiel. 7 In Section 4.4, geometry of acetamide were obtained from Ref. 7 and the geometries of 4-(N, N-dimethylamino) benzonitrile (DMABN), N-phenyl pyrrole (PP) and dipeptide were obtained from Ref. 70.

4

Results & Discussion

4.1

Distributed-Memory Canonical EOM-CCSD Implementation

The canonical closed-shell EOM-CCSD program in MPQC was implemented on top of the TiledArray tensor framework following the formalism of Bartlett and Stanton. 14 The implementation details generally follow the ground-state explicitly correlated CCSD implementation reported previously. 58 All amplitudes and intermediates are distributed in memory, and contractions are evaluated using the communication-optimal implementation of the distributed-memory scalable universal matrix multiplication algorithm (SUMMA) implemented in TiledArray. 71,72 Several variants of the EOM-CCSD have been implemented, with and without the use of density fitting, and with integrals stored explicitly in memory or evaluated lazily. The latter allows to reduce the storage requirements of EOM-CCSD to allow calculations on systems with over 1000 basis functions. The new canonical EOM-CCSD code can attain high efficiency and good parallel scalability as illustrated in Fig. 1 and Fig. 2 for realistic computations (with aug-cc-pVDZ and cc-pVTZ basis sets) on excited-states of the methylated uracil dimer with water and 11-cis-retinal protonated Schiff base, respectively.

The data in Fig. 1 corresponds to a ∼ 5 speedup from 16 to 128 nodes with the

aug-cc-pVDZ basis and a ∼ 3 speedup from 32 to 128 nodes with the cc-pVTZ basis. The data in Fig. 2 corresponds to a ∼ 5.5 speedup from 16 to 128 nodes with the aug-cc-pVDZ 13



EOM-CCSD Time per Iteration(min)

64

aDZ TZ

32 16 8 4

16

32

Number of Nodes

64

128

Figure 1: Parallel performance of EOM-CCSD on 4 states of the methylated uracil dimer with water (39 atoms) with the aug-cc-pVDZ (645) and cc-pVTZ (882) basis sets

64

Time per Iteration(min)

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 14 of 36

aDZ TZ

32

16

8

16

32

Number of Nodes

64

128

Figure 2: Parallel performance of EOM-CCSD on 4 states of the 11-cis-retinal protonated Schiff base (51 atoms) with the aug-cc-pVDZ (753) and cc-pVTZ (1050) basis sets

14


Page 15 of 36 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


basis and a ∼ 1.6 speedup from 64 to 128 nodes with the cc-pVTZ basis. The demonstrated strong scaling is not as robust as that of the ground-state CCSD program; 58 nevertheless, the performance of our code is already sufficient to be able to treat multiple states of a system with 50-100 atoms and 1000-1500 basis functions.

4.2

Accuracy of State-Averaged PNOs

To quantify the performance of state-averaged PNOs we computed errors in excitation energies relative to the canonical EOM-CCSD values introduced by the truncation of PNOs (and the corresponding truncation of OSVs). Table 1 and Table 2 list the PNO truncation errors for benzonitrile in the cc-pVTZ and aug-cc-pVTZ basis for a fixed value of the TCutPNO parameter, as a function of the number of computed states. Table 1: Truncation errors in excitation energy (eV) of benzonitrile (cc-pVTZ,V=283 a ) with respect to total number of states at TCutPNO =10−8

nStates ESnPNO S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 MAE b MAX c a b c

1 63

2 79

4 91

6 95

8 101

10 104

20 133

0.0072 0.0016 0.0014 0.0003 -0.0001 -0.0001 -0.0011 0.0043 0.0023 0.0008 0.0005 0.0004 -0.0005 0.0235 0.0039 0.0019 0.0019 0.0005 0.0263 0.0022 0.0022 0.0022 0.0009 0.0161 0.0173 0.0171 0.0040 0.0069 0.0049 0.0042 0.0004 0.0876 0.0235 0.0050 0.1834 0.0557 0.0012 0.0201 -0.0001 0.0058 0.0044 0.0072 0.0030 0.0134 0.0050 0.0372 0.0131 0.0018 0.0072 0.0043 0.0263 0.0161 0.1834 0.0557 0.0050

Total number of unoccupied orbitals Mean absolute error Max absolute error It is difficult to estimate how the increase in the number of PNOs in the state-averaged

approach will affect its cost relative to the state-specific approach. The cost of evaluating 15



Page 16 of 36

Table 2: Truncation errors in excitation energy (eV) of benzonitrile (aug-cc-pVTZ,V=456 a ) with respect to total number of states at TCutPNO =10−8

nStates ESnPNO S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 MAE b MAX c a b c

1 74

2 93

4 108

6 111

8 113

0.0159 0.0036 0.0027 0.0019 0.0084 0.0033 0.0022 0.0074 0.0062 0.0150 0.0045 0.0072 0.0250

0.0025 0.0024 0.0032 0.0036 0.0067 0.0242 0.0237 0.0024

10 118

20 118

0.0017 0.0025 0.0022 0.0033 0.0020 0.0001 0.0019 0.0009 0.0050 0.0055 0.0051 0.0059 0.0257 0.0298 0.0016 0.0003 0.0024 -0.0005 0.0050 0.0000 0.0159 0.0060 0.0071 0.0078 0.0086 0.0053 0.0049 0.0159 0.0084 0.0150 0.0250 0.0242 0.0257 0.0298

Total number of unoccupied orbitals Mean absolute error Max absolute error

integrals with k PNO indices scales with the average number of PNOs n as O nk . If we assume that the 4-external PNO integrals are the dominant expense, then as long as the increase in the number of PNOs due to state-averaging over m states is less than m1/4 the cost of the single PNO transformation in the state averaged approach will be smaller than the cost of m such transformations in the state-specific approach. The benzonitrile data suggests that increasing the number of computed states from 1 to 20 increases the number of PNOs only by a factor of ∼ 2.1 in cc-pVTZ basis and a factor of ∼ 1.6 in aug-cc-pVTZ 1

1

basis, which is smaller than m 4 (20 4 ) in this case. Clearly, the total number of state-averaged PNOs grows with the number of states sufficiently slowly relative to the linear growth of the total number of state-specific PNOs used by Hättig and co-workers. Hence, the use of stateaveraged PNOs should offer some savings in the costs of the integral transformation. Note that this comparison is actually somewhat pessimistic since the errors in excitation energies also decrease as the total number of states increases because of the concomitant increase 16


Page 17 of 36 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


in the number of PNOs. For example, the absolute truncation error for state 1 decreased from 0.0072 eV to 0.0011 eV from 1 state to 20 states, which is about 6.5 times smaller, as demonstrated in Table 1. Overall our approach to PNO construction is rather accurate: the mean absolute errors are below 0.02 eV for all cases except nStates=8 in Table 1, which has a mean absolute error of 0.037 eV. These errors are small relative to the average accuracy of the EOM-CCSD model even for states with single excitation character. The slow growth of the number of state-averaged PNOs with respect to the number of states should not be entirely surprising; since the low-energy states in molecules to zeroth order have many occupied orbitals in common; correlation effects will be largely similar among the states. State-averaging is able to combine mostly similar PNO vectors for pairs that are not involved in excitations. It is however worth mentioning that it is possible to design a situation where the average number of state-averaged PNOs will grow more rapidly (linear with number of states at worst), e.g. if there is small or no overlap between low-lying excitations and if most electron pairs are affected by excitations. Note that the mean absolute errors do not smoothly decrease as nStates increases. This is correlated with sporadic increases in the maximum absolute errors as the number of states is increased, such as in the case of states 3 and 4 when nStates is at 4 and states 7 and 8 when nStates is at 8 in Table 1. However, these errors are significantly reduced when nStates is increased to 6 and 10, respectively. This indicates that the highest excited states in the computed manifold sometimes have more significant errors with state-averaged PNOs, which can be observed from the nStates = 4,8 data in Table 1. The reason for this behavior is that the composition of N lowest EOM-CCSD states may not be similar to that of CIS, either due to pure root flipping or, more generally, due to nonperturbative effects of dynamical correlation on the excited state character and ordering. Analysis of the excited states in this example suggests that CIS states 3, 4, 5, and 6 become EOM-CCSD states 5, 6, 3, and 4 in the cc-pVTZ basis. Therefore accurate description of EOM-CCSD states 3 and 4 will require including pair densities from CIS(D) states 5 and 6. To reach faster convergence of

17



Page 18 of 36

truncation error for the targeted states, more states should be computed than the targeted states due to the root flipping issue, which increased the prefactor of the computation time. This is not a serious issue since in excited state computations to increase the probability that N lowest-energy target states have been reproduced is to compute M > N states. Hence, a slightly larger error in a few of the highest excited states would not be an issue since typically the number of computed states is always greater than the number of target states. Lastly, note that the TCutPNO threshold is kept constant in Table 1 and Table 2. Therefore the averaged errors decrease as the number of states increases, at the cost of increasing the average number of state-averaged PNOs per pair. Clearly, if we wanted to keep the average error per state constant we could loosen the TCutPNO threshold as the number of states is increased. This would further alleviate the modest increase of the total number of PNOs with the number of states. Dependence of the error on TCutPNO will be examined next. Tables 3, 4 and 5 illustrate the correlation between the TCutPNO parameter and the errors in the excitation energies of the 4 lowest singlet excited states of the phenolate form of the anionic chromophore of the photoactive yellow protein (PYPb). Since TCutPNO affects ˆ the EOM-CCSD excitation energies through both ground-state (Tˆ) and excited-state (R) operators, we examined its effects separately on the ground-state cluster operators only (Table 3), excited-state operators (Table 4) and both (Table 5). As expected (see Table 3) Table 3: Truncation error in excitation energy (eV) of PYPb (aug-cc-pVDZ, V=296 a ) by only truncating the ground-state PNOs

TCutPNO GSnPNO S1 S2 S3 S4 a

10−6 13

10−7 25

10−8 45

10−9 74

10−10 110

-0.0872 -0.0892 -0.0855 -0.0875

-0.0289 -0.0301 -0.0287 -0.0296

-0.0092 -0.0098 -0.0093 -0.0098

-0.0030 -0.0033 -0.0031 -0.0033

-0.0011 -0.0012 -0.0012 -0.0013

Total number of unoccupied orbitals

truncating the ground-state PNOs only lowers the excitation energies (the errors in excited 18


Page 19 of 36 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 4: Truncation error in excitation energy (eV) of PYPb (aug-cc-pVDZ, V=296 a ) by only truncating the excited-state PNOs

TCutPNO ESnPNO S1 S2 S3 S4 a

10−6 6

10−7 15

10−8 32

10−9 58

10−10 94

0.1323 0.1215 0.1354 0.1296

0.0564 0.0525 0.0635 0.0611

0.0211 0.0197 0.0256 0.0261

0.0053 0.0047 0.0078 0.0092

0.0015 0.0011 0.0024 0.0035


Table 5: Truncation error in excitation energy (eV) of PYPb (aug-cc-pVDZ, V=296 a ) by truncating both the ground and excited-state PNOs

TCutPNO GSnPNO ESnPNO S1 S2 S3 S4 a

10−6 13 6

10−7 25 15

10−8 45 32

0.0457 0.0332 0.0500 0.0422

0.0276 0.0226 0.0348 0.0315

0.0119 0.0098 0.0163 0.0163

10−9 74 58

10−10 110 94

0.0022 0.0003 0.0013 -0.0002 0.0046 0.0012 0.0058 0.0022


states are all negative) since the energy of the ground state becomes higher due to a decrease in the amount of correlation energy that is recovered. Similarly, truncating the excited-state PNOs raises the excitation energy since the calculated energy of the excited states is now higher as a result of recovering less of the correlation energy (Table 4). When both the ground and excited state PNOs are truncated, the opposite signs of the two sources of error partially cancel each other out, leading to smaller errors in the excitation energy, as can be seen in Table 5. However, this error cancellation may lead to occasional non-monotonic convergence.

19



4.3

Error Analysis

To further test the performance of state-averaged PNOs, we used the PNO-EOM-CCSD method to compute the excitation energies of the lowest six singlet excited states of 28 organic molecules in the benchmark dataset from Thiel et al. 7 Variation of statistical measures of the errors with the TCutPNO parameter are presented in Fig. 3 and Fig. 4 for cc-pVTZ and aug-cc-pVTZ basis sets, respectively. The corresponding average numbers of ground-state and excited-state PNOs are shown in Fig. 5 and Fig. 6, respectively. 1.0000 MAE MAX 0.1000

Error in eV

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 20 of 36

0.0100

0.0010

0.0001

6

7

8 -log(TCutPNO )

9

10

Figure 3: Mean absolute (MAE) and maximum (MAX) PNO truncation errors (in eV) of PNO-EOM-CCSD/cc-pVTZ excitation energies for the 28-molecule benchmark set. The average excitation energy errors become smaller than 0.1 eV already with TCutPNO =10−6 , further reduce to below 0.02 eV with TCutPNO =10−7 , and further decrease monotonically with TCutPNO for both basis sets. The maximum errors also decrease monotonically, but require tighter truncation, TCutPNO ={10−7 , 10−8 } for the {cc-pVTZ,aug-cc-pVTZ} basis set, to reduce below 0.1 eV. Hence, a threshold of 10−7 is suitable for general applications while a threshold of 10−8 should be sufficient for high-accuracy applications. With these thresholds, the average numbers of ground/excited state PNOs are ∼50/70 and ∼100/120 for cc-pVTZ and aug-cc-pVTZ basis, respectively, which is a significant decrease over the average number 20


Page 21 of 36

1.0000 MAE MAX

Error in eV

0.1000

0.0100

0.0010

0.0001

6

7

8 -log(TCutPNO )

9

10

Figure 4: Mean absolute (MAE) and maximum (MAX) PNO truncation errors (in eV) of PNO-EOM-CCSD/aug-cc-pVTZ excitation energies for the 28-molecule benchmark set.

350 GSnPNO ESnPNO Virtual

300 Average Number of Orbitals

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


250 200 150 100 50 0

6

7

8 -log(TCutPNO )

9

10

Figure 5: Convergence of average number of PNOs per pair per molecule in ground state (GSnPNO) and excited states (ESnPNO) of PNO-EOM-CCSD/cc-pVTZ for the 28-molecule benchmark set.

21



350 300 Average Number of Orbitals

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 22 of 36

250 200 150 100 GSnPNO ESnPNO Virtual

50 0 6

7

8 -log(TCutPNO )

9

10

Figure 6: Convergence of average number of PNOs per pair per molecule in ground state (GSnPNO) and excited states (ESnPNO) of PNO-EOM-CCSD/aug-cc-pVTZ for the 28molecule benchmark set. of total virtual orbitals. For all values of TCutPNO the number of excited state PNOs is only 20%-30% higher than the corresponding number of ground state PNOs (GSnPNO).

4.4

Rydberg and Charge Transfer States

We also studied the accuracy of state-averaged PNO-EOM-CCSD for excited states with Rydberg and charge-transfer character. We selected one example (Table 6) for Rydberg character: acetamide . We also selected four examples (Table 7) of states with charge-transfer character: ethylene-tetrafluoroethylene (C2 H4 -C2 F4 ) 73 , 4-(N,N-dimethylamino) benzonitrile (DMABN), N-phenyl pyrrole (PP) and dipeptide. 70,74 The ethylene-tetrafluoroethylene model was also used to test other PNO-based excited state methods, by Hättig and Helmich 46 and by Dutta et al. 50,51 The truncation errors for the two Rydberg states (S1 and S4 ) were found to be somewhat larger than the errors for the valence states (S2 and S3 ) with TCutPNO =10−6 but they are comparable with TCutPNO =10−7 or tighter. Overall, no significant differences in the performance of excited state PNOs is observed for Rydberg and non-Rydberg states. 22


Page 23 of 36 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 6: Truncation errors ( eV) of the aug-cc-pVTZ a PNO-EOM-CCSD excitation energies of the four lowest singlet states of acetamide . States S1 and S4 have strong Rydberg character. Excited-state PNOs were averaged over lowest 10 states.

TCutPNO ESnPNO S1 S2 S3 S4 a

10−6 32

10−7 69

10−8 125

10−9 191

10−10 245

0.0699 0.0052 -0.0017 -0.0009 -0.0002 0.0038 -0.0058 -0.0032 -0.0012 -0.0005 0.0185 -0.0001 -0.0028 -0.0011 -0.0003 0.0464 0.0012 -0.0023 -0.0004 0.0000

Total number of unoccupied orbitals = 283.

Table 7: Truncation errors (eV) of PNO-EOM-CCSD excitation energies of local (L) and charge-transfer (CT) states on various molecules. Excited-state PNOs were averaged over lowest 10 states.

Molecule

State Type

C2 H4 -C2 F4

DMABN

PP

b

b

Dipeptide

a b

b

a

10−6

10−7

10−8

10−9

10−10

0.1202 0.0166 0.0304 0.1730

0.0327 0.0015 0.0054 0.0622

0.0090 0.0001 0.0003 0.0215

0.0020 0.0000 0.0000 0.0059

0.0006 0.0000 0.0000 0.0014

S1 S3 S4 S7

CT L L CT

S1 S2

L 0.1349 0.0228 CT 0.1847 0.0340

0.0013 -0.0006 -0.0004 0.0013 -0.0014 -0.0010

S1 S2 S3 S4

L L CT CT

0.1082 0.1446 0.2159 0.2124

0.0145 0.0243 0.0494 0.0444

0.0005 -0.0003 -0.0001 0.0014 -0.0007 -0.0004 0.0064 0.0001 -0.0004 0.0097 0.0014 0.0001

S1 S2 S4 S8

L L CT CT

0.0940 0.1052 0.1436 0.1391

0.0126 -0.0002 -0.0009 0.0155 -0.0002 -0.0013 0.0229 0.0023 0.0008 0.0532 0.0126 0.0039

-0.0004 -0.0006 -0.0002 -0.0007

Computed with aug-cc-pVTZ. Computed with cc-pVTZ. The truncation errors for the charge-transfer states were found to be larger than the

valence states at loose truncation thresholds. TCutPNO =10−7 was required to reduce the errors to below 0.1 eV, and TCutPNO =10−8 is sufficient to reduce the errors to 0.05 eV for 23



Page 24 of 36

Table 8: Average number of PNOs per pair in ground state (GSnPNO) and excited states (ESnPNO) for molecules in Table 7.

Molecule

Type 10−6

10−7

10−8

10−9

10−10

0

C2 H4 -C2 F4

GSnPNO ESnPNO

14 8

30 21

54 46

81 85

112 428 136 428

DMABN

GSnPNO ESnPNO

17 17

36 42

66 85

112 148

171 431 223 431

PP

GSnPNO ESnPNO

20 19

41 45

74 87

122 144

181 418 212 418

Dipeptide

GSnPNO ESnPNO

14 16

29 39

53 77

88 130

132 375 193 375

charge-transfer excitations. This is in agreement with the findings of Hättig and Helmich, 47 who pointed out that more PNOs are needed to get the same accuracy for charge-transfer excitations relative to their valence counterparts. They attributed this to the use of the semicanonical CIS(D) amplitudes in constructing the excited-state PNOs (Eq. (18)); with localized occupied orbitals the off-diagonal matrix elements of the Fock operator are substantial and cannot be neglected. They have introduced the Laplace-transformed CIS(D), which reduces the truncation error for both local and charge-transfer states. 46 Nevertheless, the performance of semicanonical PNOs is still acceptable, hence we did not pursue the use of Laplace transform technique to for computing guess amplitudes in the local basis; this will be rectified in the future. The average number of PNOs per pair for molecules in Table 7 are listed in Table 8. The average number of PNOs per pair are around 40 and 80 at truncation threshold TCutPNO =10−7 and TCutPNO =10−8 for states in Table 7 , which contain charge-transfer character.

24


Page 25 of 36 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


5

Conclusions and Perspective

We proposed a state-averaged PNO ansatz for efficient and simple treatment of manifolds of excited states in the context of reduced-scaling excited-state many-body methods. We evaluated the performance of the state-averaged PNO ansatz in the context of PNO-EOMCCSD method for prediction of excitation energies. The PNO-EOM-CCSD implementation is based on a new distributed-memory parallel canonical implementation of EOM-CCSD in the MPQC program. The state-averaged PNO-EOM-CCSD approach has been tested on the first six excited states of 28 organic molecules, yielding an average truncation error below 0.020 eV at TCutPNO =10−7 for both the cc-pVTZ and aug-cc-pVTZ basis sets. With this truncation threshold, the number of state-averaged PNOs is reduced by more than 70% for cc-pVTZ and 80% for aug-cc-pVTZ. For charge-transfer excitations tested, a truncation threshold of 10−7 reduces the average truncation error to below 0.05 eV. Overall, the stateaveraged PNOs provide excellent accuracy for low lying valence and Rydberg states, but tighter threshold is required to achieve the same accuracy for charge-transfer states. As is the case with state averaging approaches familiar from multireference theories, the definition of state-averaged PNOs can depend strongly on the choice of averaging states; e.g. by averaging over 19 valence states and 1 Rydberg state the PNOs important for describing the Rydberg state may be underweighted. However, the numerical data in Section 4 is sufficiently promising to warrant further investigation of state-averaged PNOs. Furthermore, it should be possible to augment the state-averaged PNO definition to ensure that PNOs important for each individual state and not well represented by the state-average PNO set are included. The performance of state-averaged PNOs is sufficiently encouraging to warrant the development of a production-level PNO-EOM-CCSD code based on the state-averaged PNO definitions introduced here. The work along these lines is underway, and will be reported elsewhere.

25



6

Acknowledgments

This work was partly supported by the Exascale Computing Project (ECP), Project Number: 17-SC-20-SC, a collaborative effort of two DOE organizations: the Office of Science and the National Nuclear Security Administration; by the Virginia Tech Institute for Critical Technology and Applied Science (ICTAS) (M. C.); and by the U.S. National Science Foundation (awards 1362655 and 1550456). The authors thank the Advanced Research Computing (ARC) at Virginia Tech for making computational resources available.

26


Page 26 of 36

Page 27 of 36 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


References (1) Dreuw, A.; Head-Gordon, M. Single-reference ab initio methods for the calculation of excited states of large molecules. Chem. Rev. 2005, 105, 4009–4037. (2) Jacquemin, D.; Wathelet, V.; Perpète, E. A.; Adamo, C. Extensive TD-DFT benchmark: Singlet-excited states of organic molecules. J. Chem. Theory Comput. 2009, 5, 2420–2435. (3) Andersson, K.; Malmqvist, P.; Roos, B. O. Second–order perturbation theory with a complete active space self–consistent field reference function. J. Chem. Phys. 1992, 96, 1218–1226. (4) Angeli, C.; Cimiraglia, R.; Evangelisti, S.; Leininger, T.; Malrieu, J. P. Introduction of n-electron valence states for multireference perturbation theory. J. Chem. Phys. 2001, 114, 10252. (5) Werner, H. J.; Knowles, P. J. An efficient internally contracted multiconfigurationreference configuration interaction method. J. Chem. Phys. 1988, 89, 5803–5814. (6) Szalay, P. G.; M¨ uller, T.; Gidofalvi, G.; Lischka, H.; Shepard, R. Multiconfiguration selfconsistent field and multireference configuration interaction methods and applications. Chem. Rev. 2012, 112, 108–181. (7) Schreiber, M.; Silva-Junior, M. R.; Sauer, S. P.; Thiel, W. Benchmarks for electronically excited states: CASPT2, CC2, CCSD, and CC3. J. Chem. Phys. 2008, 128, 134110. (8) Chan, G. K.-L.; Sharma, S. The density matrix renormalization group in quantum chemistry. Annu. Rev. Phys. Chem. 2011, 62, 465–481. (9) Schollwöck, U. The density-matrix renormalization group. Rev. Mod. Phys. 2005, 77, 259–315.

27



(10) Monkhorst, H. J. Calculation of properties with the coupled–cluster method. Int. J. Quantum Chem. 1977, 12, 421–432. (11) Nakatsuji, H. Cluster expansion of the wavefunction, valence and Rydberg excitations, ionizations, and inner-valence ionizations of CO2 and N2 O studied by the SAC and SAC CI theories. Chem. Phys. 1983, 75, 425–441. (12) Nakatsuji, H. Cluster expansion of the wavefunction. Electron correlations in ground and excited states by SAC (symmetry-adapted-cluster) and SAC CI theories. Chem. Phys. Lett. 1979, 67, 329–333. (13) Sekino, H.; Bartlett, R. J. A linear response, coupled–cluster theory for excitation energy. Int. J. Quantum Chem. 1984, 26, 255–265. (14) 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. (15) 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 . (16) Pulay, P. Localizability of dynamic electron correlation. Chem. Phys. Lett. 1983, 100, 151–154. (17) Hampel, C.; Werner, H. Local treatment of electron correlation in coupled cluster theory. J. Chem. Phys. 1996, 104, 6286–6297. (18) Korona, T.; Werner, H. J. Local treatment of electron excitations in the EOM-CCSD method. J. Chem. Phys. 2003, 118, 3006–3019. (19) Daniel Crawford, T.; King, R. A. Locally correlated equation-of-motion coupled cluster theory for the excited states of large molecules. Chem. Phys. Lett. 2002, 366, 611–622. 28


Page 28 of 36

Page 29 of 36 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


(20) Russ, N. J.; Crawford, T. D. Local correlation in coupled cluster calculations of molecular response properties. Chem. Phys. Lett. 2004, 400, 104–111. (21) Edmiston, C.; Krauss, M. Configuration-Interaction Calculation of H3 and H2. J Chem Phys 1965, 42, 1119–1120. (22) Edmiston, C.; Krauss, M. Pseudonatural orbitals as a basis for the superposition of configurations. I. He+ 2 . J. Chem. Phys. 1966, 45, 1833–1839. (23) Meyer, W. Ionization energies of water from PNO-CI calculations. Int. J. Quantum Chem. 1971, 5, 341–348. (24) Ahlrichs, R.; Lischka, H.; Staemmler, V.; Kutzelnigg, W. PNO–CI (pair natural orbital configuration interaction) and CEPA–PNO (coupled electron pair approximation with pair natural orbitals) calculations of molecular systems. I. Outline of the method for closed–shell states. J Chem Phys 1975, 62, 1225–1234. (25) 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 . (26) Riplinger, C.; Neese, F. An efficient and near linear scaling pair natural orbital based local coupled cluster method. J. Chem. Phys. 2013, 138, 034106. (27) 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. (28) Sch¨ utz, M.; Werner, H.-J. Local perturbative triples correction (T) with linear cost scaling. Chem. Phys. Lett. 2000, 318, 370–378. (29) Scuseria, G. E.; Ayala, P. Y. Linear scaling coupled cluster and perturbation theories in the atomic orbital basis. J Chem Phys 1999, 111, 8330–8343. 29



(30) Pinski, P.; Riplinger, C.; Valeev, E. F.; Neese, F. Sparse maps - A systematic infrastructure for reduced-scaling electronic structure methods. I. An efficient and simple linear scaling local MP2 method that uses an intermediate basis of pair natural orbitals. J. Chem. Phys. 2015, 143, 034108. (31) 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. (32) Werner, H. J. Eliminating the domain error in local explicitly correlated second-order Møller–Plesset perturbation theory. J Chem Phys 2008, 129, 101103. (33) Tew, D. P.; Helmich, B.; Hättig, C. Local explicitly correlated second-order M?ller– Plesset perturbation theory with pair natural orbitals. J Chem Phys 2011, 135, 074107. (34) Tew, D. P.; Hättig, C. Pair natural orbitals in explicitly correlated second-order MøllerPlesset theory. Int. J. Quantum Chem. 2013, 113, 224–229. (35) Pavoˇsević, F.; Neese, F.; Valeev, E. F. Geminal-spanning orbitals make explicitly correlated reduced-scaling coupled-cluster methods robust, yet simple. J. Chem. Phys. 2014, 141, 054106. (36) Pavoˇsević, F.; Pinski, P.; Riplinger, C.; Neese, F.; Valeev, E. F. SparseMaps—A systematic infrastructure for reduced-scaling electronic structure methods. IV. Linearscaling second-order explicitly correlated energy with pair natural orbitals. J. Chem. Phys. 2016, 144, 144109. (37) Pavoˇsević, F.; Peng, C.; Pinski, P.; Riplinger, C.; Neese, F.; Valeev, E. F. SparseMaps 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. 30


Page 30 of 36

Page 31 of 36 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


(38) Ma, Q.; Werner, H. J. Scalable electron correlation methods. 2. Parallel PNO-LMP2F12 with near linear scaling in the molecular size. J. Chem. Theory Comput. 2015, 11, 5291–5304. (39) Ma, Q.; Werner, H.-J. Scalable electron correlation methods. 5. Parallel perturbative triples correction for explicitly correlated local coupled cluster with pair natural orbitals. J. Chem. Theory Comput. 2017, 14, acs.jctc.7b01141. (40) Schmitz, G.; Helmich, B.; Hättig, C. A scaling PNO–MP2 method using a hybrid OSV– PNO approach with an iterative direct generation of OSVs. Molecular Physics 2013, 111, 2463–2476. (41) Schmitz, G.; Hättig, C.; Tew, D. P. Explicitly correlated PNO-MP2 and PNO-CCSD and their application to the S66 set and large molecular systems. Phys Chem Chem Phys 2014, 16, 22167–22178. (42) Schmitz, G.; Hättig, C. Perturbative triples correction for local pair natural orbital based explicitly correlated CCSD(F12) using Laplace transformation techniques. J. Chem. Phys. 2016, 145, 234107. (43) Schwilk, M.; Ma, Q.; Köppl, C.; Werner, H.-J. Scalable electron correlation methods. 3. Efficient and accurate parallel local coupled cluster with pair natural orbitals (PNOLCCSD). J. Chem. Theory Comput. 2017, 13, 3650–3675. (44) Ma, Q.; Werner, H.-J. Scalable electron correlation methods. 5. Parallel perturbative triples correction for explicitly correlated local coupled cluster with pair natural orbitals. J. Chem. Theory Comput. 2018, 14, 198–215. (45) Saitow, M.; Becker, U.; Riplinger, C.; Valeev, E. F.; Neese, F. A new near-linear scaling, efficient and accurate, open-shell domain-based local pair natural orbital coupled cluster singles and doubles theory. J Chem Phys 2017, 146, 164105.

31



(46) Helmich, B.; Hättig, C. A pair natural orbital implementation of the coupled cluster model CC2 for excitation energies. J. Chem. Phys. 2013, 139, 084114. (47) Helmich, B.; Hättig, C. Local pair natural orbitals for excited states. J. Chem. Phys. 2011, 135, 214106. (48) Helmich, B.; Hättig, C. A pair natural orbital based implementation of ADC(2)-x: Perspectives and challenges for response methods for singly and doubly excited states in large molecules. Comput. Theor. Chem. 2014, 1040-1041, 35–44. (49) Frank, M. S.; Hättig, C. A pair natural orbital based implementation of CCSD excitation energies within the framework of linear response theory. J. Chem. Phys. 2018, 148, 134102. (50) Dutta, A. K.; Neese, F.; Izsák, R. Towards a pair natural orbital coupled cluster method for excited states. J. Chem. Phys. 2016, 145, 034102. (51) Dutta, A. K.; Nooijen, M.; Neese, F.; Izsák, R. Exploring the accuracy of a low scaling similarity transformed equation of motion method for vertical excitation energies. J. Chem. Theory Comput. 2018, 14, 72–91. (52) Kats, D.; Korona, T.; Sch¨ utz, M. Local CC2 electronic excitation energies for large molecules with density fitting. J. Chem. Phys. 2006, 125, 104106. (53) Freundorfer, K.; Kats, D.; Korona, T.; Sch¨ utz, M. Local CC2 response method for triplet states based on Laplace transform: Excitation energies and first-order properties. J. Chem. Phys. 2010, 133, 244110. (54) Kats, D.; Sch¨ utz, M. A multistate local coupled cluster CC2 response method based on the Laplace transform. J. Chem. Phys. 2009, 131, 124117. (55) Guo, Y.; Sivalingam, K.; Valeev, E. F.; Neese, F. SparseMaps - A systematic infrastructure for reduced-scaling electronic structure methods. III. Linear-scaling multireference 32


Page 32 of 36

Page 33 of 36 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


domain-based pair natural orbital N-electron valence perturbation theory. J. Chem. Phys. 2016, 144, 094111. (56) Menezes, F.; Kats, D.; Werner, H. J. Local complete active space second-order perturbation theory using pair natural orbitals (PNO-CASPT2). J. Chem. Phys. 2016, 145, 124115. (57) Peng, C.; Lewis, C.; Wang, X.; Clement, M.; Pavosevic, F.; Zhang, J.; Rishi, V.; Teke, N.; Pierce, K.; Calvin, J.; Kenny, J.; Seidl, E.; Janssen, C.; Valeev, E. F. The Massively Parallel Quantum Chemistry Program (MPQC). 2018; available under an open source license from www.mpqc.org. (58) 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 TiledArray framework. J. Phys. Chem. A 2016, 120, 10231–10244. (59) Yang, J.; Kurashige, Y.; Manby, F. R.; Chan, G. K. Tensor factorizations of local second-order Møller-Plesset theory. J. Chem. Phys. 2011, 134, 044123. (60) 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. (61) McAlexander, H. R.; Crawford, T. D. A comparison of three approaches to the reducedscaling coupled cluster treatment of non-Resonant molecular response properties. J. Chem. Theory Comput. 2016, 12, 209–222. (62) Foster, J. M.; Boys, S. F. Canonical configurational interaction procedure. Rev. Mod. Phys. 1960, 32, 300–302. (63) Boys, S. F. Construction of some molecular orbitals to be approximately invariant for changes from one molecule to another. Rev. Mod. Phys. 1960, 32, 296–299. 33



(64) 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. (65) Vahtras, O.; Almlöf, J.; Feyereisen, M. W. Integral approximations for LCAO-SCF calculations. Chem. Phys. Lett. 1993, 213, 514–518. (66) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (67) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of the first–row atoms revisited. Systematic basis sets and wave functions. J. Chem. Phys. 1992, 96, 6796–6806. (68) Weigend, F.; Köhn, A.; Hättig, C. Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations. J. Chem. Phys. 2002, 116, 3175–3183. (69) 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. (70) Peach, M. J. G.; Benfield, P.; Helgaker, T.; Tozer, D. J. Excitation energies in density functional theory: An evaluation and a diagnostic test. J. Chem. Phys. 2008, 128, 044118. (71) Calvin, J. A.; Valeev, E. F. Task-based algorithm for matrix multiplication: A step towards block-sparse tensor computing. 2015, 9. (72) Calvin, J. A.; Lewis, C. A.; Valeev, E. F. Scalable task-based algorithm for multiplication of block-rank-sparse matrices. Proc. 5th Work. Irregul. Appl. Archit. Algorithms IA3 ’15 2015, 1–8. 34


Page 34 of 36

Page 35 of 36 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


(73) Dreuw, A.; Weisman, J. L.; Head-Gordon, M. Long-range charge-transfer excited states in time-dependent density functional theory require non-local exchange. J. Chem. Phys. 2003, 119, 2943–2946. (74) Rohrdanz, M. A.; Martins, K. M.; Herbert, J. M. A long-range-corrected density functional that performs well for both ground-state properties and time-dependent density functional theory excitation energies, including charge-transfer excited states. J. Chem. Phys. 2009, 130, 054112.

35



Page 36 of 36

ES3 ES3 Density Average

PNO Truncation

ES2 ES2 Density

State-Averaged Density

State-Averaged PNOs

ES1 ES1 Density

PNO Truncation

GS GS Density


GS PNOs

State-Averaged Pair Natural Orbitals for Excited States: A Route

Recommend Documents