Efficient Implementation of the Second-Order Quasidegenerate

Jul 18, 2019 - The high cost of common multireference second-order perturbation theory ... second-order quasidegenerate perturbation theory (MC-QDPT2)...
14 downloads 0 Views 847KB Size
Subscriber access provided by UNIV OF SOUTHERN INDIANA

Quantum Electronic Structure

An Efficient Implementation of the Second-Order Quasidegenerate Perturbation Theory with Density-Fitting and Cholesky Decomposition Approximations: Is It Possible to Use HartreeFock Orbitals for a Multiconfigurational Perturbation Theory? U#ur Bozkaya J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.9b00378 • Publication Date (Web): 18 Jul 2019 Downloaded from pubs.acs.org on July 18, 2019

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

An Efficient Implementation of the Second-Order Quasidegenerate Perturbation Theory with Density-Fitting and Cholesky Decomposition Approximations: Is It Possible to Use Hartree–Fock Orbitals for a Multiconfigurational Perturbation Theory? Uğur Bozkaya∗ Department of Chemistry, Hacettepe University, Ankara 06800, Turkey E-mail: [email protected]



To whom correspondence should be addressed

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

Abstract The high cost of common multireference second-order perturbation theory (MRPT2) methods compared with the single-reference variant (MP2) arises from the expensive complete active space self-consistent field (CASSCF) orbital optimization step. Furthermore, the use of conventional four-index electron repulsion integrals (ERIs) prevents their application to larger molecular systems due to expensive I/O procedures. To address these bottlenecks of the multiconfigurational second-order quasidegenerate perturbation theory (MC-QDPT2), an efficient implementation of QDPT2 with the density-fitting (DF) and Cholesky decomposition (CD) approximations, denoted by DF-QDPT2 and CD-QDPT2, is reported. For the DF/CD-QDPT2 methods, the Hose–Kaldor approach is used. The DF-QDPT2 method, with the cc-pwCVTZ basis set, dramatically reduces the computational cost compared to conventional multiconfigurational QDPT2 (MC-QDPT2, from the Gamess 2017.R2 package), with a more than 122-fold reduction for the largest member of the diradical test set considered. The DF approximation enables substantially accelerated energies to be obtained for the QDPT2 approach due to the significantly reduced I/O time. The performance of the DF-QDPT2 and CD-QDPT2 methods is compared with that of CASSCF, the multireference second-order perturbation theory (MRMP2), MC-QDPT2, and CASSCF-based second-order perturbation theory (CASPT2) methods for singlet-triplet energy splitting (EST ) in O2 and C2 molecules and for the dissociation energy of F2 . For the O2 and C2 molecules, the performance of the DF-QDPT2 and CD-QDPT2 methods is significantly better than that of CASSCF, MRMP2, MC-QDPT2, and CASPT2; while for the F2 case, the results of DF-QDPT2, CD-QDPT2, MRMP2, MC-QDPT2, and CASPT2 are similar and remarkably better than that of CASSCF, which fails dramatically. Moreover, the DF-QDPT2, CASSCF, CASPT2, and MRCI+Q methods are applied to potential energy curves (PECs) for N2 , CH4 , and F2 molecules. Our results demonstrate that the performance of DF-QDPT2 is substantially better than that of CASSCF, and comparable with that of CASPT2 for the molecules considered. Overall, the present application results demonstrate that the DF-QDPT2 and CD-QDPT2

2

ACS Paragon Plus Environment

Page 2 of 61

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

methods are very promising for electronically challenging molecular systems suffering from (quasi)degeneracy problems, where the single-reference methods cannot provide an accurate electronic description, but the DF-QDPT2 and CD-QDPT2 methods can do so at significantly reduced computational costs.

3

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

The multiconfigurational self-consistent field (MCSCF) approach is frequently used to treat chemical systems that exhibit quasidegeneracies and therefore require a multireference wave function. 1–5 The most common variant of the MCSCF method is the complete active space self-consistent field (CASSCF) method. 5–7 Quasidegeneracies occur, for example, in the vicinity of conical intersections, during bond breaking, in free radical chemistry, excited states, and some transition metal compounds. Although MCSCF provides qualitatively correct zeroth-order wave functions, the insufficient treatment of dynamic correlation generally prevents quantitative agreement with experiment. Hence the need for multireference methods that include dynamic correlations, such as multireference configuration interaction (MRCI), 3,8–15 multireference coupled-cluster (MRCC), 16–22 and multireference perturbation theory (MRPT). 23 MRCI methods, such as MRCISD, are expensive, and they suffer from the size extensivity problem. Furthermore, even though MRCC methods are size extensive, their high computational cost limits their applications to relatively small chemical systems. However, MRPT methods are relatively inexpensiv,e and they are also size extensive. Many variants of MRPT methods are described in the literature. Among these variants, one of the most common and systematic models is the quasidegenerate perturbation theory (QDPT). 24–30,30–45 Especially, the multistate structure of QDPT makes it promising for the study of excited electronic states as well as ground states. In the QDPT2 approach, the Møller–Plesset partitioning is used, and it reduces to the MP2 method in the case of a single reference. The QDPT2 approach fulfills all the desirable features of the non-degenerate perturbation theory, such as being non-iterative and size extensive. In addition, the QDPT2 formulation can be used with the Hartree–Fock (HF) reference as well as the MCSCF reference. Other commonly used MRPT variants are the CASPT2 method of Roos and coworkers 46,47 and the MRMP2 method of Hirao. 48–51 Other variants are Mukherjee’s state-specific MRPT method (Mk-MRPT2), 52–56 n-electron valence perturbation theory (NEVPT2), 57–59 multireference Brillouin-Wigner perturbation theory (MRBWPT), 60–64 and the generalized 4

ACS Paragon Plus Environment

Page 4 of 61

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

van Vleck perturbation theory (GVVPT). 65–69 Tensor decomposition of the electron repulsion integrals (ERIs) has attracted significant interest in modern computational chemistry. 70–93 The most common ERI factorization method is density fitting (DF). 70–77,84–93 In the DF approach, the four-dimensional ERI tensor is expressed in terms of three-dimensional tensors. Another popular integral approximation method is the partial Cholesky decomposition (CD). 80–85,88,89 The DF and CD approximations are very useful in reducing the storage requirements of molecular integrals as well as the cost of integral transformations. In the DF approach, ERIs are expanded in terms of a predefined auxiliary basis set, while in the CD method, the Cholesky factors are formed from the atomic orbital (AO) integrals. Although second-order MRPT methods (MRPT2) are not as expensive as the MRCI and MRCC methods, their computational costs are still significantly higher than that of singlereference MP2. The high cost of the common MRPT2 approaches compared with MP2 arises from the expensive CASSCF orbital optimization step. Further, the usage of conventional four-index ERIs prevents their application to larger molecular systems due to expensive I/O procedures. To address these bottlenecks of the conventional MC-QDPT2 approach, in this study, an efficient implementation of QDPT2 with the DF and CD approximations, denoted by DF-QDPT2 and CD-QDPT2, is reported. For the DF/CD-QDPT2 methods, the Hose– Kaldor approach is used. Furthermore, the present implementation employs Hartree–Fock orbitals directly to avoid the expensive CASSCF orbital optimization procedure. In addition to DF/CD-QDPT2, we also implement the DF-CAS-QDPT2 approach, for which DFCASSCF orbitals (from Psi4) are employed instead of HF orbitals. The presented equations have been coded in the C++ language, using unrestricted formalism, as an external plugin to the Psi4 package. 94 The computational time of the DF-QDPT2 energies is compared with that of the conventional QDPT2 method (from Gamess 2017.R2). 95,96 The DF/CD-QDPT2 methods are applied to chemical systems with challenging electronic structures.

5

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

2

Page 6 of 61

Theoretical Approach

2.1

Quasidegenerate Perturbation Theory

Let us assume that we know the solution of the Schrodinger equation (SE) for the unperturbed Hamiltonian.

ˆ 0 Φn = E (0) Φn , H n

(1)

ˆΨ e n = En Ψ e n. H

(2)

and the SE for the exact problem is

The Hamiltonian operator can be partitioned as follows:

ˆ = H ˆ 0 + Vˆ , H

(3)

where Vˆ is the perturbation operator. In QDPT, the space spanned by the reference functions Φν is called the model space (MS), which is the quasidegenerate set of zero-order states. To obtain working equations of QDPT, we define the projection operators in terms of the orthonormal zero-order functions 43 as follows:

Pˆ =

X

|Φν ihΦν | =

ν

X

Pˆν ,

(4)

ν

where Pˆ is the projector onto the MS. The orthogonal complement of Pˆ is defined as follows:

ˆ = ˆ1 − Pˆ = Q

X

|ΦI ihΦI |,

(5)

I ∈M / S

where the indices µ, ν, . . . sum over the MS functions and the indices I, J, . . . over all others. 6

ACS Paragon Plus Environment

Page 7 of 61 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

ˆ is the wave operator (WO), which generates the space spanned by the In the QDPT, Ω perturbed wave functions when operating on the MS.

ˆ ν = Ω ˆ Pˆ Φν . Ψν = ΩΦ

(6)

ˆ 0 . However, eigenfunctions of the The model state functions Φν are eigenfunctions of H ˆ are expected to contain sizable contributions from more than one model full Hamiltonian H function because of their (quasi)degeneracy. Using intermediate normalization, we may express the WO as follows:

ˆ = Pˆ + Q ˆ Ω. ˆ Ω

(7)

e ν of H ˆ because it just adds contribuThe WO cannot directly generate the eigenfunctions Ψ tions from the orthogonal space to each model function. Therefore, a final transformation eν: among the perturbed functions Ψν is required to generate Ψ

eν = Ψ

X

X

cµν Ψµ =

ˆ µ = Ω ˆΦ eν, cµν ΩΦ

(8)

µ

µ

where

eν = Φ

X

cµν Φµ .

(9)

µ

The energies and the transformation matrix c are obtained by a final diagonalization of an effective Hamiltonian. In QDP,T we can write the following SE for the exact problem:

ˆ ef f Φ e ν = Eν Φ eν, H

7

ACS Paragon Plus Environment

(10)

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

Page 8 of 61

ˆ ef f is the effective Hamiltonian: where H

ˆ ef f = Pˆ H ˆ Ω. ˆ H

(11)

The effective Hamiltonian can be partitioned as follows:

ˆ ef f = H ˆ 0 Pˆ + W ˆ, H

(12)

ˆ is called the level-shift operator (LSO). where W We can define the cumulative second-order effective Hamiltonian matrix as follows:

ef f (2) Hµν = Eν(0) δνµ + Vµν +

X

VµI RνI VIν ,

(13)

I

where

RνI =

1 1 = (0) (0) DνI Eν − EI

(2)

(14)

(2)

The energies (Eν ) and the transformation coefficients (cµν ) for the second-order QDPT are obtained by solving the non-Hermitian matrix eigenvalue problem.

Heff (2) c(2) = c(2) E(2) .

2.2

(15)

Normal Ordering for Perturbation Operator and LSO

To take advantage of the generalized Wick’s theorem and the diagrammatic techniques, we use the normal-ordered perturbation operator, which is defined as,

VˆN = Vˆ − h0|Vˆ |0i.

8

ACS Paragon Plus Environment

(16)

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

It is also convenient to define a modified LSO as follows: 43

ˆ0 ≡ W ˆ − h0|Vˆ |0i, W

(17)

ˆ , only the first-order part is modified. 43 where, in the perturbative series expansion for W The effective Hamiltonian can then be written as.

ˆ ef f = H ˆ 0 + h0|Vˆ |0i + W ˆ 0(1) + W ˆ (2) + . . . H

(18)

ˆ (1) are With these definitions, the generalized Block equation remains valid if Vˆ and W ˆ 0(1) . 43 replaced by VˆN and W

2.3

The Hose–Kaldor Approach

An alternative formalism to Brandow’s original formulation was introduced by Hose and Kaldor. 27–29 In the Hose–Kaldor (HK) approach, different Fermi vacuum states (reference states) are used for the computation of the different columns of the WO and LSO matrices. The ket state is used as the Fermi vacuum for all matrix elements in the same column. Hence, different Fermi vacuums are used for the different columns of the WO and LSO matrices. In a 1981 study, Jeziorski and Monkhorst 16 suggested a similar approach for the multireference coupled-cluster theory. In presenting this formalism, we use the diagrammatic approach of Shavitt and Bartlett. 43 In the HK approach, the molecular orbital (MO) space is divided into core, valence (active), and virtual spaces. However, for each different Fermi vacuum, the active orbitals are classified as active holes and active particles. For each reference state, the hole orbitals consist of the core plus the active holes of that reference, and the particle orbitals consist of the virtual orbitals plus the active particles. This classification of the MO space is depicted in Figure 1. However, the active orbitals in this figure do not correspond to the global energy ordering of the relevant MOs, but the ordering of the active MOs is based on the MS used as 9

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

Page 10 of 61

the reference. In the diagrammatic approach, active orbitals are shown by double arrows. 43 In Figure 1, the following notation is used for the orbital labels: i, j, k, . . . for holes; a, b, c, . . . for particles; i0 , j 0 , k 0 , . . . for valence (active) holes; a0 , b0 , c0 , . . . for valence (active) particles; i00 , j 00 , k 00 , . . . for core (inactive holes); and a00 , b00 , c00 , . . . for virtual (inactive particles).

2.4

The One-Electron Interaction

The cancellation between the two-electron part of the normal-ordered Fock matrix (UˆN ) and VˆN operators that simplifies the expression for the single reference perturbation theory and coupled-cluster methods does not hold in the HK formalism of QDPT. 40,43 The reason for this is that due to the varying definition of the reference state for the different columns of the WO and LSO matrices, a different definition of the occupied orbitals is used to form the Fock and perturbation matrices. The two-electron part of the Fock matrix is defined P as upq = i hpi||qii, where the summation over i runs over the occupied spin-orbitals in the determinant used as the reference state for that Fock operator. The common choice for the reference determinant used to build the Fock matrix is just the core determinant, which includes only core orbitals. However, one of the components of the two-electron part of the P 0 = i hpi||qii, where the summation over i runs over the perturbation matrix is given by vpq occupied spin-orbitals in the current Fermi vacuum state, which includes the present active holes as well as the core orbitals. The one-electron part of the normal-ordered perturbation operator in the HK formalism is written as. 43

ˆ N = Fˆ o + Vˆ 0 − UˆN = G N N

X

gpq {ˆ p† qˆ},

(19)

p,q

where FˆNo is the off-diagonal part of the normal-ordered Fock matrix. If we use the core state

10

ACS Paragon Plus Environment

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

to build the Fock matrix, then gpq can be written as. 43 o gpq = fpq +

X

hpi0 ||qi0 iDF ,

(20)

i0

where hpq||rsiDF are the antisymmetrized two-electron integrals within the DF/CD approximations, and the sum over i0 is running over the active occupied orbitals for the current Fermi vacuum. It should be noted that different gpq matrices are used for the different columns of ˆ operator will be denoted by the vertex • − − − ⊗ in the the WO and LSO matrices. The G diagrams (see Supporting Information).

2.5

First-Order Diagrams

The zeroth-order energy for a model state can be written as. X

E (0) =

εi ,

(21)

i

where εi are the orbital energies, which are equal to the diagonal elements of the Fock matrix (fii ). However, the range of summation depends on the Fermi vacuum used. Hence, the zeroth-order energies are different for different reference states. The first-order energy can be written as follows:

E (1) = h0|Vˆ |0i =

X i

gii −

1X hij||ijiDF . 2 i,j

(22)

ˆ 0(1) , which is equal The remaining part of the first-order LSO is the modified operator W to the normal-ordered perturbation operator (VˆN ). Hence, its diagonal elements are zero. The off-diagonal part of the first-order LSO is represented by two diagrams. The first one corresponds to the following algebraic equation: 0 ˆ 0(1) |0i = ga0 i0 , hΦai0 |W

11

ACS Paragon Plus Environment

(23)

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

Page 12 of 61

which corresponds to a single replacement with respect to the Fermi vacuum (the ket state). Similarly, the second diagram can be written as 0 0

ˆ 0(1) |0i = ha0 b0 ||i0 j 0 iDF , hΦia0 jb0 |W

(24)

which corresponds to a double replacement with respect to the Fermi vacuum (the ket state). 0

In these equations, |0i represents the current Fermi vacuum (the ket state), and |Φai0 i and 0 0

|Φai0 jb0 i are written with respect to the ket state.

2.6

Second-Order Diagrams

Skeleton diagrams for the second-order LSO were provided by Shavitt and Bartlett. 43 However, they provide neither explicit diagrams nor explicit algebraic equations. Starting from the skeletons provided by Shavitt and Bartlett, we generate the explicit diagrams and convert them to algebraic equations. Our second-order LSO diagrams are provided in Supporting Information. In our diagrams, indices presented in parentheses, such as (a), indicate dummy indices, while the other indices are free indices. Further, because the WO includes the resolˆ vent operator, which includes the Q-space functions, summations in the WO diagrams must exclude MS functions. Therefore, in the LSO diagrams, the summations over intermediate states must exclude MS functions. 43 In other words, all indices for the denominators Dia ab cannot be active orbitals simultaneously; at least one of them should be inactive. and Dij

Moreover, because we are using a complete model state (CMS) approach, only connected diagrams are considered. Before presenting the second-order diagrams, let us introduce the first-order amplitudes as follows:

tai Dia = gai , ab tab = hab||ijiDF , ij Dij

12

ACS Paragon Plus Environment

(25) (26)

Page 13 of 61 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

where

Dia = fii − faa − ∆, ab Dij = fii + fjj − faa − fbb − ∆,

(27) (28)

where ∆ is the level-shift parameter, which should be added to the denominators if so-called intruder state problems occur. In our implementation, we use ∆ = 0 as the default value. However, users may provide a positive value for ∆, typically between 0.2 and 0.5 a.u., to enforce level shifting. We have the following terms, starting with the diagonal term:

ˆ (2) |0i = h0|W

X 1 X ab tij hij||abiDF + tai gia . 4 ijab ia

(29)

Single replacements: 0

1 X ab 0 t 0 ha i||abiDF 2 iab i i 1 X a0 a − t hij||i0 aiDF 2 ija ij X 0 + tai0 ia gia

ˆ (2) |0i = hΦai0 |W

ia

+

X

+

X



X

tai ha0 i||i0 aiDF

ia

tai0 ga0 a

a 0

tai gii0 .

i

13

ACS Paragon Plus Environment

(30)

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

Page 14 of 61

Double replacements: 0 0

1 X ab 0 0 t 0 0 ha b ||abiDF 2 ab i j 1 X a0 b0 t hij||i0 j 0 iDF + 2 ij ij X 0 0 0 − P− (i0 j 0 )P− (a0 b0 ) tab i0 i ha i||aj iDF

ˆ (2) |0i = hΦai0 jb0 |W

ia 0 0

X

0 0

X

0 0

X

0 0

X

+ P− (a b )

0 tai0 ja0 gb0 a

a

− P− (i j )

0 0

tai0 ib gij 0

i

+ P− (i j )

taj0 ha0 b0 ||i0 aiDF

a

− P− (a b )

0

tbi ha0 i||i0 j 0 iDF .

(31)

i

Triple replacements: 0 0 0 ˆ (2) |0i = P (i0 j 0 /k 0 )P (a0 /b0 c0 ) hΦia0 jb0 kc0 |W

X

0

tai0 ja0 hb0 c0 ||ak 0 iDF

a 0

0 0

0 0

0

− P (i /j k )P (a b /c )

X

0 0

tai0 ib hic0 ||j 0 k 0 iDF ,

(32)

i

where

P− (ij) = 1 − P(ij),

(33)

P (ij/k) = 1 − P(ik) − P(jk),

(34)

P (i/jk) = 1 − P(ij) − P(ik),

(35)

where P(ij) permutes the indices i and j.

14

ACS Paragon Plus Environment

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

2.7

Approximate hSˆ2 i Value for the QDPT2 Wave Function

Because in our QDPT approach we employ the SD basis, different roots of the effective Hamiltonian may include different spin states. For example, in a CMS consisting of two orbitals, ˆ ef f will be one α-electron, and one β-electron, which is denoted by CMS(211), one root of H approximately a triplet state, while others will be approximately singlet states. Therefore, monitoring the spin multiplicity for the considered roots is important. Computation of the exact hSˆ2 i value for the QDPT2 wave function is possible, but it is expensive. Instead, we may compute an approximate value. For the ith root of the effective Hamiltonian, hSˆ2 i can be approximated as follows:

hΨi |Sˆ2 |Ψi i =

MS X

cµi cνi hΦµ |Sˆ2 |Φν i.

(36)

µν

In the HK approach, we always choose the ket state as the reference state. Using the secondquantization formalism for the Sˆ2 operator, 87 we obtain the following matrix elements, between the Slater determinants, for the diagonal and off-diagonal terms, respectively: occ X   1 α 1 α β 2 β 2 ˆ N −N + N +N − SIj SjI , hΦν |S |Φν i = 4 2 Ij

hΦai |Sˆ2 |0i = −

X

(37)

SaJ SJi ,

(38)

SjI SAj ,

(39)

J

ˆ2 hΦA I |S |0i = −

X j

ˆ2 hΦbA jI |S |0i

= −SbI SAj ,

where S is the overlap matrix.

15

ACS Paragon Plus Environment

(40)

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

2.8

Page 16 of 61

Formation of the MS

To generate the MS functions, we employ Rettrup’s graphical representation, 97 which provides a straightforward and efficient way to address the MS functions. In the graphical representation, the Slater determinants are represented by α and β strings. 10,98,99 Rettrup’s method is similar to Olsen’s graphical method. 99 However, in Rettrup’s method, one needs to compute only the vertex weights (VWs), while in Olsen’s case, one needs to compute both arc weights and VWs. Another attractive feature of Rettrup’s method is that it is possible to readily compute the index number from the configuration and vice versa. We have slightly modified Rettrup’s algorithm to adapt it to the C++ language. We compute the VWs as binomial coefficients:

Wk,i =

   

k! (i+1)!(k−i−1)!

  0

if k > i else

where the first dimension of the matrix W is the number of active orbitals (for example, α-orbitals) and the second dimension is the number of active electrons (for example, αelectrons). The index number for a given configuration (Slater determinant) K can then be computed as follows:

I(K) =

Ne X

Wk,i .

(41)

i

2.9

ˆ ef f Diagonalization of H

For the diagonalization of the effective Hamiltonian, we use two algorithms. If there is enough memory, we choose full diagonalization, which is feasible for relatively small active spaces. ˆ ef f in the core memory, then we use the iterative If there is not enough memory to store H Davidson algorithm. 10,100–102 However, we slightly modify the original Davidson algorithm.

16

ACS Paragon Plus Environment

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

In the original algorithm, we use the CI Hamiltonian and define

Gij = hbi |H|bj i = hbi |σj i,

(42)

where {bi } is a set of orthonormal guess eigenvectors, and σj is the transformed Hamiltonian. However, in the case of QDPT2, because we use the ket state as the reference state (HK ˆ ef f is a non-symmetric matrix, approach), we would like to sum over bra indices. Because H we slightly modify the definition of σ to the following:

Gij = hbi |Hef f |bj i = hσi |bj i.

2.10

(43)

Comparison of QDPT and MCSCF Methods

At first, let us consider the matrix form of the effective Hamiltonian, which can be written as follows:

ef f (2) Hµν = Hµν + Wµν + ....

(44)

In the MCSCF method, the Hµν matrix is diagonalized, while in the QDPT method the ef f Hµν matrix is diagonalized. If we assume that the MCSCF wave function includes only (2)

non-dynamical correlation, then we can state that the presence of Wµν and higher-order matrices introduces dynamical correlation into the wave function. Therefore, we may define the dynamical correlation operator as follows:

ˆ dc = W ˆ (2) + . . . ; W

(45)

ˆ ef f = Pˆ H ˆ Pˆ + W ˆ dc . H

(46)

then

17

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

Page 18 of 61

Further, when all the orbitals and electrons are chosen to be active, which corresponds to full configuration interaction (FCI), the effective Hamiltonian reduces to the full Hamiltonian ˆ dc vanishes: because FCI gives the exact solution. Hence, the operator W

ˆ ef f = Pˆ H ˆ Pˆ , (if all active) H ˆ dc = 0, (if all active). W

(47) (48)

ˆ dc operator introduces the missing dynamical correlation This situation indicates that the W due to the limited active space, and it vanishes when all the MOs and electrons are active. Further, the FCI energy is invariant with respect to any orbital rotations, and hence one can use the HF orbitals directly; orbital optimization is not necessary for the FCI method. If all the dynamical correlation is completely described in the Hamiltonian that will be diagonalized then there is no need for orbital optimization, and the HF orbitals can be directly used. Hence, in the case of limited active spaces, the QDPT wave function will be less sensitive to the orbital relaxation effects compared to the MCSCF case, because the QDPT effective Hamiltonian includes the dynamical correlation operator. Further, we expect that the amount of dynamical correlation is increased in the effective Hamiltonian that is used, the orbital effects will be less important. On the basis of this discussion, in this study, we will use the HF orbitals directly to generate our multiconfigurational wave function. However, we note that the use of the HF orbitals instead of the CASSCF orbitals may lead to higher absolute energies. Nevertheless, in most cases, we interest in energy differences rather than absolute energies.

18

ACS Paragon Plus Environment

Page 19 of 61 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

2.11

ERI Decomposition Techniques

The AO basis integrals can be cast into the following form with the help of the DF and CD approaches:

(µν|λσ)DF =

N aux X

Q bQ µν bλσ ,

(49)

Q

In the CD method, the CD tensors bQ µν are generated during the partial CD procedure, and Q is a Cholesky vector, while the DF tensors bQ µν are defined as follows:

bQ µν

N aux X

=

(µν|P )[J−1/2 ]P Q ,

(50)

P

where Z Z (µν|P ) =

χµ (r1 )χν (r1 )

1 ϕP (r2 ) dr1 dr2 , r12

(51)

and Z Z JP Q =

ϕP (r1 )

1 ϕQ (r2 ) dr1 dr2 , r12

(52)

where χµ (r) and ϕP (r) are the primary and auxiliary basis set members, respectively. Similarly, the two-electron integrals in the MO basis may be written as follows:

(pq|rs)DF =

N aux X

Q bQ pq brs .

(53)

Q

where bQ pq is an MO basis DF/CD tensor. In our present study, all the necessary molecular integrals are generated on the fly from the DF/CD factors, which significantly reduces the I/O time. Further, in the HK approach, one needs to sort the two-electron integrals (TEIs) when looping over the ket states, because

19

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

Page 20 of 61

different Fermi vacuums are used for the different columns of the WO and LSO matrices. In the case of conventional TEIs, the sorting procedure is performed on the disk, which substantially slow downs the QDPT2 method. Hence, conventional TEIs limit application of the HK approach to small-sized chemical systems. However, with the DF/CD approaches, the sorting procedure is performed in the core memory, and it scales only as O(N 3 ), where N is the number of basis functions. Therefore, with the help of DF/CD approach, one can apply the HK approach to large-sized chemical systems. Consequently, the DF/CD approaches have been utilized for integral evaluation, transformation, sorting procedures, and diagonal elements evaluation, but not for defining new intermediates for off-diagonal elements. With the help of the DF/CD approach, we can evaluate the first-order amplitudes as follows:

ab tab ij Dij

=

N aux X



 Q Q Q bQ b − b b ia jb ib ja ,

(54)

Q

and hence, we directly compute the tab ij amplitudes from the DF/CD factors. Furthermore, the diagonal term can be written as follows: Naux X X 1X Q (2) ˆ bQ t + tai gia , h0|W |0i = 2 Q ia ia ia ia

(55)

where

tQ ia

=

occ X vir X j

20

Q tab ij bjb .

b

ACS Paragon Plus Environment

(56)

Page 21 of 61 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

3

Results and Discussion

Results from the DF-QDPT2, CASSCF (from Psi4), and CAS-QDPT2 methods were obtained for a set of diradicals for assessment of the computational cost for single-point energy computations. For the diradicals, the cc-pwCVTZ basis set was employed. 103–105 The DF-QDPT2, CASSCF, MRMP2, 48–51 MC-QDPT2, 41 and CASPT2 46,47 methods were applied to the singlet-triplet energy splitting for the O2 and C2 molecules. Further, the dissociation energy of the F2 molecule was also considered to assess the performance of our DF-QDPT2 method. These computations were performed with Dunning’s correlationconsistent polarized core and valence triple-ζ basis sets (cc-pwCVTZ and aug-cc-pwCVTZ ). 103–105 For the cc-pwCVTZ/aug-cc-pwCVTZ primary basis sets, def2-QZVPP-JKFIT 75,106 and cc-pwCVTZ-RI/aug-cc-pwCVTZ-RI 107 auxiliary basis set pairs were used for reference and correlation energies, respectively. For the CASSCF and CASPT2 methods, geometry optimizations and frequency computations were performed via analytic derivatives, while for the other methods, computations were performed via numerical derivatives. For diatomic molecules, bond distances were obtained using the diatomic 108 program. For each method five single point energies, uniformly distributed about the equilibrium bond length (±0.005 Å, ±0.010 Å), were provided to the diatomic program. Moreover, the DF-QDPT2, CASSCF, CASPT2, and MRCISD with Davidson correction (MRCI+Q) 3,8,9 methods were applied to potential energy surfaces (PESs) for the N2 , CH4 , and F2 molecules. In addition, to demonstrate the applicability of the DF-QDPT2 approach for large-sized chemical systems, we consider the singlet-triplet energy splittings for the diarylethene derivatives, which are well-known photochromophores. 109–111 DF-QDPT2, CD-QDPT2, CAS-QDPT2, and DF-CAS-QDPT2 computations were performed with our present code, while the conventional MC-QDPT2 and MRMP2 computations were carried out with the Gamess 2017.R2 program, 95,96 and CASPT2 with the Molpro2012.1 program. 112,113 Finally, we note that with our determinantal multistate QDPT2 approach, we can study singlet and triplet states in the same subspace because both states 21

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

have MS = 0 components. Thus, throughout this study for our DF/CD-QDPT2 methods, all triplet energies are obtained as low-spin triplets. Hence, HF orbitals for singlet states are always employed for DF/CD-QDPT2.

3.1

Efficiency of the DF-QDPT2 Method

We consider a set of diradicals to assess the efficiency of the DF-QDPT2 method. The chemical structures of the diradicals considered are depicted in Figure 2. The total computational (wall) times for single-point energies using the CASSCF (from Psi4), DF-QDPT2, and CASQDPT2 methods, where conventional CASSCF orbitals are used for our QDPT2 instead of the HF orbitals, are presented in Figure 3. These computations were performed on a single core of an Intel(R) Xeon(R) CPU E5-2620 v3 @ 2.40 GHz computer. For each method, ∼64 GB memory is provided. For all diradical computations, a CMS consisting of two orbitals, one α-electron, and one β-electron [CMS(211)], which is commonly denoted by CAS(2,2), is considered. For the largest member of the diradical set, DR13, the numbers of occupied and virtual orbitals are O = 35 and V = 466. As shown in Figure 3, the DF-QDPT2 method dramatically reduces the total computational cost compared to conventional CAS-QDPT2; there is a 72-fold reduction in wall time compared to CAS-QDPT2 for the largest member of the diradical test set. Most of this cost difference arises from the CASSCF orbital optimization procedure. Our DF-QDPT2 method is 47-fold faster than the CASSCF method for the largest member. Hence, due to the dramatically reduced I/O time, the DF approach enables substantially accelerated QDPT2 computations. The number of CASSCF iterations is 18 for DR7, while it is 8–12 for the remaining diradicals. The average number of CASSCF iterations is ∼ 10 for the test set considered. Furthermore, we compare the computational time of DF-QDPT2 with the MC-QDPT2 (from Gamess 2017.R2) methods for the largest member of the diradicals considered, DR13. For the diradical DR13, the computational times are 10.53 (DF-QDPT2) and 1288.77 min (MC-QDPT2). Hence, there is a 122-fold reduction in wall time compared to MC-QDPT2. In 22

ACS Paragon Plus Environment

Page 22 of 61

Page 23 of 61 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

our DF-QDPT2 computations, we first perform an HF computation, after which we generate the MS (active space) determinants and compute the QDPT2 effective Hamiltonian. Finally, we obtain the determinant coefficients and the state energies from the diagonalization of ˆ ef f . However, in the MC-QDPT2 method of Nakano, 41 a CASSCF computation is first H ˆ ef f matrix is then formed in the MCSCF basis, and it is diagonalized. performed. The H ˆ ef f in the Hence, our DF-QDPT2 code directly uses the HF orbitals and diagonalizes the H space spanned by the Slater determinants forming the MS. Another important difference is that the MC-QDPT2 program uses configuration state functions (CSFs), whereas we use the determinant basis. In addition, we use the DF integrals, whereas the MC-QDPT2 code uses the conventional integrals.

3.2 3.2.1

Accuracy of the DF-QDPT2 Method Singlet-Triplet Energy Splitting in O2

To assess the accuracy of the DF/CD-QDPT2 method, we start with the oxygen molecule, for which the X 3 Σg and a 1 ∆g states are considered. Table 1 presents the total energies, equilibrium bond lengths, and singlet-triplet energy splittings (EST ) for the O2 molecule from the DF-QDPT2, CD-QDPT2, CASSCF, MRMP2, MC-QDPT2, and CASPT2 methods with the cc-pwCVTZ basis set. In addition, we also consider the DF-CAS-QDPT2 approach, for which DF-CASSCF orbitals (from Psi4) are employed instead of HF orbitals. For each method, a CMS consisting of two orbitals, one α-electron, and one β-electron is used [CMS(211)]. One may consider a larger CMS; however, on the basis of our previous discussion of our DF/CD-QDPT2 approaches, we may employ small CMSs. In defining the CMS, our strategy is to always start with CMS(211), which uses the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). We then investigate the HOMO-1 and LUMO+1 orbitals and determine if there is a degeneracy/quasidegeneracy with HOMO and LUMO, respectively. For O2 at the triplet geometry, re = 1.254 Å, the HF orbital energies are −0.630186 (HOMO-1), −0.485193 (HOMO), 0.009543 (LUMO), and 23

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

Page 24 of 61

0.364511 (LUMO+1) hartree. Similarly, at the singlet geometry, at re = 1.271 Å, the HF orbital energies are −0.622051 (HOMO-1), −0.491863 (HOMO), 0.002017 (LUMO), and 0.347611 (LUMO+1) hartree. Hence, a CMS(211) MS appears to be an appropriate choice for our purpose. For the ground state (X 3 Σg ) of the O2 molecule, the absolute errors in bond lengths (∆re ) are 0.046 (DF-QDPT2), 0.046 (CD-QDPT2), 0.045 (DF-CAS-QDPT2), 0.055 (CASSCF), 0.108 (MRMP2), 0.089 (MC-QDPT2), and 0.030 (CASPT2) Å. Hence, the CASPT2 method provides the lowest errors with respect to the experimental value of 1.208 Å. 114 The errors of DF-QDPT2, CD-QDPT2, and DF-CAS-QDPT2 are quite reasonable, while those of MRMP2 and MC-QDPT2 are noticeably higher compared with DF/CD-QDPT2: there is almost a twofold difference in the errors. Both methods (MRMP2 and MC-QDPT2) fail to improve upon CASSCF for the bond length. For the lowest singlet state (a 1 ∆g ) of the O2 molecule, the ∆re values are 0.055 (DF-QDPT2), 0.055 (CD-QDPT2), 0.113 (DF-CAS-QDPT2), 0.063 (CASSCF), 0.174 (MRMP2), 0.025 (MC-QDPT2), and 0.136 (CASPT2) Å. In this case, the MC-QDPT2 method provides the lowest error, while MRMP2 and CASPT2 provide the highest errors. Further, in both the singlet and triplet cases, the CASSCF method provides shorter bond lengths, while the other methods generally provide longer values compared with the experiment. 114 Next, we consider EST values for the DF-QDPT2, CD-QDPT2, DF-CAS-QDPT2, CASSCF, MRMP2, MC-QDPT2, and CASPT2 methods (Table 1). The absolute errors in the EST values (∆EST ) are 0.2 (DF-QDPT2), 0.1 (CD-QDPT2), 2.7 (DF-CAS-QDPT2), 6.7 (CASSCF), 1.0 (MRMP2), 2.6 (MC-QDPT2), and 8.0 (CASPT2) kcal mol-1 . Hence, the DF-QDPT2 and CD-QDPT2 methods perform significantly better than DF-CAS-QDPT2, MRMP2, and MCQDPT2. There are 5-, 15-, and 15-fold errors in the results of MRMP2, MC-QDPT2, and DF-CAS-QDPT2, respectively, compared with DF-QDPT2. Further, the results of CASSCF and CASPT2 exhibit significant errors of 6.7 and 8.0 kcal mol-1 , respectively.

24

ACS Paragon Plus Environment

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

3.2.2

Singlet-Triplet Energy Splitting in C2

As the second step of our assessment, we consider the dicarbon molecule, for which the X 1 Σg and a 3 Πu states are considered. Table 2 presents the total energies, equilibrium bond lengths, and the singlet-triplet energy splittings for the C2 molecule from the DFQDPT2, CD-QDPT2, DF-CAS-QDPT2, CASSCF, MRMP2, MC-QDPT2, and CASPT2 methods. In the case of C2 , the singlet-triplet energy splitting is very small (−2.0 kcal mol-1 ); therefore, for a better description of electronic states, one needs to use diffuse functions. Hence, we employ the aug-cc-pwCVTZ basis set in this case. For C2 at the singlet geometry, re = 1.259 Å, the HF orbital energies are −0.454215 (HOMO-1), −0.454215 (HOMO), −0.116027 (LUMO), and 0.056418 (LUMO+1) hartree. Similarly, at the triplet geometry, re = 1.284 Å, the HF orbital energies are −0.448726 (HOMO-1), −0.448726 (HOMO), −0.115955 (LUMO), and 0.056219 (LUMO+1) hartree. Hence, for each method, a CMS consisting of three orbitals, two α-electrons, and two β-electrons is used [CMS(322)] because there is a degeneracy between HOMO and HOMO-1. For the ground state (X 1 Σg ) of the C2 molecule, the ∆re values are 0.016 (DF-QDPT2), 0.017 (CD-QDPT2), 0.007 (DF-CAS-QDPT2), 0.045 (CASSCF), 0.028 (MRMP2), 0.026 (MC-QDPT2), and 0.144 (CASPT2) Å. Hence, the DF-CAS-QDPT2 method provides the lowest errors with respect to the experimental value of 1.243 Å. 114 The errors of MRMP2, MC-QDPT2, and especially CASPT2 are remarkably higher compared with DF/CD-QDPT2. However, in this case, both MRMP2 and MC-QDPT2 significantly improve upon CASSCF with respect to the bond length. For the triplet state (a 3 Πu ) of the C2 molecule, the ∆re values are 0.028 (DF-QDPT2), 0.028 (CD-QDPT2), 0.001 (DF-CAS-QDPT2), 0.140 (CASSCF), 0.092 (MRMP2), 0.091 (MC-QDPT2), and 0.056 (CASPT2) Å. For the excited state, the DF-CAS-QDPT2 method again provides the lowest error compared with MRMP2, MC-QDPT2, and CASPT2, whereas the results of DF/CD-QDPT2 are reasonably accurate. Next, we consider EST values for the DF-QDPT2, CD-QDPT2, DF-CAS-QDPT2, CASSCF, MRMP2, MC-QDPT2, and CASPT2 methods (Table 2). The ∆EST values are 1.4 (DF25

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

QDPT2), 1.4 (CD-QDPT2), 9.9 (DF-CAS-QDPT2), 11.5 (CASSCF), 27.5 (MRMP2), 26.4 (MC-QDPT2), and 18.9 (CASPT2) kcal mol-1 . Hence, the DF-QDPT2 and CD-QDPT2 methods perform dramatically better than MRMP2, MC-QDPT2, and CASPT2. There are 14-, 20-, and 19-fold errors in the results of CASPT2, MRMP2, and MC-QDPT2, respectively, compared with DF-QDPT2. Further, the result of CASSCF also exhibits significant errors of 11.5 kcal mol-1 .

3.2.3

Dissociation Energy of F2

As the final step of our assessment, we consider the dissociation energy of the F2 molecule, for which the ground state X 1 Σg is considered. It is well known the dissociation energy of F2 is a very challenging problem for computational methods, and most of the common approaches fail to produce the correct result. 56,115 Table 3 presents the total energies, equilibrium bond lengths, and dissociation energy (De ) for the F2 molecule from the DFQDPT2, CD-QDPT2, DF-CAS-QDPT2, CASSCF, MRMP2, MC-QDPT2, and CASPT2 methods with the cc-pwCVTZ basis set. For the F2 molecule, at re = 1.423 Å, the HF orbital energies are −0.806829 (HOMO-4), −0.806829 (HOMO-3), −0.745506 (HOMO-2), −0.666890 (HOMO-1), −0.666890 (HOMO), 0.085496 (LUMO), and 0.833879 (LUMO+1) hartree. Hence, HOMO and HOMO-1 should be included in the MS because they are degenerate. However, at the dissociation limit, at re = 100.0 Å, the HF orbital energies are −0.728218 (HOMO-4), −0.728218 (HOMO-3), −0.728218 (HOMO-2), −0.728218 (HOMO1), −0.374739 (HOMO), −0.369447 (LUMO), and 0.894299 (LUMO+1) hartree. At the dissociation limit, HOMO-1 is degenerate with HOMO-2, HOMO-3, and HOMO-4. Hence, to satisfy the consistency of the MS, for the dissociation problem of F2 a CMS consisting of six orbitals, five α-electrons, and five β-electrons is used [CMS(655)]. Further, for the F2 molecule, we employ a tighter CD tolerance, 10−10 , instead of the default value of 10−6 . Unlike the DF-QDPT2 approach, the auxiliary basis set of CD-QDPT2 is not a pre-optimized basis set; instead, it is generated on the fly from the primary basis set, and the quality of the

26

ACS Paragon Plus Environment

Page 26 of 61

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

CD auxiliary basis set is controlled by the CD tolerance. We observe that the default CD tolerance is not tight enough to obtain the desired accuracy. Hence, we employed a tighter CD tolerance for the dissociation problem of F2 . For the ground state of the F2 molecule, the ∆re values are 0.011 (DF-QDPT2), 0.012 (CD-QDPT2), 0.008 (DF-CAS-QDPT2), 0.051 (CASSCF), 0.011 (MRMP2), 0.010 (MCQDPT2) Å, and 0.009 (CASPT2) Å. Hence, the performances of DF-QDPT2, CD-QDPT2, DF-CAS-QDPT2, MRMP2, MC-QDPT2, and CASPT2 are very similar, and all of these methods provide reasonable errors with respect to the experimental value of 1.412 Å. 114 Further, the error of the CASSCF method is remarkably higher compared with the other methods considered. Next, we consider the De values for the DF-QDPT2, CD-QDPT2, DF-CAS-QDPT2, CASSCF, MRMP2, MC-QDPT2, and CASPT2 methods (Table 3). The absolute errors in the De values (∆De ) are 1.9 (DF-QDPT2), 1.9 (CD-QDPT2), 0.3 (DFQDPT2), 20.5 (CASSCF), 1.7 (MRMP2), 1.8 (MC-QDPT2), and 1.6 (CASPT2) kcal mol-1 . Hence, the result of DF-CAS-QDPT2 is in very good agreement with the experimental value of 38.3 kcal mol-1 . Further, the performances of DF-QDPT2, CD-QDPT2, MRMP2, MCQDPT2, and CASPT2 are again very similar and significantly better than that of CASSCF, which fails dramatically.

3.3 3.3.1

PESs N2

To assess the performance of DF-QDPT2 for PESs, we start with the N2 molecule. The total energies from the FCI, 116 DF-QDPT2, CASSCF, CASPT2, and MRCI+Q methods with the cc-pVDZ basis set are reported in Table 4. These energies are presented graphically in Figure 4. For each multireference method, a CMS consisting of four orbitals, two α-electrons, and two β-electrons is used [CMS(422)]. The results of the MRCI+Q method are very close to the FCI energies; the mean absolute error is 2.03 millihartree (mh) throughout the potential energy curve (PEC). The DF-QDPT2 27

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

and CASPT2 curves also closely follow the FCI curve. At the minimum energy geometry, RN N = 1.1208 Å, the DF-QDPT2 method covers 99.991 % of the CASPT2 energy, and the DF-QDPT2 and CASPT2 curves almost overlap. Further, because we do not have a full PES for N2 , we investigate the maximum absolute error (∆max ) for the methods considered. The ∆max values are 25.58 (DF-QDPT2), 400.43 (CASSCF), 16.42 (CASPT2), and 5.52 (MRCI+Q) mh. The nonparallelity error (NPE) results show that the performance of DFQDPT2 is significantly better than that of CASSCF, and its performance approaches that of CASPT2. Therefore, we conclude that the DF-QDPT2 method provides results with similar quality as those from CASPT2 at a significantly lower cost.

3.3.2

CH4

As the second step of our assessment of the performance of DF-QDPT2 for PESs, we consider the CH4 molecule. The PEC for the CH4 molecule was obtained by altering the distance of a single C–H bond (RCH ), while the remaining three C–H bonds are fixed to the equilibrium value of 1.086 Å, and the H–C–H angles are kept at the tetrahedral value of 109.471◦ . The total energies from the FCI, 117 DF-QDPT2, CASSCF, CASPT2, and MRCI+Q methods with the 6-31G(d) basis set are reported in Table 5. These energies are presented graphically in Figure 5. For each multireference method, a CMS consisting of four orbitals, three αelectrons, and three β-electrons is used [CMS(433)]. The energies of the MRCI+Q method are comparable with the FCI energies, and the mean absolute error is 4.19 mh throughout the PEC. The PECs of DF-QDPT2 and CASPT2 are also reasonably accurate. At the minimum energy geometry, RCH = 1.1 Å, the DFQDPT2 method covers 99.996 % of the CASPT2 energy, and the DF-QDPT2 and CASPT2 curves almost overlap. Further, because we have a large enough interval for RCH values, we may compute the C–H bond dissociation energies (BDEs) to assess the performance of the methods considered. The C–H BDE values are 108.3 (FCI), 107.9 (DF-QDPT2), 98.4 (CASSCF), 107.0 (CASPT2), and 108.4 (MRCI+Q) kcal mol-1 . Hence, the absolute

28

ACS Paragon Plus Environment

Page 28 of 61

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

errors with respect to FCI are 0.42 (DF-QDPT2), 9.8 (CASSCF), 1.3 (CASPT2), and 0.1 (MRCI+Q) kcal mol-1 . These results demonstrate that the performances of DF-QDPT2 and CASPT2 are very similar, that of DF-QDPT2 is slightly better, while the BDE of CASSCF exhibits significant errors.

3.3.3

F2

As the final step of our assessment of the performance of DF-QDPT2 for PESs, we consider the F2 molecule. For this molecule, the FCI energies are not available, and hence we consider MRCI+Q as our reference method. The total energies from the DF-QDPT2, CASSCF, CASPT2, and MRCI+Q methods with the cc-pwCVTZ basis set are reported in Table 6. These energies are presented graphically in Figure 6. For each multireference method, a CMS consisting of six orbitals, five α-electrons, and five β-electrons is used [CMS(655)]. The PECs of CASPT2 and MRCI+Q are reasonably close to each other, while that of DF-QDPT2 is parallel to the MRCI+Q curve. At the minimum energy geometry, RF F = 1.4 Å, the DF-QDPT2 method covers 99.973 % of the CASPT2 energy. Further, because we have a large enough interval for the RF F values, we may compute the dissociation energies to assess the performance of the methods considered. The computed dissociation energies are 36.0 (DF-QDPT2), 16.8 (CASSCF), 35.1 (CASPT2), and 34.8 (MRCI+Q) kcal mol-1 . Hence, the absolute errors with respect to MRCI+Q are 1.2 (DF-QDPT2), 18.0 (CASSCF), and 0.2 (CASPT2) kcal mol-1 . These results demonstrate that the performances of both DFQDPT2 and CASPT2 are significantly accurate, with CASPT2 being more accurate than DF-QDPT2, while that of CASSCF again exhibits significant errors.

3.4

Large-Scale Applications

To demonstrate the applicability of the DF-QDPT2 approach for large-sized chemical systems, we consider the singlet-triplet energy splittings for the diarylethene derivatives, which are well-known photochromophores. 109–111 We consider the open and closed forms of 1,229

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

Page 30 of 61

bis(2-methyl-5-phenyl-3-thienyl)perfluorocyclopentene 109 and 1,2-bis(3-methyl-5-phenyl-2-thienyl)perfluoro The former is denoted by N-diarylethene (C27 H18 S2 F6 ) and the latter by I-diarylethene (C27 H18 S2 F6 ). Their chemical structures are shown in Figure 7. The optimized structures of N-diarylethene and I-diarylethene were obtained from Yanai et al. 111 The 6-311G(2d,2p) basis set is employed for diarylethenes, which includes 1053 basis functions (Cartesian functions are used). Table 7 presents the singlet-triplet energy splittings for the closed and open forms of the diarylethene derivatives from the DF-QDPT2 method. For the diarylethene derivatives considered, a CMS consisting of two orbitals, one α-electron, and one β-electron is used [CMS(211)]. The computed singlet-triplet energy splittings are 65.0 [N-diarylethene (open)], 24.6 [N-diarylethene (closed)], 54.1 [I-diarylethene (open)], and 38.1 [I-diarylethene (closed)] kcal mol-1 . These computations were completed in just ∼5 h in our Linux cluster. The same computations take several days with CASPT2, while it is not feasible to accomplish it with the Gamess MC-QDPT2 code. Hence, these results demonstrate that with our DF-QDPT2 approach, one can routinely study large-sized chemical systems.

4

Conclusions

In this research, an efficient implementation of the second-order quasidegenerate perturbation theory with the DF and CD approximations, which are denoted by DF-QDPT2 and CDQDPT2, has been reported. The computational time of the DF-QDPT2 single point energies has been compared with that of the conventional MC-QDPT2 method (from the Gamess 2017.R2 package). DF-QDPT2 dramatically reduces the computational cost compared to MC-QDPT2, with a more than 122-fold reduction for the largest member of the diradical test set considered, in a cc-pwCVTZ basis set. Due to the dramatically reduced I/O time, the DF approach enables substantially accelerated energies for QDPT2. However, we note that the use of the HF orbitals instead of the CASSCF orbitals may lead to higher absolute

30

ACS Paragon Plus Environment

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

energies. Nevertheless, in most cases, we interest in energy differences rather than absolute energies. DF-QDPT2 and CD-QDPT2 have been applied for singlet-triplet energy splitting (EST ) in the O2 and C2 molecules, and they have been applied for the dissociation problem of F2 . For the O2 molecule, absolute errors in EST values (∆EST ) are 0.2 (DF-QDPT2), 0.1 (CD-QDPT2), 6.7 (CASSCF), 1.0 (MRMP2), 2.6 (MC-QDPT2), and 8.0 (CASPT2) kcal mol-1 . Hence, the performances of DF-QDPT2 and CD-QDPT2 are significantly better than those of MRMP2 and MC-QDPT2. There are 5- and 15-fold errors in the results of MRMP2 and MC-QDPT2 compared with DF-QDPT2. Further, the results of CASSCF and CASPT2 exhibit significant errors of 6.7 and 8.0 kcal mol-1 , respectively. For the C2 molecule, the ∆EST values are 1.4 (DF-QDPT2), 1.4 (CD-QDPT2), 11.5 (CASSCF), 27.5 (MRMP2), 26.4 (MC-QDPT2), and 18.9 (CASPT2) kcal mol-1 . Hence, the performances of DF-QDPT2 and CD-QDPT2 are dramatically better than those of MRMP2, MC-QDPT2, and CASPT2. There are 14-, 20-, and 19-fold errors in the results of CASPT2, MRMP2, and MC-QDPT2, respectively, compared with DF-QDPT2. Further, the result of CASSCF also exhibits significant errors of 11.5 kcal mol-1 . For the F2 molecule, the absolute errors in dissociation energies are 1.9 (DF-QDPT2), 1.9 (CD-QDPT2), 20.5 (CASSCF), 1.7 (MRMP2), 1.8 (MC-QDPT2), and 1.6 (CASPT2) kcal mol-1 . Hence, the performances of DF-QDPT2, CDQDPT2, MRMP2, MC-QDPT2, and CASPT2 are very similar and significantly better than that of CASSCF, which fails dramatically. Moreover, DF-QDPT2, CASSCF, CASPT2, and MRCI+Q have been applied to PECs for the N2 , CH4 , and F2 molecules. Our results demonstrate that the performance of DF-QDPT2 is substantially better than that of CASSCF for the molecules considered. Further, the PEC results also show that the DF-QDPT2 method provides similar quality results as CASPT2, but at a significantly lower cost. The presented methods, DF-QDPT2 and CD-QDPT2, avoid the CASSCF orbital optimization procedure and build a CAS-CI space with HF orbitals. The computational costs of our methods are dominated by the evaluation of the diagonal elements of the effective

31

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

matrix, which is equal to a single-reference DF-MP2 energy computation. Hence, the cost diagonal elements scale as O(N )5 , where N is the number of basis functions. However, one needs to compute the diagonal elements for each MS (active space) determinant. The total cost of the methods considered then scales as Ndets ∗ O(N )5 . Therefore, the computational cost of DF/CD-QDPT2 is identical to the cost of DF-MP2 times the number of determinants [Cost(DF − QDP T 2) ≈ Ndets ∗ Cost(DF − M P 2)]. For the molecules considered in this study, the computational cost of the DF-QDPPT2 method is significantly lower than that of CASPT2 because the CASSCF orbital optimization procedure is avoided. However, in the case of very large active spaces, the CASPT2 method may provide lower costs compared to DF-QDPT2 due to its internally contracted feature. Our results demonstrate that the DF-QDPT2 and CD-QDPT2 methods are very promising for electronically challenging molecular systems suffering from (quasi)degeneracy problems, where the single-reference methods cannot provide an accurate electronic description. Even though the CASSCF method considers non-dynamical correlation, it may still fail dramatically due to insufficient dynamical correlation. Further, we conclude that as more dynamical correlation is considered in the effective Hamiltonian, multireference approaches become less sensitive to orbital relaxation effects. Hence, one may keep using the HF orbitals and skip the expensive orbital optimization procedure. Consequently, considering both the computational efficiency and the accuracy of the DF-QDPT2 and CD-QDPT2 approaches, we conclude that these methods emerge as very useful tools for computational chemistry.

Supporting Information Available (1) Second-order diagrams. (2) Spin-factored equations for the DF/CD-QDPT2 method. This material is available free of charge via the Internet at http://pubs.acs.org/.

32

ACS Paragon Plus Environment

Page 32 of 61

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

Acknowledgement This research was supported by the Scientific and Technological Research Council of Turkey (TÜBİTAK-116Z506).

References (1) Dalgaard, E.; Jørgensen, P. Optimization of Orbitals for Multiconfigurational Reference States. J. Chem. Phys. 1978, 69, 3833–3844. (2) Jørgensen, P.; Olsen, J.; Yeager, D. L. Optimization and Characterization of a Multiconfigurational Self–Consistent Field (MCSCF) State. Adv. Chem. Phys. 1983, 54, 1–176. (3) Werner, H.-J. Matrix-Formulated Direct Multiconfiguration Self-Consistent Field and Multiconfiguration Reference Configuration-Interaction Methods. Adv. Chem. Phys. 1987, 69, 1–62. (4) Shepard, R. The Multiconfiguration Self-Consistent Field Method. Adv. Chem. Phys. 1987, 69, 63–200. (5) Roos, B. O. A Complete Active Space Scf Method (CASSCF) Using a Density Matrix Formulated Super-CI Approach. Adv. Chem. Phys. 1987, 69, 399–446. (6) Roos, B. O.; Taylor, P.; Siegbahn, P. E. A Complete Active Space SCF Method (CASSCF) Using a Density Matrix Formulated Super-CI Approach. Chem. Phys. 1980, 48, 157–173. (7) Schmidt, M. W.; Gordon, M. S. The Construction and Interpretation of MCSCF Wavefunctions. Annu. Rev. Phys. Chem. 1998, 49, 233–266. (8) Werner, H.-J.; Reinsch, E. The Self-Consistent Electron Pairs Method for Multiconfiguration Reference State Functions. J. Chem. Phys. 1982, 76, 3144–3156. 33

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

(9) Werner, H.-J.; Knowles, P. J. An Efficient Internally Contracted MulticonfigurationReference Configuration Interaction Method. J. Chem. Phys. 1988, 89, 5803–5814. (10) Sherrill, C. D. The Configuration Interaction Method: Advances in Highly Correlated Approaches. Adv. Quantum Chem. 1998, 34, 143–269. (11) Venkatnathan, A.; Szilva, A. B.; Walter, D.; Gdanitz, R. J.; Carter, E. A. Size Extensive Modification of Local Multireference Configuration Interaction. J. Chem. Phys. 2004, 120, 1693–1703. (12) Shamasundar, K. R.; Knizia, G.; Werner, H.-J. A New Internally Contracted MultiReference Configuration Interaction Method. J. Chem. Phys. 2011, 135, 054101. (13) Szalay, P. G.; Muller, T.; Gidofalvi, G.; Lischka, H.; Shepard, R. Multiconfiguration Self-Consistent Field and Multireference Configuration Interaction Methods and Applications. Chem. Rev. 2012, 112, 108–181. (14) Krisiloff, D. B.; Carter, E. A. Approximately Size Extensive Local Multireference Singles and Doubles Configuration Interaction. Phys. Chem. Chem. Phys. 2012, 14, 7710–7717. (15) Saitow, M.; Kurashige, Y.; Yanai, T. Multireference Configuration Interaction Theory Using Cumulant Reconstruction with Internal Contraction of Density Matrix Renormalization Group Wave Function. J. Chem. Phys. 2013, 139, 044118–044133. (16) Jeziorski, B.; Monkhorst, H. J. Coupled-Cluster Method for Multideterminantal Reference States. Phys. Rev. A 1981, 24, 1668–1681. (17) Mahapatra, U. S.; Datta, B.; Mukherjee, D. A State-Specific Multi-Reference Coupled Cluster Formalism with Molecular Applications. Mol. Phys. 1998, 94, 157–171.

34

ACS Paragon Plus Environment

Page 34 of 61

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

(18) Mahapatra, U. S.; Datta, B.; Mukherjee, D. A Size-Consistent State-Specific Multireference Coupled Cluster Theory: Formal Developments and Molecular Applications. J. Chem. Phys. 1999, 110, 6171–6188. (19) Evangelista, F. A.; Allen, W. D.; Schaefer, H. F. High-Order Excitations in StateUniversal and State-Specific Multireference Coupled Cluster Theories: Model Systems. J. Chem. Phys. 2006, 125, 154113–154129. (20) Evangelista, F. A.; Allen, W. D.; Schaefer, H. F. Coupling Term Derivation and General Implementation of State-Specific Multireference Coupled Cluster Theories. J. Chem. Phys. 2007, 127, 024102–024119. (21) Evangelista, F. A.; Simmonett, A. C.; Allen, W. D.; Schaefer, H. F.; Gauss, J. Triple Excitations in State-Specific Multireference Coupled Cluster Theory: Application of Mk-MRCCSDT and Mk-MRCCSDT-n Methods to Model Systems. J. Chem. Phys. 2008, 128, 124104–124117. (22) Kállay, M.; Szalay, P. G.; Surjan, P. R. A General State-Selective Multireference Coupled-Cluster Algorithm. J. Chem. Phys. 2002, 117, 980–990. (23) Roskop, L.; Gordon, M. S. Quasi-Degenerate Second-Order Perturbation Theory for Occupation Restricted Multiple Active Space Self-Consistent Field Reference Functions. J. Chem. Phys. 2011, 135, 044101–044112. (24) Brandow, B. H. Linked-Cluster Expansions for the Nuclear Many-Body Problem. Rev. Mod. Phys. 1967, 39, 771–828. (25) Lindgren, I. A Coupled-Cluster Approach to the Many-Body Perturbation Theory for Open-Shell Systems. Int. J. Quantum Chem. 1978, 14, 33–58. (26) Brandow, B. H. Formal Theory of Effective Π-Electron Hamiltonians. Int. J. Quantum Chem. 1979, 15, 207–242. 35

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

(27) Hose, G.; Kaldor, U. Diagrammatic Many-Body Perturbation Theory for General Model Spaces. J. Phys. B. 1979, 12, 3827–3855. (28) Hose, G.; Kaldor, U. A General-Model-Space Diagrammatic Perturbation Theory. Phys. Scripta 1980, 21, 357–361. (29) Hose, G.; Kaldor, U. Quasidegenerate Perturbation Theory. J. Phys. Chem. 1982, 86, 2133–2140. (30) Freed, K. F.; Sheppard, M. G. Analysis of Ab Initio Effective Valence Shell Hamiltonian Calculations Using Third Order Quasidegenerate Many-Body Perturbation Theory. J. Chem. Phys. 1981, 75, 4507–4524. (31) Primas, H. Generalized Perturbation Theory in Operator Form. Rev. Mod. Phys. 1963, 35, 710–712. (32) Kirtman, B. Variational Form of Van Vleck Degenerate Perturbation Theory with Particular Application to Electronic Structure Problems. J. Chem. Phys. 1968, 49, 3890–3894. (33) Certain, P. R.; Hirschfelder, J. O. New Partitioning Perturbation Theory. I. General Formalism. J. Chem. Phys. 1970, 52, 5977–5987. (34) Hirschfelder, J. O. Almost Degenerate Perturbation Theory. Chem. Phys. Lett. 1978, 54, 1–3. (35) Hegarty, D.; Robb, M. A. Calculation of Effective Hamiltonians Using QuasiDegenerate Rayleigh-SchrÃűdinger Perturbation Theory (QD-RSPT). Mol. Phys. 1979, 37, 1455–1468. (36) Kirtman, B. Simultaneous Calculation of Several Interacting Electronic States by Generalized Van Vleck Perturbation Theory. J. Chem. Phys. 1981, 75, 798–808.

36

ACS Paragon Plus Environment

Page 36 of 61

Page 37 of 61 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

(37) Redmon, L. T.; Shavitt, I. Quasidegenerate Perturbation Theories. A Canonical Van Vleck Formalism and Its Relationship to Other Approaches. J. Chem. Phys. 1980, 73, 5711–5717. (38) Freed, K. F.; Sheppard, M. G. Convergence Studies of the Effective Valence Shell Hamiltonian for Correlation Energies of the Fluorine Atom and Its Ions Using Third Order Quasidegenerate Many-Body Perturbation Theory. J. Chem. Phys. 1981, 75, 4525–4538. (39) Freed, K. F.; Sheppard, M. G. Third-Order Quasidegenerate Many-Body Perturbation Theory Calculations for Valence State Correlation Energies of the Nitrogen and Oxygen Atoms and Their Ions. Int. J. Quantum Chem. 1981, 15, 21–31. (40) Bartlett, R. J.; Kucharski, S. A. Multireference Many-Body Perturbation Theory. Int. J. Quantum Chem. 1988, 22, 383–405. (41) Nakano, H. Quasidegenerate Perturbation Theory with Multiconfigurational SelfConsistent-Field Reference Functions. J. Chem. Phys. 1993, 99, 7983–7992. (42) Hirao, K.; Nakano, H.; Nakatani, J. Second-Order Quasi-Degenerate Perturbation Theory with Quasi-Complete Active Space Self-Consistent Field Reference Functions. J. Chem. Phys. 2001, 114, 1133–1141. (43) Shavitt, I.; Bartlett, R. J. Many-Body Methods in Chemistry and Physics, 1st ed.; Cambridge Press: New York, 2009; pp 185–250. (44) Roskop, L.; Gordon, M. S. Quasi-Degenerate Second-Order Perturbation Theory for Occupation Restricted Multiple Active Space Self-Consistent Field Reference Functions. J. Chem. Phys. 2011, 135, 044101. (45) Sharma, S.; Jeanmairet, G.; Alavi, A. Quasi-Degenerate Perturbation Theory Using Matrix Product States. J. Chem. Phys. 2016, 144, 034103. 37

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

(46) Andersson, K.; Malmqvist, P. A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. Second-Order Perturbation Theory with a CASSCF Reference Function. J. Phys. Chem. 1990, 94, 5484–5488. (47) Andersson, K.; Malmqvist, P. A.; Roos, B. O. Second-Order Perturbation Theory with a Complete Active Space Self-Consistent Field Reference Function. J. Phys. Chem. 1992, 96, 1218–1226. (48) Hirao, K. Multireference Møller–Plesset Method. Chem. Phys. Lett. 1992, 190, 374– 380. (49) Hirao, K. Multireference Møller–Plesset Perturbation Theory for High-Spin OpenShell Systems. Chem. Phys. Lett. 1992, 196, 397–403. (50) Hirao, K. Multireference Møller–Plesset Perturbation Treatment of Potential Energy Curve of N2 . Int. J. Quant. Chem. 1992, 26, 517–526. (51) Hirao, K. State-Specific Multireference Møller–Plesset Perturbation Treatment for Singlet and Triplet Excited States, Ionized States and Electron Attached States of H2 O. Chem. Phys. Lett. 1993, 201, 59–66. (52) Datta, B.; Mahapatra, U. S.; Mukherjee, D. Development of a Size-Consistent StateSpecific Multireference Perturbation Theory with Relaxed Model-Space Coefficients. Chem. Phys. Lett. 1999, 299, 42–50. (53) Datta, B.; Mahapatra, U. S.; Mukherjee, D. Molecular Applications of a SizeConsistent State-Specific Multireference Perturbation Theory with Relaxed ModelSpace Coefficients. J. Phys. Chem. A 1999, 103, 1822–1830. (54) Chattopadhyay, S.; Ghosh, P.; Jana, D.; Mukherjee, D. State-Specific Multi-Reference Perturbation Theories with Relaxed Coefficients: Molecular Applications. Int. J. Mol. Sci. 2002, 3, 733–754. 38

ACS Paragon Plus Environment

Page 38 of 61

Page 39 of 61 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

(55) Chattopadhyay, S.; Chaudhuri, R. K.; Mahapatra, U. S. Application of Improved Virtual Orbital Based Multireference Methods to N2 , LiF, and C4 H6 Systems. J. Chem. Phys. 2008, 129, 244108–244117. (56) Allen, W. D.; Evangelista, F. A.; Mukherjee, D.; Schaefer, H. F.; Simmonett, A. C. A Companion Perturbation Theory for State-Specific Multireference Coupled Cluster Methods. Phys. Chem. Chem. Phys. 2009, 11, 4728–4741. (57) 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–10264. (58) Angeli, C.; Borini, S.; Cestari, M.; Cimiraglia, R. A Quasidegenerate Formulation of the Second Order N-Electron Valence State Perturbation Theory Approach. J. Chem. Phys. 2004, 121, 4043–4049. (59) Angeli, C.; Borini, S.; Cavallini, A.; Cimiraglia, R. Third-Order Multireference Perturbation Theory: The N-Electron Valence State Perturbation-Theory Approach. J. Chem. Phys. 2006, 124, 054108–054116. (60) Hubac, I.; Mach, P.; Wilson, S. On the Application of Brillouin–Wigner Perturbation Theory to Multireference Configuration Mixing. Mol. Phys. 2002, 100, 859–863. (61) Hubac, I.; Mach, P.; Papp, P.; Wilson, S. Multireference Second-Order Brillouin– Wigner Perturbation Theory. Mol. Phys. 2004, 102, 701–709. (62) Hubac, I.; Mach, P.; Papp, P.; Pittner, J.; Wilson, S. Many-Body Brillouin–Wigner Second-Order Perturbation Theory Using a Multireference Formulation: An Application to Bond Breaking in the Diatomic Hydrides BH and FH. Mol. Phys. 2006, 104, 2367–2386.

39

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

(63) Hubac, I.; Mach, P.; Papp, P.; Wilson, S. Many-Body Brillouin–Wigner Second-Order Perturbation Theory: A Robust and Efficient Approach to the Multireference Correlation Problem. Int. J. Quantum Chem. 2007, 107, 2622–2631. (64) Hubac, I.; Mach, P.; Neogrady, P.; Papp, P.; Pittner, J.; Wilson, S. Many-Body Brillouin–Wigner Second-Order Perturbation Theory: An Application to the Autoaromatisation of Hex-3-Ene-1,5-Diyne (the Bergman Reaction). Mol. Phys. 2008, 106, 57–74. (65) Hoffmann, M. R. Canonical Van Vleck Quasidegenerate Perturbation Theory with Trigonometric Variables. J. Phys. Chem. 1996, 100, 6125–6130. (66) Hoffmann, M. R.; Khait, Y. G.; Song, J. Explication and Revision of Generalized Van Vleck Perturbation Theory for Molecular Electronic Structure. J. Chem. Phys. 2002, 117, 4133–4145. (67) Dudley, T. J.; Hoffmann, M. R.; Khait, Y. G. Molecular Gradients for the SecondOrder Generalized Van Vleck Variant of Multireference Perturbation Theory. J. Chem. Phys. 2003, 119, 651–660. (68) Hoffmann, M. R.; Jiang, W.; Khait, Y. G. Third-Order Generalized Van Vleck Perturbation Theory Variant of Multireference Perturbation Theory for Electron Correlation. Int. J. Quantum Chem. 2009, 109, 1855–1873. (69) Helgaker, T.; Hoffmann, M. R. Use of Density Functional Theory Orbitals in the GVVPT2 Variant of Second-Order Multistate Multireference Perturbation Theory. J. Phys. Chem. A 2015, 119, 1548–1553. (70) Whitten, J. L. Coulombic Potential Energy Integrals and Approximations. J. Chem. Phys. 1973, 58, 4496–4501.

40

ACS Paragon Plus Environment

Page 40 of 61

Page 41 of 61 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

(71) 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. (72) 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. (73) Vahtras, O.; Almlöf, J.; Feyereisen, M. W. Integral Approximations for LCAO-SCF Calculations. Chem. Phys. Lett. 1993, 213, 514–518. (74) 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. (75) 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. (76) Sodt, A.; Subotnik, J. E.; Head-Gordon, M. Linear Scaling Density Fitting. J. Chem. Phys. 2006, 125, 194109–194118. (77) Werner, H.-J.; Manby, F. R.; Knowles, P. J. Fast Linear Scaling Second-Order Møller– Plesset Perturbation Theory (MP2) Using Local and Density Fitting Approximations. J. Chem. Phys. 2003, 118, 8149–8160. (78) Schütz, M.; Manby, F. R. Linear Scaling Local Coupled Cluster Theory with Density Fitting. Part I: 4-External Integrals. Phys. Chem. Chem. Phys. 2003, 5, 3349–3358. (79) Werner, H.-J.; Schütz, M. An Efficient Local Coupled Cluster Method for Accurate Thermochemistry of Large Systems. J. Chem. Phys. 2011, 135, 144116–144131. (80) Beebe, N. H. F.; Linderberg, J. Simplifications in the Generation and Transformation

41

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

of Two-Electron Integrals in Molecular Calculations. Int. J. Quant. Chem. 1977, 12, 683–705. (81) Roeggen, I.; Wisloff-Nilssen, E. On the Beebe-Linderberg Two-Electron Integral Approximation. Chem. Phys. Lett. 1986, 132, 154–160. (82) Koch, H.; de Meras, A. S.; Pedersen, T. B. Reduced Scaling in Electronic Structure Calculations Using Cholesky Decompositions. J. Chem. Phys. 2003, 118, 9481–9484. (83) 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–194117. (84) DePrince, A. E.; Sherrill, C. D. Accuracy and Efficiency of Coupled-Cluster Theory Using Density Fitting/Cholesky Decomposition, Frozen Natural Orbitals, and a T1Transformed Hamiltonian. J. Chem. Theory Comput. 2013, 9, 2687–2696. (85) Bozkaya, U. Orbital-Optimized Second-Order Perturbation Theory with DensityFitting and Cholesky Decomposition Approximations: An Efficient Implementation. J. Chem. Theory Comput. 2014, 10, 2371–2378. (86) Bozkaya, U. Derivation of General Analytic Gradient Expressions for DensityFitted Post-Hartree-Fock Methods: An Efficient Implementation for the DensityFitted Second-Order Møller–Plesset Perturbation Theory. J. Chem. Phys. 2014, 141, 124108–124122. (87) Bozkaya, U. Analytic Energy Gradients and Spin Multiplicities for Orbital-Optimized Second-Order Perturbation Theory with Density-Fitting Approximation: An Efficient Implementation. J. Chem. Theory Comput. 2014, 10, 4389–4399. (88) Bozkaya, U. Orbital-Optimized MP3 and MP2.5 with Density-Fitting and Cholesky Decomposition Approximations. J. Chem. Theory Comput. 2016, 12, 1179–1188. 42

ACS Paragon Plus Environment

Page 42 of 61

Page 43 of 61 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

(89) Bozkaya, U. Orbital-Optimized Linearized Coupled-Cluster Doubles with DensityFitting and Cholesky Decomposition Approximations: An Efficient Implementation. Phys. Chem. Chem. Phys. 2016, 18, 11362–11373. (90) Bozkaya, U.; Sherrill, C. D. Analytic Energy Gradients for the Coupled-Cluster Singles and Doubles Method with the Density-Fitting Approximation. J. Chem. Phys. 2016, 144, 174103–174117. (91) Bozkaya, U. A Noniterative Asymmetric Triple Excitation Correction for the DensityFitted Coupled-Cluster Singles and Doubles Method: Preliminary Applications. J. Chem. Phys. 2016, 144, 144108–144120. (92) Bozkaya, U.; Sherrill, C. D. Analytic Energy Gradients for the Coupled-Cluster Singles and Doubles with Perturbative Triples Method with the Density-Fitting Approximation. J. Chem. Phys. 2017, 147, 044104–044114. (93) Bozkaya, U. Analytic Energy Gradients for Orbital-Optimized MP3 and MP2.5 with the Density-Fitting Approximation: An Efficient Implementation. J. Comp. Chem. 2018, 39, 351–360. (94) Parrish, R. M.; Burns, L. A.; Smith, D. G. A.; Simmonett, A. C.; DePrince, A. E.; Hohenstein, E. G.; Bozkaya, U.; Sokolov, A. Y.; Remigio, R. D.; Richard, R. M.; Gonthier, J. F.; James, A. M.; McAlexander, H. R.; Kumar, A.; Saitow, M.; Wang, X.; Pritchard, B. P.; Verma, P.; Schaefer, H. F.; Patkowski, K.; King, R. A.; Valeev, E. F.; Evangelista, F. A.; Turney, J. M.; Crawford, T. D.; Sherrill, C. D. Psi4 1.1: An OpenSource Electronic Structure Program Emphasizing Automation, Advanced Libraries, and Interoperability. J. Chem. Theory Comput. 2017, 13, 3185–3197. (95) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.;

43

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

Dupuis, M.; Montgomery, J. A. General Atomic and Molecular Electronic Structure System. J. Comp. Chem. 1993, 14, 1347–1363. (96) Gordon, M. S.; Schmidt, M. W. In Theory and Applications of Computational Chemistry: The First Forty Years, 1st ed.; Dykstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005; pp 1167–1189. (97) Rettrup, S. Alternative Graphical Representation of Determinantal Many-Electron States. Int. J. Quant. Chem. 2006, 106, 2511–2517. (98) Knowles, P. J.; Handy, N. C. A New Determinant-Based Full Configuration Interaction Method. Chem. Phys. Lett. 1984, 111, 315–321. (99) Olsen, J.; Roos, B. J.; Jorgensen, P.; Jensen, H. J. A. Determinant Based Configuration Interaction Algorithms for Complete and Restricted Configuration Interaction Spaces. J. Chem. Phys. 1988, 89, 2185–2192. (100) Davidson, E. R. The Iterative Calculation of a Few of the Lowest Eigenvalues and Corresponding Eigenvectors of Large Real-Symmetric Matrices. J. Comput. Phys. 1975, 17, 87–94. (101) Liu, B. In Numerical Algorithms in Chemistry: Algebraic Methods; Moler, C., Shavitt, I., Eds.; Technical Report LBL-8158, Lawrence Berkeley Laboratory, University of California: Berkeley, 1978; pp 49–53. (102) Leininger, M. L.; Sherrill, C. D.; Allen, W. D.; Schaefer, H. F. Systematic Study of Selected Diagonalization Methods for Configuration Interaction Matrices. J. Comp. Chem. 2001, 22, 1574–1589. (103) 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.

44

ACS Paragon Plus Environment

Page 44 of 61

Page 45 of 61 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

(104) Woon, D. E.; Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. V. Core-Valence Basis Sets for Boron Through Neon. J. Chem. Phys. 1995, 103, 4572–4585. (105) Peterson, K. A.; Dunning, T. H. Accurate Correlation Consistent Basis Sets for Molecular Core-Valence Correlation Effects: The Second Row Atoms Al-Ar, and the First Row Atoms B-Ne Revisited. J. Chem. Phys. 2002, 117, 10548–10560. (106) Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. (107) 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. (108) diatomic is a C++ program written for geometry optimization, frequency analysis, and spectroscopic constants computation for diatomic molecules via numerical derivatives by J. M. Turney, Center for Computational Quantum Chemistry, University of Georgia, Athens, GA, 30602, USA, (2010). (109) Irie, M.; Lifka, T.; Kobatake, S.; Kato, N. Photochromism of 1,2-Bis(2-methyl-5phenyl-3-thienyl)perfluorocyclopentene in a Single-Crystalline Phase. J. Am. Chem. Soc. 2000, 122, 4871–4876. (110) Uchida, K.; Matsuoka, T.; Kobatake, S.; Yamaguchi, T.; Irie, M. Substituent Effect on the Photochromic Reactivity of Bis(2-thienyl)perfluorocyclopentenes. Tetrahedron 2001, 57, 4559–4565. (111) Yanai, T.; Saitow, M.; Xiong, X.-G.; Chalupsky, J.; Kurashige, Y.; Guo, S.; Sharma, S. Multistate Complete-Active-Space Second-Order Perturbation Theory Based on Den-

45

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

sity Matrix Renormalization Group Reference States. J. Chem. Theory Comput. 2017, 13, 4829–4840. (112) H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, T. Korona, R. Lindh, A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, E. Goll, C. Hampel, A. Hesselmann, G. Hetzer, T. Hrenar, G. Jansen, C. Köppl, Y. Liu, A. W. Lloyd, R. A. Mata, A. J. May, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklass, D. P. O’Neill, P. Palmieri, D. Peng, K. Pflüger, R. Pitzer, M. Reiher, T. Shiozaki, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, M. Wang, MOLPRO, version 2012.1, a package of ab initio programs, 2012, see http://www.molpro.net. (113) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: A General-Purpose Quantum Chemistry Program Package. WIREs Comput. Mol. Sci. 2012, 2, 242–253. (114) Huber, K. P.; Herzberg, G. Molecular Spectra and Molecular Structure; Constants of Diatomic Molecules; Van Nostrand: Princeton, 1979; Vol. 4. (115) Csontos, B.; Nagy, B.; Csontos, J.; Kállay, M. Dissociation of the Fluorine Molecule. J. Phys. Chem. A 2013, 117, 5518–5528. (116) Larsen, H.; Olsen, J.; Jørgensen, P.; Christiansen, O. Full Configuration Interaction Benchmarking of Coupled-Cluster Models for the Lowest Singlet Energy Surfaces of N2 . J. Chem. Phys. 2000, 113, 6677–6686. (117) Dutta, A.; Sherrill, C. D. Full Configuration Interaction Potential Energy Curves for Breaking Bonds to Hydrogen: An Assessment of Single-Reference Correlation Methods. J. Chem. Phys. 2003, 118, 1610–1619.

46

ACS Paragon Plus Environment

Page 46 of 61

Page 47 of 61 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

Table 1: Total energies for the ground state, ET (in hartrees), equilibrium bond lengths, re (in Å), and singlet-triplet energy splittings (EST , in kcal mol-1 ) between the states X 3 Σg and a 1 ∆g for the O2 molecule from the DF-QDPT2, CD-QDPT2 (with a CD tolerance of 10−6 ), CASSCF, MRMP2, MC-QDPT2, and CASPT2 methods along with the cc-pwCVTZ basis set. For each method a CMS consisting of two orbitals, one α-electron, and one β-electron is used [CMS(211)]. Method DF-QDPT2 CD-QDPT2 DF-CAS-QDPT2a CASSCF MRMP2 MC-QDPT2 CASPT2 Expb a

ET −150.211 −150.211 −150.237 −149.658 −150.179 −150.278 −150.265

166 284 811 912 573 923 933

re (S) re (T) 1.271 1.254 1.271 1.254 1.329 1.253 1.153 1.152 1.389 1.315 1.240 1.297 1.351 1.178 1.216 1.208

DF-CASSCF MOs were used for DF-QDPT2. b Huber and Herzberg. 114

47

ACS Paragon Plus Environment

EST 22.8 22.7 25.3 29.3 21.6 25.2 14.6 22.6

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

Page 48 of 61

Table 2: Total energies for the ground state, ES (in hartrees), equilibrium bond lengths, re (in Å), and singlet-triplet energy splittings (EST , in kcal mol-1 ), between the states X 1 Σg and a 3 Πu for the C2 molecule from the DF-QDPT2, CD-QDPT2 (with a CD tolerance of 10−6 ), CASSCF, MRMP2, MC-QDPT2, and CASPT2 methods along with the aug-ccpwCVTZ basis set. For each method, a CMS consisting of three orbitals, two α-electrons, and two β-electrons is used [CMS(322)]. Method DF-QDPT2 CD-QDPT2 DF-CAS-QDPT2a CASSCF MRMP2 MC-QDPT2 CASPT2 Expb a

ES −75.790 −75.790 −75.814 −75.493 −75.772 −75.860 −75.799

129 139 767 405 554 456 472

re (S) re (T) 1.259 1.284 1.259 1.284 1.250 1.311 1.197 1.171 1.271 1.220 1.269 1.221 1.386 1.368 1.243 1.312

DF-CASSCF MOs were used for DF-QDPT2. b Huber and Herzberg. 114

48

ACS Paragon Plus Environment

EST −3.4 −3.4 −11.9 −13.6 −29.6 −28.5 16.9 −2.0

Page 49 of 61 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

Table 3: Total energies, E(F2 ) (in hartrees), equilibrium bond lengths, re (in Å), and dissociation energies, De (in kcal mol-1 ) for the F2 molecule from the DF-QDPT2, CD-QDPT2 (with a CD tolerance of 10−10 ), CASSCF, MRMP2, MC-QDPT2, and CASPT2 methods along with the cc-pwCVTZ basis set. For each method, a CMS consisting of six orbitals, five α-electrons, and five β-electrons is used [CMS(655)]. Method DF-QDPT2 CD-QDPT2 DF-CAS-QDPT2a CASSCF MRMP2 MC-QDPT2 CASPT2 Expv a

E(F2 ) −199.345 644 −199.345 730 −199.356 543 −198.830 252 −199.295 172 −199.402 361 −199.400 172

re 1.423 1.423 1.420 1.463 1.423 1.422 1.421 1.412

De 36.4 36.4 38.6 17.8 36.6 36.5 36.7 38.3

DF-CASSCF MOs were used for DF-QDPT2. b Huber and Herzberg. 114

49

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

Page 50 of 61

Table 4: Total energies (in hartrees) from the FCI,a DF-QDPT2, CASSCF, CASPT2, and MRCI+Q methods for the N2 molecule with the cc-pVDZ basis set. For each multireference method a CMS consisting of 4 orbitals, 2 α-electrons, and 2 β-electrons is used [CMS(422)]. RN N (Å) FCIa 0.7938 −108.624 700 0.9525 −109.167 573 1.0679 −109.270 384 1.0943 −109.276 528 1.1113 −109.278 135 1.1208 −109.278 339 1.1473 −109.276 583 1.1737 −109.271 915 1.2700 −109.238 397 1.4288 −109.160 305 a

DF-QDPT2 −108.599 885 −109.141 989 −109.251 776 −109.259 146 −109.261 669 −109.262 353 −109.261 794 −109.258 437 −109.228 329 −109.159 345

CASSCF −108.386 889 −108.895 328 −108.970 742 −108.970 190 −108.967 413 −108.965 131 −108.956 320 −108.944 462 −108.883 186 −108.759 874

Larsen, Olsen, Jørgensen, and Christiansen 116

50

ACS Paragon Plus Environment

CASPT2 −108.613 762 −109.156 849 −109.262 364 −109.269 504 −109.271 861 −109.272 512 −109.272 140 −109.269 076 −109.243 464 −109.143 887

MRCI+Q −108.628 943 −109.169 285 −109.270 102 −109.275 689 −109.276 933 −109.276 922 −109.274 536 −109.269 212 −109.232 879 −109.160 661

Page 51 of 61 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

Table 5: Total energies (in hartrees) from the FCI,a DF-QDPT2, CASSCF, CASPT2, and MRCI+Q methods for the CH4 molecule with the 6-31G(d) basis set. For each multireference method a CMS consisting of 4 orbitals, 3 α-electrons, and 3 β-electrons is used [CMS(433)]. RCH (Å) 0.8 0.9 1.0 1.1 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 a

FCIa −40.253 342 −40.320 513 −40.349 369 −40.356 202 −40.350 579 −40.322 605 −40.289 114 −40.258 549 −40.233 555 −40.214 618 −40.201 257 −40.192 439 −40.186 932 −40.183 629

DF-QDPT2 −40.235 153 −40.302 407 −40.330 453 −40.337 430 −40.330 373 −40.299 594 −40.266 475 −40.236 318 −40.211 075 −40.192 577 −40.179 472 −40.171 000 −40.167 150 −40.165 521

CASSCF −40.106 498 −40.172 622 −40.200 252 −40.205 659 −40.206 713 −40.181 080 −40.149 652 −40.120 644 −40.096 710 −40.078 473 −40.065 593 −40.057 126 −40.051 882 −40.048 770

Dutta and Sherrill 117

51

ACS Paragon Plus Environment

CASPT2 −40.237 681 −40.304 428 −40.332 779 −40.338 971 −40.334 146 −40.306 344 −40.272 995 −40.242 591 −40.217 820 −40.199 137 −40.185 979 −40.177 261 −40.171 766 −40.168 431

MRCI+Q −40.258 003 −40.325 009 −40.353 744 −40.360 484 −40.354 717 −40.326 695 −40.293 185 −40.262 615 −40.237 626 −40.218 702 −40.205 358 −40.196 557 −40.191 065 −40.187 773

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

Page 52 of 61

Table 6: Total energies (in hartrees) from the DF-QDPT2, CASSCF, CASPT2, and MRCI+Q methods for the F2 molecule with the cc-pwCVTZ basis set. For each method a CMS consisting of 6 orbitals, 5 α-electrons, and 5 β-electrons is used [CMS(655)]. RF F (Å) 1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00 3.20 3.60 4.00 5.00 6.00 7.00 8.00

DF-QDPT2 −199.086 337 −199.296 660 −199.345 338 −199.334 637 −199.314 372 −199.299 205 −199.290 587 −199.287 003 −199.286 734 −199.287 394 −199.287 969 −199.287 761 −199.288 263 −199.287 587 −199.288 343 −199.287 998 −199.288 432 −199.287 898

CASSCF −198.548 876 −198.780 944 −198.828 665 −198.826 045 −198.815 622 −198.808 296 −198.804 513 −198.802 776 −198.802 046 −198.801 783 −198.801 718 −198.801 723 −198.801 770 −198.801 813 −198.801 876 −198.801 897 −198.801 904 −198.801 907

52

CASPT2 −199.125 869 −199.356 499 −199.399 939 −199.389 675 −199.371 304 −199.358 055 −199.350 715 −199.347 061 −199.345 333 −199.344 553 −199.344 217 −199.344 075 −199.343 993 −199.343 995 −199.344 035 −199.344 052 −199.344 058 −199.344 061

ACS Paragon Plus Environment

MRCI+Q −199.139 898 −199.371 420 −199.414 405 −199.403 223 −199.384 512 −199.371 583 −199.364 718 −199.361 421 −199.359 908 −199.359 246 −199.358 972 −199.358 864 −199.358 810 −199.358 820 −199.358 864 −199.358 881 −199.358 887 −199.358 890

Page 53 of 61 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

Table 7: Singlet-triplet energy splittings (EST , in kcal mol-1 ) between the ground state S0 and the lowest triplet state T1 for the diarylethene derivatives from the DF-QDPT2 method along with the 6-311G(2d,2p) basis set. For each molecule, there are 1053 basis functions in the 6-311G(2d,2p) basis set, and a CMS consisting of 2 orbitals, 1 α-electrons, and 1 β-electrons is used [CMS(211)]. Molecule N-diarylethene (open) N-diarylethene (closed) I-diarylethene (open) I-diarylethene (closed)

53

EST 65.0 24.6 54.1 38.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

Virtual (a00 , b00 , . . . ) Particles (a, b, . . . ) Valence particles (a0 , b0 , . . . )

Valence holes (i0 , j 0 , . . . ) Holes (i, j, . . . ) Core (i00 , j 00 , . . . )

Figure 1: Classification of the molecular orbital space in the Hose–Kaldor approach.

54

ACS Paragon Plus Environment

Page 54 of 61

Page 55 of 61 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

Figure 2: Chemical structures of the diradicals considered.

55

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

800 700

CAS-QDPT2 CASSCF DF-QDPT2

600 Wall Time (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

500 400 300 200 100 0 1

2

3

4

5

6

7

8

9 10 11 12 13

Set Entry

Figure 3: Total wall time (in min) for computations of single-point energies for the diradical set considered from the DF-QDPT2, CASSCF (from Psi4), and CAS-QDPT2 methods with the cc-pwCVTZ basis set. For the largest member, the number of basis functions is 501. All computations were performed on a single node (1 core) Intel(R) Xeon(R) CPU E5-2620 v3 @ 2.40 GHz computer (memory ∼64 GB).

56

ACS Paragon Plus Environment

Page 56 of 61

Page 57 of 61

FCI CASSCF CASPT2 MRCI+Q DF-QDPT2

-108.4000

-108.6000 Energy (hartree)

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

-108.8000

-109.0000

-109.2000

-109.4000 0.7

0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

R (Angstrom)

Figure 4: Potential energy curves for N2 in a cc-pVDZ basis set using various approximate correlation methods. For each multireference method, a CMS consisting of four orbitals, two α-electrons, and two β-electrons is used [CMS(422)].

57

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

-40.0000 -40.0500

FCI CASSCF CASPT2 MRCI+Q DF-QDPT2

-40.1000

Energy (hartree)

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 58 of 61

-40.1500 -40.2000 -40.2500 -40.3000 -40.3500 -40.4000 0.5

1.0

1.5

2.0

2.5

3.0

R (Angstrom)

Figure 5: Potential energy curves for CH4 in a 6-31G(d) basis set using various approximate correlation methods. For each multireference method, a CMS consisting of four orbitals, three α-electrons, and three β-electrons is used [CMS(433)].

58

ACS Paragon Plus Environment

Page 59 of 61

-198.5000

CASSCF CASPT2 MRCI+Q DF-QDPT2

-198.6000 -198.7000 -198.8000 Energy (hartree)

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

-198.9000 -199.0000 -199.1000 -199.2000 -199.3000 -199.4000 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 R (Angstrom)

Figure 6: Potential energy curves for F2 in a cc-pwCVTZ basis set using various approximate correlation methods. For each multireference method, a CMS consisting of six orbitals, five α-electrons, and five β-electrons is used [CMS(655)].

59

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

I-diarylethene (closed)

I-diarylethene (open)

N-diarylethene (closed)

N-diarylethene (open)

Figure 7: Chemical structures of the diarylethene derivatives considered.

60

ACS Paragon Plus Environment

Page 60 of 61

Page 61 of 61 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

TOC Graphic

DF-QDPT2/6-311G(2d,2p) 1053 basis functions 5.1 hours

61

ACS Paragon Plus Environment