Large-Scale Variational Two-Electron Reduced ... - ACS Publications

Subscriber access provided by McMaster University Library

Article

Large-scale v2RDM-driven CASSCF methods Jacob Fosso-Tande, Truong-Son Nguyen, Gergely Gidofalvi, and A. Eugene DePrince J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00190 • Publication Date (Web): 11 Apr 2016 Downloaded from http://pubs.acs.org on April 12, 2016

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

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

Page 1 of 39

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

Journal of Chemical Theory and Computation

Large-scale v2RDM-driven CASSCF methods Jacob Fosso-Tande,† Truong-Son Nguyen,‡ Gergely Gidofalvi,‡ and A. Eugene DePrince III∗,† Department of Chemistry and Biochemistry, Florida State University, Tallahassee, FL 32306-4390, and Department of Chemistry and Biochemistry, Gonzaga University, Spokane, Washington 99258-0089 E-mail: [email protected]

Abstract A large-scale implementation of the complete active space self-consistent field (CASSCF) method is presented. The active space is described using the variational two-electron reduced-density-matrix (v2RDM) approach, and the algorithm is applicable to much larger active spaces than can be treated using configuration-interaction-driven methods. Density fitting or Cholesky decomposition approximations to the electron repulsion integral tensor allow for the simultaneous optimization of large numbers of external orbitals. We have tested the implementation by evaluating singlet-triplet energy gaps in the linear polyacene series and two dinitrene biradical compounds. For the acene series, we report computations that involve active spaces consisting of as many as 50 electrons in 50 orbitals and the simultaneous optimization of 1892 orbitals. For the dinitrene compounds, we find that the singlet-triplet gaps obtained from v2RDMdriven CASSCF with partial three-electron N -representability conditions agree with those obtained from configuration-interaction-driven approaches to within one third of ∗

To whom correspondence should be addressed Department of Chemistry and Biochemistry, Florida State University, Tallahassee, FL 32306-4390 ‡ Department of Chemistry and Biochemistry, Gonzaga University, Spokane, Washington 99258-0089 †

1

ACS Paragon Plus Environment


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

one kcal mol−1 . When enforcing only the two-electron N -representability conditions, v2RDM-driven CASSCF yields less accurate singlet-triplet energy gaps in these systems, but the quality of the results is still far superior to those obtained from standard single-reference approaches.

1

Introduction

The standard approach for capturing nondynamical correlation effects in quantum chemical computations is the complete active space self-consistent field (CASSCF) method. 1–4 In the traditional approach to CASSCF, the electronic structure of the active space is represented by a full configuration interaction (CI) expansion of the wave function, the complexity of which increases exponentially with the size of the active space. Hence, the application of CIdriven CASSCF to large active spaces is nontrivial. In order to describe large active spaces efficiently, one must either (i) abandon full CI in favor of some other wave function expansion that scales polynomially with system size or (ii) abandon the wave function altogether. The density-matrix renormalization group (DMRG) approach 5–14 offers an alternative expansion of the wave function that facilitates the design of polynomially-scaling CASSCF. 15–18 Indeed, DMRG-driven CASSCF methods represent the current state of the art for describing static electron correlation in systems with large active spaces. Alternatively, methods that employ the two-electron reduced-density matrix (2-RDM) as the central variable, such as the variational 2-RDM (v2RDM) method, 19–34 also have the potential to overcome the scaling limitations that plague CI. In principle, v2RDM-driven CASSCF should be competitive with the DMRG-based approach, and v2RDM methods offer the potential advantage over DMRG that a v2RDM-based optimization should be less sensitive to the physical dimensionality of the system than DMRG. Further, v2RDM-driven CASSCF is slightly more “black-box” in the sense that one need not localize orbitals or seek an optimal orbital ordering, as is recommended for DMRG computations. 35 Additional polynomially scaling 2-RDM-based methods can be achieved through the solution of the contracted Schrödiger equation, 36–38 the antiher2


Page 2 of 39

Page 3 of 39

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


mitian contracted Schrödinger equation, 39–41 or the particle-hole hypervirial equations. 42,43 Despite recent theoretical 44–49 and algorithmic 29,50–55 advances, the chemistry community has resisted adopting v2RDM-based approaches to the static correlation problem. One possible explanation for this reluctance is a practical one: few implementations of v2RDMdriven CASSCF and related methods have been reported, 49,56–59 and, to our knowledge, only one code, which performs a v2RDM-driven optimization corresponding to a seniority-zero CI wave function, 49 is publically available. The present work represents our efforts to remedy this situation. We have developed a free and open-source implementation of a v2RDM-driven CASSCF method as a plugin to the Psi4 electronic structure package. 60 The active-space 2-RDM is obtained from a semidefinite optimization procedure, without the use of the N electron wave function. The active-space 2-RDM computation is generalized to treat both open- and closed-shell systems and to enforce several well-known conditions to ensure that the 2-RDM represents a physically meaningful wave function. 19,29,61 Density fitting 62–65 or Cholesky decomposition 66–69 approximations to the electron repulsion integral (ERI) tensor facilitate the optimization of large numbers of external orbitals. This paper is organized as follows. Section 2.1 provides an overview of the v2RDM approach and conditions for the feasibility of the 2-RDM, Section 2.2 describes the semidefinite optimization procedure used to obtain the active-space 2-RDM, and Section 2.3 summarizes the main features of the orbital optimization procedure. A quasi one-step v2RDM-CASSCF procedure is discussed in Section 2.4, and in Section 3 we explore the numerical performance of the present approach for several challenging systems that require the consideration of strong correlation effects. For the linear polyacene series, we find that singlet-triplet energy gaps from CASSCF computations using 2-RDMs that satisfy two-electron N -representability conditions 19 are in good agreement with those from DMRG. 70 The largest member of the series considered consists of 12 fused benzene rings and requires an active space of 50 active orbitals with 50 electrons as well as the simultaneous optimization of 1892 orbitals. Additionally, we explore the singlet-triplet energy gap in a family of dinitrene biradical compounds

3



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

and find that excellent agreement with CI-driven CASSCF can be obtained by imposing partial three-electron N -representability conditions. 29,61

2

Theory

It has long been recognized that the ground-state electronic energy for a many-electron system is an exact functional of the 2-RDM. 71–73 However, the variational determination of the ground-state 2-RDM without knowledge of the wave function is complicated by the fact that not all 2-RDMs correspond to realistic N -electron wave functions. A 2-RDM that does correspond to an N -electron wave function is said to be N -representable, 74 and substantial effort has been dedicated to the practical application of known N -representability conditions. 19–24 If the N -representability conditions are exactly satisfied, then a variational 2-RDM (v2RDM) computation is numerically equivalent to the full CI; that is, it is exact within a given basis set. Similarly, if one performs an exact v2RDM computation within an active space, the result is equivalent to a complete active space (CAS) CI computation. By simultaneously minimizing the electronic energy with respect to rotations between inactive, active, and external orbitals, one achieves a v2RDM-driven CASSCF procedure. Below, we provide the details of such an approach. Throughout, we follow the usual conventions for orbital index labels in CASSCF: the indices i, j, k, and l represent inactive (doubly occupied) orbitals; the indices t, u, v, w, x, and y represent active (partially occupied) orbitals; the indices a, b, c, and d correspond to external (empty) orbitals; and the indices p, q, r, and s represent general orbitals.

2.1

N -representability conditions

The active-space electronic Hamiltonian in second-quantized notation can be expressed as X X X′ † X ′ † 1X † ˆ H= a ˆtσ a ûτ , (1) htu + [2(tu|ii) − (ti|ui)] (tv|uw) a ˆtσ a ûτ a ˆwλ a ˆvκ + 2 tuvw στ tu i στ κλ 4


Page 4 of 39

Page 5 of 39

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


where (tv|uw) represents an ERI in Mulliken notation, htu represents the sum of electron kinetic energy and electron-nucleus potential energy integrals, and the symbols a ˆ† and a ˆ represent second-quantized creation and annihilation operators, respectively. The Greek P′ labels run over α and β spin, and the symbol indicates that the only nonvanishing contributions to the sums are those for which the total number of creation and annihilation operators for a given spin are equal. The CAS electronic energy is then given by X X X′ X ′ 1X † † hΨ|ˆ a†tσ a ûτ |Ψi. E= htu + [2(tu|ii)−(ti|ui)] (tv|uw) hΨ|ˆ atσ a ûτ a ˆwλ a ˆvκ |Ψi+ 2 tuvw στ tu i στ κλ (2) where the expectation values define elements of the 2-RDM

2

ˆ†uτ a ˆwλ a ˆvκ |Ψi, a†tσ a Dvtσκuwτλ = hΨ|ˆ

(3)

ûτ |Ψi. a†tσ a Dutστ = hΨ|ˆ

(4)

and the 1-RDM 1

Hence, the CAS electronic energy is a known functional of the active-space 1- and 2-RDM. If our goal is to determine of the elements of the 1- and 2-RDM without knowledge of the wave function, we must establish some necessary conditions for the N -representability of the RDMs. Given a ground-state N -electron wave function, |Ψi, one can generate a set of basis functions {|ΨI i} through the application of a set of operators {CÎ } as |ΨI i = CÎ |Ψi.

(5)

The Gram matrix formed from the inner product of any set of these functions is a positive semidefinite matrix, and N -representability conditions can be organized into a hierarchy of

5



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

Page 6 of 39

conditions according the number of creation/annihilation operators that comprise {CÎ }. 26 If {CÎ } contains a single annihilation operator, we obtain the definition of the 1-RDM in Eq. (4), and if we take {CÎ } to be comprised of one creation operator, we obtain the definition of the 1-hole RDM 1

atσ a ˆ†uτ |Ψi. Qtuστ = hΨ|ˆ

(6)

Because the Gram matrix is positive semidefinite, 1 D ≥ 0 and 1 Q ≥ 0, where the notation M ≥ 0 denotes a positive semidefinite matrix M with non-negative eigenvalues. Elements of the 1-RDM and the 1-hole RDM are related to one another by linear mappings that are implied by the anticommutation relations of fermionic creation/annihilation operators. At the two-body level, if we limit {CÎ } to products of two creation/annihilation operators, we obtain a set of two-body N -representability conditions known as the P, Q, and G (or PQG) conditions. 19 First, consider the set of operators that generates (N -2)-electron wave functions: CÎ ∈ {ˆ awλ a ˆvκ }. The elements of the Gram matrix for this vector space are the elements of the 2-RDM, given by Eq. (3); hence, an N -representable 2-RDM must be positive ˆ†vκ }, we find that the two-hole semidefinite (the P-condition). Next, by choosing CÎ ∈ {ˆ a†wλ a RDM, whose elements are defined by

2

ˆ†vκ |Ψi, atσ a ûτ a ˆ†wλ a Qtvσκuwτλ = hΨ|ˆ

(7)

ˆvκ }, we should be positive semidefinite (the Q-condition). Lastly, by selecting CÎ ∈ {ˆ a†wλ a find that the particle-hole RDM, the elements of which are

2

ˆvκ |Ψi, a†tσ a ûτ a ˆ†wλ a Gtvσκuwτλ = hΨ|ˆ

(8)

should be positive semidefinite (the G-condition). Each of these matrices are related to one another according to the linear mappings implied by the anticommutation relations of fermionic creation/annihilation operators.

6


Page 7 of 39

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


At the three-body level, four unique reduced-density matrices can be defined ˆ†uτ a ˆ†vκ a ˆyν a ˆxµ a ˆwλ |Ψi, a†tσ a Dwtσλuxτµvyκν = hΨ|ˆ

(9)

a†tσ a ˆ†uτ a ˆvκ a ˆ†yν a ˆxµ a ˆwλ |Ψi, Ewtσλuxτµvyκν = hΨ|ˆ

(10)

3

ûτ a ˆvκ a ˆ†yν a ˆ†xµ a ˆwλ |Ψi, Fwtσλuxτµvyκν = hΨ|ˆ a†tσ a

(11)

3

atσ a ûτ a ˆvκ a ˆ†yν a ˆ†xµ a ˆ†wλ |Ψi, Qtwσλuxτµvyκν = hΨ|ˆ

(12)

3

3

each of which should be positive semidefinite and related to each other and the 1- and 2electron RDMs according to the anticommutation relations for fermionic creation/annihilation operators. Because the sum of any two positive semidefinite matrices is also positive semidefinite, the positivity of the RDMs in Eqs. (9) - (12) implies a weaker set of partial threeelectron conditions. The T1 and T2 conditions 29,61 maintain the positivity of

T1 = 3 D + 3 Q,

(13)

T2 = 3 E + 3 F.

(14)

and

The right hand sides of Eqs. (13) and (14) depend on only the 1- and 2-RDMs. Our present semidefinite solver, outlined in Sec. 2.2, can enforce the PQG and T1/T2 conditions for both closed- and open-shell systems. We note that the strength of the T2 condition presented here is not invariant to the order of the creation and annihiliation operators contained in CÎ for the E and F conditions. A strengthened T2 condition (sometimes denoted T2′ 75 ) that does display this property has also been reported. 61,75,76 Our present software does not yet enforce the T2′ condition. Additional constraints may be derived by considering the number of electrons in the active space and the total spin of the system. For non-relativistic Hamiltonians, the wave function has well-defined total spin (S) and projection of spin (MS ). As a result, the only 7



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

Page 8 of 39

nonzero reduced-density matrix elements are those for which the total number of creation and annihilation operators for a given spin are equal. Thus, the 1-RDM is comprised of the 1

Dαα block with elements 1

ûα |Ψi, Dutαα = hΨ|ˆ a†tα a

1

ûβ |Ψi, Duββ = hΨ|ˆ a†tβ a

(15)

and the 1 Dββ block with elements t

(16)

and the 1 Dαβ off-diagonal block is zero. Similarly, the 2-RDM is expressed in terms of 2 ββ 2 αβ uncoupled 2 Dαα αα , Dββ , and Dαβ blocks. The three-particle matrix T1 consists of four nonβββ αββ ααβ redundant uncoupled spin blocks (T1ααα ααα , T1ααβ , T1αββ , and T1βββ ), while T2 has a much

more complicated spin structure: 

ααα T2ααα



T2αββ ααα

0 0 0 0     T2ααα T2αββ 0 0 0 0   αββ αββ     βββ βαα 0 0   0 0 T2βαα T2βαα . T2 =    βββ βαα   0 0 0 T2 0 T2 βββ βββ     ααβ  0 0  0 0 0 T2ααβ     ββα 0 0 0 0 0 T2ββα

(17)

Expressions for the spin blocks of the remaining RDMs are given elsewhere. 33,77 A spin-pure active-space 2-RDM must satisfy the expression 77,78 X tu

2

1 t u Duαα tββ = (Nα + Nβ ) + MS2 − S(S + 1), 2

(18)

where Nα and Nβ denote the number of active electrons with α and β spin, respectively. It has been shown previously that the extremal spin projection states are numerically the most well-constrained; 79 hence, the present work treats all open-shell systems as the maximal spin

8


Page 9 of 39

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


projection. Because the Hamiltonian conserves the number of electrons, several trace and contraction conditions must also be obeyed by N -representable RDMs. For example, the traces of each spin block of the active-space 2-RDM are X

2

X

2

Dttααuuαα = Nα (Nα − 1),

(19)

tu

t u

Dtββ uββ = Nβ (Nβ − 1),

(20)

tu

and X

2

t u

Dtααuββ = Nα Nβ .

(21)

tu

Further, spin blocks of the 2-RDM are also related to spin blocks of the 1-RDM according to the contractions X

2

X

2

=

X

2

=

X

2

(Nα − 1)1 Dutαα =

Dutααvvαα ,

(22)

Duββ vββ ,

t v

(23)

Duαα vββ ,

t v

(24)

v t

(25)

v

t

(Nβ − 1)1 Duββ =

v

Nβ Dutαα 1

v

t Nα Duββ 1

Dvααuββ .

v

Lastly, the block structure of each of the RDMs can be refined beyond the spin considerations discussed above to include the spatial symmetry of the one-, two-, or three-particle basis functions. We fully exploit Abelian point group symmetry, which can dramatically reduce the number of non-zero elements in each RDM and, consequently, the number of N -representability constraints that must be explicitly enforced. We can further reduce the number of non-zero elements of each of the RDMs by recognizing that an N -representable RDM must be antisymmetric with respect to particle exchange (e.g. 2 Dvtααuwαα = −2 Dvuααwtαα = −2 Dwtααuvαα = 2 Dwuααtvαα ). We enforce this trivial constraint through the use of antisymmetrized 9



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

Page 10 of 39

two-particle and three-particle basis functions. For example, the 2 Dαα αα block of the 2-RDM is spanned by n2A antisymmetric functions rather than n2A symmetric functions. Similarly, ααβ nA and n2A nA basis functions, respecthe T1ααα ααα and T1ααβ blocks of T1 are spanned by 3 tively. Here, nA denotes the number orbitals in the active space. The hermiticity condition

(2 Dvtσκuwτλ = 2 Dtvσκuwτλ ) is also enforced, but not in a way that reduces the number of elements of the RDM that must be stored.

2.2

Semidefinite optimization algorithm

The variational determination of the ground-state active-space 2-RDM is a semidefinite optimization problem, 24–34 the primal formulation of which is Eprimal = cT · x

minimize such that and

Ax

=b

M (x)

≥ 0.

(26)

Here, x represents the primal solution vector, the vector c contains all information defining the quantum system (the 1- and 2-electron integrals), and the mapping M (x) maps the primal solution onto the set of positive semidefinite RDMs 

1

0

D

   0   M (x) =   0    0  0

1

0

0

0



  Q 0 0 0    0 2D 0 0   ≥ 0.   2 0 0 Q 0   2 0 0 0 G

(27)

When enforcing the partial three-electron N -representability conditions, M (x) contains additional blocks corresponding to T1 and T2. The action of the constraint matrix, A, on x is a compact representation of the N -representability conditions. For example, three rows of 10


Page 11 of 39

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


A correspond to the three trace conditions given by Eqs. (19)-(21). Upon convergence, the numerical result of the action of those rows of A on x will equal the right-hand sides of Eqs. (19)-(21), which are each elements of the constraint vector b. Thus, A maintains the appropriate mappings between each block of M (x) implied by the fermionic creation/annihilation operator anticommutation relations and enforces the appropriate spin, trace, and contraction conditions. Alternatively, one could consider the dual formulation of the semidefinite problem, expressed as

maximize such that and

Edual = bT · y z

(28)

= c − AT y

M (z) ≥ 0

where y and z are the dual solutions, and M (z) is constrained to be positive semidefinite. To compute the primal and dual solutions, we use a boundary-point semidefinite optimization algorithm 53–55 that reduces to the following iterative two-step procedure: 1. Solve AAT y = A(c − z) + τ µ(b − Ax) for y by conjugate gradient methods. 2. Diagonalize Ω = M (µx + AT y − c), and separate the positive and negative eigenvalues and corresponding eigenvectors, P+ and P− . Update x and z according to M (x) = µ−1 P+ Ω+ PT+ and M (z) = −P− Ω− PT− , where Ω+/− are diagonal matrices with the positive and negative eigenvalues of Ω on the diagonal. Here, τ is a step-length parameter that lies in the interval [1.0,1.6]; 55 we choose τ = 1.6 in all of our computations. The penalty parameter µ controls how strictly the primal or dual constraints are enforced and is updated every 500 macroiterations, according to the protocol outlined in Ref. 55. The algorithm is considered converged when the primal error ||Ax − b||, the dual error ||AT y − c + z||, and the primal/dual energy gap |Eprimal − Edual | are sufficiently small. We note that the convergence properties of the boundary-point approach 11



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

Page 12 of 39

to the v2RDM problem differ substantially from many familiar algorithms, such as the Davidson diagonalization procedure in CI-driven CASSCF. Each macroiteration typically involves 50-100 conjugate gradient steps, and the overall algorithm may involve thousands or even tens of thousands of macroiterations to reach convergence. The computational cost of the semidefinite optimization of the active-space 2-RDM scales polynomially with the size of the active space. The first step requires O(n4A ) or O(n6A ) floating-point operations when enforcing the PQG or T1/T2 conditions, respectively, where nA represents the number of active orbitals. The second step requires O(n6A ) or O(n9A ) floating-point operations when enforcing the PQG or T1/T2 conditions. The storage requirements for the two-electron and three-electron matrices are O(n4A ) and O(n6A ), respectively. These favorable floating-point costs and storage requirements allow for the consideration of much larger active spaces than can be treated with CI-based algorithms.

2.3

Orbital optimization

Here we provide the specific details of our orbital optimization procedure. For more details on the optimization of the orbitals in CASSCF, the reader is referred to Refs. 4 and 80–82. The quasi-Newton orbital optimization is based on the exponential parameterization of the orbital transformation 83,84 U = eK ,

(29)

where U is a unitary matrix and the updated orbitals (C′ ) are obtained from the expansion of the current orbitals (C) in the molecular orbital basis according to C′ = CU. The unique elements of the skew-symmetric matrix K may be collected in the vector κ, where the element κpq (p < q) corresponds to the orbital mixing coefficient between orbitals p and q. The CASSCF energy is invariant to unitary transformations among inactive or external orbitals, and the active-space v2RDM energy is invariant to rotations among active orbitals. Thus, only the non-redundant inactive-active, inactive-external, and active-external orbital

12


Page 13 of 39

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


rotations need to be considered. The expression for the CASSCF energy truncated at secondorder in κ is 1 E(κ) = E(0) + κT g + κT Bκ, 2

(30)

where g and B are the orbital gradient and Hessian, respectively. The variational condition, g = 0, is solved iteratively using the approximate Newton step ˜ −1 gk , κk+1 = κk − B k

(31)

˜ k is the diagonal, approximate Hessian matrix at iteration k. 85 Thus, the orbital where B optimization requires the repeated evaluation of the energy, gradient, and diagonal Hessian, computation of the matrix exponential in Eq. (29), and transformation of the integrals. Expressions for the gradient and approximate diagonal Hessian matrix elements have been presented elsewhere, 85 and here we only summarize the expressions employed in our implementation. The gradient elements for non-redundant orbital rotations, expressed in ¯ and Z ¯ terms of the inactive (I F) and active (A F) Fock matrices as well as the auxiliary Y matrices, are Fti + A Fti − 2 Y¯ti + Z¯ti , = 4 I Fai + A Fai , = 2 Y¯ta + Z¯ta .

git = 4 gia gta

I

(32) (33) (34)

The computational cost for evaluating the preconditioner in Eq. (31) is considerably reduced by the use of approximate expressions ˜it,it = 4 B ˜ia,ia B ˜ta,ta B

I

Ftt + A Ftt − I Fii − A Fii + 2γtt = 4 I Faa + A Faa − I Fii − A Fii , = 2γtt I Faa + A Faa − 2 Y¯tt + Z¯tt , 13

I

Fii + A Fii − 2 Y¯tt + Z¯tt ,


(35) (36) (37)


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

Page 14 of 39

for the diagonal elements of the Hessian matrix. The spatial 1-RDM is computed from the t

spin blocks of the 1-RDM as γtu = 1 Dutαα +1 Duββ . Provided the auxiliary matrices are available, ˜ scales as the number of non-redundant the computational cost for evaluating both g and B orbital rotation pairs. Thus, the overall cost of updating κ in Eq. (31) is determined by the ¯ and Z. ¯ effort for evaluating the auxiliary matrices I F, A F, Y, Although our implementation allows the use of conventional four-index integrals, the use of three-index ERIs

N aux X

(pq|rs) ≈

bPpq bPrs ,

(38)

P

significantly reduces the I/O cost for large-scale applications. The factors bPpq depend on the type of integral decomposition employed, and Naux represents the number of auxiliary basis functions or Cholesky vectors in the density fitting 62–65 or Cholesky decomposition 66–69 approximations, respectively. The matrix elements for the inactive Fock matrix are given by I

Fpq = hpq + 2

X

bPpq

X i

P

bPii

−

X Pi

bPpi bPqi

,

(39)

and the active Fock matrix elements are contractions of the spatial 1-RDM with the ERIs A

Fpq =

X P

bPpq

X

γtu bPtu

tu

X 1X P P − bpt bqu . γtu 2 tu P

(40)

γtu I Fpu ,

(41)

The remaining auxiliary matrices are given by Z¯tp =

X u

and Y¯tp =

X uP

bPpu

X vw

Γtuvw bPvw

,

(42)

where the symmetrized spatial 2-RDM (Γ) displays the same eightfold permutational symmetry as the ERIs. Parentheses denote the order of operations we have chosen for the

14


Page 15 of 39

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


evaluation of intermediate quantities. For a sequence of calculations with basis sets of increasing size and the same active space, as would be encountered when the CASSCF energy is extrapolated to the infinite basis set limit, the transformation of the three-index integral factors, bPpq , dominates the cost of the orbital optimization procedure. For each auxiliary function P , the transformation of bPpq for all p and q requires two n × n matrix multiplications. Therefore, the integral transformation step requires O(n3 Naux ) effort, where n is the total number of orbitals. The next rate-limiting step is the evaluation of the matrix exponential in Eq. (29) after each κ update. To avoid complex arithmetic, the matrix exponential is computed using trigonometric functions; 83,84 the necessary operations involve a dense diagonalization, as well as matrix-matrix multiplications with a computational cost that scales as O(n3 ). The remaining cost of the orbital optimization is in the construction of the Fock matrix and other auxiliary matrices. Since the Fock matrix elements for which p 6= q and both p and q are external orbitals are not needed, the number of matrix elements that need to be computed scales as O(nnocc ). Here, nocc = nI + nA is the total number of occupied orbitals, and nI and nA are the number of inactive and active orbitals, respectively. Using three-index integrals, the number of ¯ floating-point operations required for the rate limiting steps in evaluating I F, A F, and Y are O(nocc nI nNaux ), O(nocc n2A nNaux ), and O(n4A Naux + n2A nNaux ), respectively. Computing ¯ requires negligible O(n2 n) computational effort. Z A

2.4

Overall v2RDM CASSCF procedure

Prior implementations of the v2RDM-CASSCF method fully converged the semidefinite optimization problem before each orbital optimization step. 56 For large active spaces, the active-space v2RDM optimization should dominate the cost of the overall computation, and thus the total computation time of a truly two-step v2RDM-CASSCF algorithm is expected to increase linearly with the number of orbital optimization (OO) steps. A onestep procedure could reduce the overall number of the active-space v2RDM optimization 15



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

steps by updating the orbitals at each step of the v2RDM procedure. Such an approach may be impractical, however. First, it is not immediately clear how important the coupling would be between the OO and CAS-v2RDM steps. Second, because the v2RDM procedure typically involves thousands or sometimes even tens of thousands of iterations, the number of orbital optimization steps would drastically increase. A good compromise between these approaches is a quasi one-step algorithm where the orbital optimization step outlined above is carried out once every M v2RDM iterations. This approach has the desirable feature that it retains the simplicity of the two-step approach, and, as we demonstrate in Sec. 3, it leads to significant computational savings. Figure 1 illustrates the overall algorithm for the quasi one-step v2RDM-CASSCF procedure. Overall convergence is achieved when both the v2RDM optimization and OO satisfy their respective convergence criteria. Unless otherwise noted, the v2RDM optimization is considered converged when both the primal and dual errors, ||Ax − b|| and ||AT y − c + z||, fall below 1 × 10−5 and the primal-dual energy gap, |Eprimal − Edual |, falls below 1 × 10−4 Eh . The OO procedure is considered converged when the change in the energy due to orbital rotations is less than 1 × 10−8 Eh and the magnitude of the orbital gradient falls below 1 × 10−4 Eh . Figure 1: Overall v2RDM-CASSCF procedure.

16


Page 16 of 39

Page 17 of 39

3

Results and Discussion

To facilitate v2RDM-CASSCF computations with large numbers of external orbitals, we employ DF/CD approximations to the ERI tensor. Here, we establish that the errors due to the use of the DF/CD approximations are negligible relative to the errors stemming from the incompleteness of the PQG, T1, and T2 N -representability conditions. Figure 2(a) illustrates the potential energy curve for the dissociation of N2 in a cc-pVQZ basis set; we present results for both CI- and v2RDM-driven CASSCF approaches. The calculations employ the D2h point group and a (6,6) active space, where the notation (ne , no ) denotes an active space with ne electrons and no orbitals. The 1ag , 2ag , 1b1u , and 2b1u orbitals are doubly occupied, and the 3ag , 3b1u , 1b2u , 1b3u , 1b2g , and 1b3g orbitals comprise the active space. Within the full D∞h molecular point group, these active orbitals correspond to the σg , σu , πu , and πg orbitals, respectively. The qualitative features of the CI-CASSCF curve are reproduced well Figure 2: The (a) potential energy curves for the dissociation of molecular nitrogen using CIand v2RDM-driven CASSCF approaches, and the (b) error in the v2RDM-CASSCF energy, relative to CI-CASSCF. All CASSCF methods use the cc-pVQZ basis set and a (6,6) active space. −108.6

0.000 CI PQG PQG+T1T2

−108.7

−0.005 error (Eh)

−108.8 energy (Eh)

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


−108.9

−0.010

PQG PQG+T1 PQG+T2 PQG+T1T2

−109.0 −0.015 −109.1

−109.2

(b)

(a) 1.0

1.5

2.0 2.5 N−N bond length (Å)

3.0

−0.020

3.5

1.0

1.5


3.0

3.5

by the v2RDM-CASSCF approach using either the PQG or PQG+T1T2 conditions, but the PQG conditions alone fail to provide quantitative agreement with CI-driven CASSCF over

17



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

the entire potential energy curve [Fig. 2(b)]. This discrepancy highlights the importance of the partial three-electron conditions, especially in the vicinity of the spin-recoupling region. Near equilibrium (at RN −N = 1.1 ˚ A), the energy error in the PQG results (-4.5 mEh ) are comparable to the errors of PQG+T1T2 conditions (-0.9 mEh ). The largest-magnitude energy error using the PQG conditions is almost -20 mEh (at RN −N = 1.7 ˚ A), while the T1 and T2 conditions result in a fourfold reduction in the error (-5 mEh ). Also included in Fig. 2(b) are curves that illustrate the error in the v2RDM-CASSCF energy using the PQG+T1 and PQG+T2 conditions. The T1 condition improves upon the PQG potential energy curves but does not improve the quality of the v2RDM-CASSCF energy with PQG+T2 conditions. Hence, in practice, it may be reasonable to ignore this condition when enforcing the much stronger T2 condition. Figure 3 illustrates the error introduced by the use of the DF or CD approximations in a v2RDM-CASSCF computation using the PQG+T1T2 N -representability conditions. The DF fitting basis was the cc-pVQZ-JKFIT basis, and the tolerance for the Cholesky decomposition was 1×10−4 . The DF approximation introduces errors between 4 and 63 µEh , while the largest error introduced by the CD approximation is less than 3 µEh . Neither of these errors are significant relative to the magnitude of the errors resulting from incomplete N -representability conditions. For the remainder of this work, all v2RDM-CASSCF computations employ the DF approximation, and the fitting basis is the JKFIT basis set that corresponds to the primary basis set. Figure 4 illustrates the occupations of the highest occupied natural orbital (HONO), lowest unoccupied natural orbital (LUNO), the HONO-1, and the LUNO+1 for N2 at the v2RDM-CASSCF(6,6)/cc-pVQZ level of theory. Here, the HONO refers to natural orbital number nA /2 when the natural orbitals are ordered in terms of decreasing occupancy. Similarly, the LUNO represents natural orbital number nA /2 + 1. Near equilibrium, the HONO (πu ) and HONO-1 (σg ) are doubly occupied, and the LUNO (πg ) and LUNO+1 (σu ) orbitals are unoccupied. All six frontier orbitals become degenerate at the dissociation limit, and, as would be expected based on the relative strength of σ and π interactions, the onset of

18


Page 18 of 39

Page 19 of 39

Figure 3: Error in the potential energy curve for the dissociation of N2 introduced by the DF/CD approximations. Results are presented for the v2RDM-CASSCF method in a ccpVQZ basis set, using a (6,6) active space and enforcing the PQG+T1T2 N -representability conditions. 70.0 DF CD

60.0

error (µEh)

50.0 40.0 30.0 20.0 10.0 0.0

1.0

1.5


3.0

3.5

degeneracy for the πu /πg orbitals occurs at slightly shorter bond lengths than for the σg /σu orbitals. Figure 4: Natural orbital occupation numbers convey the onset of degeneracy of the πu /πg orbitals and the σg /σu orbitals for the dissociation of the N2 molecule. 2.0

HONO (πu) LUNO (πg) HONO−1 (σg)

1.5 occupation

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


LUNO+1 (σu)

1.0

0.5

0.0 1.0

1.5


3.0

3.5

To demonstrate the applicability of the v2RDM-CASSCF approach to the description 19



Figure 5: Singlet-triplet energy gaps as a function of the number of fused benzene rings (k) in linear polycyclic aromatic hydrocarbons. Results are presented at the v2RDM-CASSCF level of theory using the cc-pVDZ and cc-pVTZ basis sets and enforcing the PQG N representability conditions. The v2RDM-CASSCF gaps agree well with those obtained from CAS-DMRG 70 and experiment. 86–89 70 singlet-triplet gap (kcal mol-1)

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 39

PQG (cc-pVDZ) PQG (cc-pVTZ) DMRG (DZ) experiment

60 50 40 30 20 10 0 2

3

4

5

6 7 8 k-acene

9

10

11

12

of large active spaces and the optimization of large numbers of external orbitals, we next consider the linear polyacene series. The ground states for the linear polyacene molecules are known to display complex singlet polyradical behavior, 70 and the series has become a popular test case for demonstrating the ability of an electronic structure method to capture strong electron correlation effects. 33,56,57,70,90–92 A proper description of the electronic structure of the larger members of this series requires much larger active spaces than can be used within conventional CI-driven CASSCF approaches. Figure 5 illustrates the singlet-triplet energy gap for the linear acene series computed at the v2RDM-CASSCF level of theory using the PQG two-electron N -representability conditions. The D2h molecular geometries for the singlet states were optimized using the Becke three-parameter Lee-Yang Par (B3LYP) functional and the 6-31G* basis set, as implemented in the Psi4 electronic structure package. Geometries for the triplet states were optimized at the unrestricted B3LYP (UB3LYP)/631G* level of theory. These optimized geometries, as well as absolute v2RDM-CASSCF energies are available in the Supporting Information. For a linear polyacene with k fused 20


Page 21 of 39

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


six-membered rings, v2RDM-CASSCF calculations employ a (4k + 2, 4k + 2) active space (2k + 1 π bonding orbitals and 2k + 1 antibonding π ∗ orbitals). The active orbitals are constructed from pz atomic orbitals, and there are k orbitals of b2g symmetry, k orbitals of au symmetry, k + 1 orbitals of b3g symmetry and k + 1 orbitals of b1u symmetry. The largest system considered, 12-acene, requires an active space with 50 electrons in 50 orbitals. In the cc-pVTZ basis set, the total number of orbitals (inactive + active + external) for 12-acene is 1892; all of these orbitals are optimized in the v2RDM-CASSCF procedure. With the exception of 11-acene, the singlet-triplet energy gap decreases monotonically with increasing k; we suspect that this unusual behavior for 11-acene stems from the mismatch in levels of theory used for geometry optimization and the evaluation of the singlet-triplet energy gaps. Also included in Fig. 5 are DMRG and experimental results taken from Refs. 70 and 86–89. The DMRG computations in Ref. 70 differ slightly from the present computations in that they employ a smaller Dunning-Hay double-ζ basis set 93,94 and a larger active space comprised of the same number of π orbitals used in this study, but twice as many π ∗ antibonding orbitals. Additionally, we note that the DMRG computations in Ref. 70 do not include orbital rotations and thus correspond to a DMRG-driven CASCI computation, rather than DMRG-driven CASSCF. Nonetheless, the present v2RDM-CASSCF singlet-triplet gaps agree quite well with those obtained from DMRG. The singlet-triplet gaps at the v2RDMCASSCF level of theory remain essentially unchanged from the cc-pVDZ to the cc-pVTZ basis set; the largest difference, observed for naphthalene, is 0.7 kcal mol−1 . Complete-basisset limit gaps (available in the Supporting Information) obtained using Helgaker’s two-point extrapolation procedure 95 are slightly larger than either the cc-pVDZ or cc-pVTZ results, a trend that agrees with that reported in Ref. 96 using CI-driven CASSCF methods and significantly more limited active spaces. The v2RDM-CASSCF/cc-pVTZ natural orbitals presented in Fig. 6(a) clearly convey the onset of di- and even polyradical behavior in the singlet ground state of the linear polyacene series for large k. Similar trends were reported previously using both CAS-DMRG 70

21



and CAS-v2RDM 33,56 approaches. The present results indicate that the occupation numFigure 6: The (a) occupation numbers for all of the natural orbitals and the (b) HONO and LUNO for the singlet ground state of linear acene molecules. (a) 2.0

(b) 2.0

2−acene 3−acene 4−acene 5−acene 6−acene 7−acene 8−acene 9−acene 10−acene 11−acene 12−acene

1.6 1.4 1.2 1.0

b1u

1.8

occupation

1.8

occupation

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

Page 22 of 39

0.8

b3g

1.4

au

1.2 1.0 0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.0

b2g

1.6

0.0 0

5

10

15 20 25 30 35 natural orbital number

40

45

50

1

2

3

4

5

6 7 k−acene

8

9

10

11

12

bers for the HONO and LUNO become nearly equal at 10-acene but do not remain so for larger polyacenes. Instead, the occupations of the au /b3g and b1u /b2g orbitals cross at 11-acene [see Fig. 6(b)]. This behavior was recently reported using the thermally-assisted occupation density-functional theory (TAO-DFT) method 97 but has not been observed using CAS-DMRG or CAS-v2RDM approaches. Previous v2RDM-CASCI 33 and DMRG-CASCI 70 computations in minimal basis sets suggested that no such crossing occurs, even at 12-acene. Our results thus indicate that the larger basis sets and orbital optimization are necessary for correctly predicting the qualitative behavior of the natural orbital occupation numbers in these molecules. For the linear polyacene series, the PQG N -representability conditions appear to be sufficient to achieve good agreement between v2RDM- and DMRG-driven CASSCF methods. In other cases, however, one must consider the application of stronger conditions, such as the T1/T2 conditions, for some assurance of quantitative accuracy. Here, we demonstrate that v2RDM-CASSCF with partial three-electron conditions can achieve sub kcal mol−1 accuracy in the singlet-triplet energy gaps in a family of dinitrene biradical molecules, relative to CI22


Page 23 of 39

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


driven CASSCF. Table 1 provides the singlet-triplet energy gaps for 1,4-phenylenedinitrene and biphenyl-4,4’-dinitrene. The CASSCF computations employ (10,10) and (16,16) active spaces for 1,4-phenylenedinitrene and biphenyl-4,4’-dinitrene, respectively. For a dinitrene compound with k rings, the active space includes the 3k + 1 occupied π molecular orbitals, the 3k + 1 π ∗ antibonding molecular orbitals, as well as the singly-occupied σ orbitals on the nitrogen atoms that lie in the molecular plane. More details on the active spaces are provided in the Supporting Information. The geometries for the dinitrene molecules were optimized using the NWChem package 98 at the CASSCF/cc-pVDZ level of theory, using a (10,10) active space for 1,4-phenylenedinitrene and an (8,8) active space for biphenyl-4,4’-dinitrene. Both the singlet and triplet energies were computed using the triplet geometry, and the singlettriplet gaps reported herein thus correspond to vertical gaps. Results for several standard single-reference methods, including second-order Møller-Plesset perturbation theory (MP2), coupled-cluster with single and double excitations (CCSD), 99 CCSD with perturbative triple excitations [CCSD(T)], 100 and B3LYP are also included for comparison. Additionally, we report the singlet-triplet gap computed using equations-of-motion spin-flip CCSD (EOM-SFCCSD) approach, 101 which was designed for the balanced description of singlet and triplet biradicals, such as the dinitrene compounds in this study. Energies from the single-reference methods were obtained using the Q-Chem electronic structure package. 102 Electron paramagnetic resonance (EPR) measurements indicate that the singlet-triplet gaps in 1,4-phenylenedinitrene and biphenyl-4,4’-dinitrene are in the range 0.572–0.722 and 0.578–0.750 kcal mol−1 , respectively. 104,105 More recent magneto-optical measurements suggest that the singlet-triplet gap in 1,4-phenylenedinitrene may be as small as 0.179 kcal mol−1 . 103 These small gaps are well-reproduced by methods that account for strong electron correlation effects; for 1,4-phenylenedinitrene and biphenyl-4,4’-dinitrene, CI-driven CASSCF predicts that the singlet states lie 0.919 and 0.395 kcal mol−1 below the high-spin states, respectively. The v2RDM-CASSCF method with PQG+T1T2 N -representability conditions yields similar predictions; v2RDM-CASSCF gaps differ from those obtained by

23



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

Page 24 of 39

Table 1: Absolute energies (Eh ) and relative energies (kcal mol−1 ) for the lowest-energy singlet and triplet states for 1,4-phenylenedinitrene and biphenyl-4,4’-dinitrene.a,b molecule

1,4-phenylenedinitrene

singlet Expt. (Ref. 103) – Expt. (Ref. 104) – Expt. (Ref. 105) – CASSCF −338.4826 PQG+T1T2 −338.4829 PQG −338.5024 EOM-SF-CCSD −339.4285 MP2 −339.0179 CCSD −339.0641 CCSD(T) −339.3457 B3LYP −340.3389 a

triplet gap – 0.179 – 0.572 – 0.722 −338.4811 0.919 −338.4810 1.181 −338.4951 4.528 −339.2693 99.896 −339.3534 −210.527 −339.4331 −231.523 −339.4806 −84.679 −340.4191 −50.379

biphenyl-4,4’-dinitrene

singlet – – – −568.0936 −568.0949 −568.1348 −569.6272 −569.5176 −569.5894 −569.6849 −571.3806

triplet gap – – – 0.578 – 0.750 −568.0930 0.395 −568.0938 c 0.705 −568.1286 3.848 −569.5283 62.041 −569.5293 −7.344 −569.6383 −30.694 −569.7491 −40.338 −571.4609 −50.409

All v2RDM-CASSCF energies reported correspond to the primal formulation of the energy. b All v2RDM-CASSCF computations employ the density fitting approximation. c The primal/dual energy gap was converged to 0.2 mE . h

CI-driven methods by less than one third of one kcal mol−1 for both molecules. When enforcing the PQG conditions alone, v2RDM-CASSCF overstabilizes the singlet state, resulting in larger gaps of 4.528 and 3.848 kcal mol−1 for 1,4-phenylenedinitrene and biphenyl-4,4’dinitrene, respectively. At first, these larger errors seem to indicate a failure of the v2RDMdriven approach, but the quality of the results is still far superior to those obtained from standard single-reference approaches. Neither MP2, CCSD, CCSD(T), nor B3LYP are capable of correctly predicting the most stable spin state. Additionally, the magnitude of the spin gap varies wildly depending on the level of theory employed, and only B3LYP correctly predicts that the magnitude of the singlet-triplet gaps should be similar for the two compounds. The only single-reference method to predict the correct sign for the spin gaps is the EOM-SF-CCSD approach, but the magnitude of these gaps is two orders of magnitude too large. Further, for the B3LYP functional, we obtain similar results with similar trends using triplet geometries optimized at the B3LYP level of theory; hence, our conclusions here are 24


Page 25 of 39

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


independent of the quality of the geometry. We note that the absolute energy for the triplet state of 1,4-phenylenedinitrene obtained from v2RDM-CASSCF using three-particle N -representability conditions is slightly higher than that obtained from CI-driven CASSCF. However, the energy obtained from a v2RDMdriven approach that employs incomplete N -representability conditions should always be a lower-bound to the exact result (in this case that obtained from CI-driven CASSCF). The density fitting approximation is the source of this discrepancy; the v2RDM-CASSCF computations employ DF, while the CI-CASSCF computations use conventional 4-index integrals. We also note that, as was observed for the dissociation of N2 in Fig. 2(b), the T1 condition does not significantly improve upon the T2 condition. When enforcing the PQG+T2 conditions, the v2RDM-CASSCF approach predicts that the singlet-triplet energy gap for 1,4-phenylenedinitrene is 1.142 kcal/mol, which is within 0.04 kcal mol−1 of the PQG+T1T2 result. Lastly, we evaluate the effect of the frequency of the orbital optimization (OO) procedure on the overall time-to-solution and the distribution of effort between the OO procedure and the semidefinite optimization that constitutes the CAS portion of a v2RDM-CASSCF calculation. Figure 7 illustrates the total wall time for the two- and quasi one-step v2RDMCASSCF procedures for the singlet states of the linear acene series in the cc-pVDZ and cc-pVTZ basis sets. All computations enforce the PQG N -representability conditions, and all timings are evaluated using a single core of a 2.4 GHz Intel Xeon E5-2630 CPU. In the quasi one-step procedure, the orbitals are updated at intervals of either every 250 or 500 v2RDM optimization steps. These frequencies are somewhat arbitrary, and we are currently working to establish reliable protocols for dynamically determining the update interval. The present v2RDM-CASSCF computations are seeded with self-consistent field (SCF) orbitals and the corresponding uncorrelated 1- and 2-RDMs. In the early portion of the semidefinite optimization, the 1- and 2-RDMs may stray far from these “reasonable” values while the boundary-point optimization improves the quality of the dual solution (which is seeded with

25



zeros, in general). As a check on the quality of the 1- and 2-RDM, we require that the v2RDM energy fall below the SCF energy before initiating the orbital optimization procedure. In both the cc-pVDZ [Fig. 7(a)] and cc-pVTZ [Fig. 7(b)] basis sets, the quasi one-step procedure clearly outperforms the two-step procedure for all system sizes. An illustration of the cost of the quasi one-step procedure relative to the two-step procedure is provided in the insets of Fig. 7. The advantage of the one-step procedure is slightly more apparent in the smaller basis set, where the increase in the number of orbital optimization steps does not significantly impact the overall time-to-solution. In the cc-pVDZ basis, the average cost of the quasi onestep procedure is 47% or 57% that of the two step procedure, when using orbital update frequencies of 250 or 500, respectively. The increase in the percentages (to 53% and 61%, respectively) upon increasing the basis set to cc-pVTZ is not significant. An orbital update frequency of 250 results in shorter computation times for almost all systems considered, but these differences are small. Figure 7: Total wall time (hours) for two- and quasi one-step v2RDM-CASSCF for k-acene in the (a) cc-pVDZ and (b) cc-pVTZ basis sets. The ratio of the total time required for the quasi one-step procedures to the total time required for the two-step algorithm is provided in the insets. 45.0

45.0 0.8 0.7

30.0

0.6

35.0

0.5

30.0

time (hours)

ratio

35.0

0.4

25.0

0.3

20.0

2 3 4 5 6 7 8

15.0

1−step−250 1−step−500 2−step

10.0 5.0

0.8

40.0

0.7 ratio

40.0

time (hours)

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 26 of 39

0.6 0.5 0.4

25.0

0.3

20.0

2 3 4 5 6 7 8

15.0 10.0 5.0

(a)

0.0

(b)

0.0 2

3

4

5 k−acene

6

7

8

2

3

4

5 k−acene

6

7

8

Figure 8 illustrates the percentage of the total time-to-solution spent in the OO and CAS procedures for the linear polyacene series. In a cc-pVDZ basis set [Fig. 8(a)], the CAS 26


Page 27 of 39

portion of the calculation dominates the cost of the computation for all system sizes. For active spaces larger than 14 electrons in 14 orbitals (anthracene and beyond), the orbital optimization constitutes less than 10% of the total wall time, regardless of the orbital update frequency. In general, the amount of time spent in the OO procedure reflects the chosen orbital update frequency; the percentage is smallest for the two-step procedure and largest for the quasi one-step procedure with and update frequency of 250. In the cc-pVTZ basis set [Fig. 8(b)], the OO procedure constitutes a larger portion of the total time-to-solution, and the percentage of time spent in the OO procedure generally decreases with system size. The v2RDM-CASSCF optimization for 8-acene involves a (34,34) active space and the simultaneous optimization of 1300 orbitals; for this calculation, the OO step constitutes between 11 and 15% of the total time, depending on the update frequency. Figure 8: Percentage of total wall time spent within the CAS (solid lines) and OO (dashed lines) portions of the v2RDM-CASSCF procedure. Results are presented for the linear acene series in the (a) cc-pVDZ and (b) cc-pVTZ basis sets using 2-step (black lines) and quasi one-step procedures where the orbitals are updated every 250 (red lines) or 500 (green lines) iterations. 100.0

100.0

90.0

90.0

80.0

50.0

70.0 % total time

60.0

(b)

80.0

CAS(250) CAS(500) CAS OO(250) OO(500) OO

70.0 % total time

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


40.0 30.0

60.0 50.0 40.0 30.0

20.0

20.0

(a)

10.0

10.0

0.0

0.0 2

3

4

5 k−acene

6

7

8

2

27

3


4

5 k−acene

6

7

8


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

4

Conclusions

We have presented a v2RDM-driven implementation of the CASSCF method. The code is a free and open-source plugin to the Psi4 electronic structure that can be downloaded from GitHub (https://github.com/edeprince3/v2rdm casscf). The present implementation improves upon previously reported v2RDM-driven CASSCF implementations in three ways: (i) integral transformations involving DF/CD integrals are less computationally intensive than those involving conventional four-index integrals; (ii) the present quasi-Newton orbital optimization procedure is less computationally expensive than the Jacobi procedure utilized in Ref. 56; and (iii) ignoring the cost of the orbital optimiztion procedure itself, the quasi one-step procedure outlined in Section 2.4 yields an overall algorithm that would be roughly twice as efficient as the two-step procedure used in Ref. 56. Using two-particle N -representability conditions, we explored the singlet-triplet energy gaps in the linear acene series and demonstrated that v2RDM-CASSCF provides gaps in good agreement with both experimentally obtained values and those from DMRG. The favorable scaling of the v2RDM method with two-particle conditions allows for the consideration of much larger active spaces than can be treated using more conventional CI-driven methods; the largest computation in the polyacene series, 12-acene, involved 50 active electrons in 50 orbitals. For the same system, the use of the DF approximation enabled the simultaneous optimization of 1892 orbitals. We also explored singlet-triplet energy gaps in a family of dinitrene biradicals. CI-driven CASSCF accurately reproduces the small energy gaps inferred from EPR and magneto-optical spectroscopy experiments. Using partial three-particle N representability conditions, v2RDM-driven CASSCF energy gaps agree with those from CICASSCF to within one third of one kcal mol−1 . The v2RDM approach with only two-particle N -representability conditions correctly predicts that the singlet is the more stable spin state, but the energy gaps are too large by 3-4 kcal mol−1 . The only single-reference method able to predict the correct sign of the spin gap was the EOM-SF-CCSD method, but the EOMSF-CCSD energy gaps are two orders of magnitude too large. 28


Page 28 of 39

Page 29 of 39

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


While the formal scaling properties of the v2RDM-CASSCF method are quite appealing, the convergence of the boundary-point approach to the semidefinite problem is less than optimal, especially when enforcing the three-particle N -representability conditions. In this case, it is not unreasonable to expect the procedure to take many thousands or tens of thousands of iterations to reach convergence. We stress that these convergence issues only become problematic when considering the application of three-particle N -representability conditions. Hence, future work will consider strategies to improve the convergence of the boundary-point semidefinite programming solver. Additional future work will consider methods to capture dynamical electron correlation within the framework of v2RDM-driven CASSCF; several good candidates have already been proposed. Chan and coworkers have developed a canonical transformation theory 106–108 to incorporate dynamic electron correlation into DMRG-CASSCF. 17 Yanai and coworkers have developed internally contracted second-order perturbation theory 109 and configuration interaction 110 methods based on DMRG references. These computations express the dynamical correlation energy in terms of the 2-, 3-, and 4-RDM. More recently, Li and Evangelista presented a multireference version of the driven similarity renormalization group (DSRG) approach 111 to calculate a perturbative estimate of the dynamical correlation energy. 112 This last method is quite appealing in the context of v2RDM-driven CASSCF, as the energy expression depends only on the 3-RDM, which can be constrained explicitly in the semidefinite optimization or approximated using a cumulant decomposition. 37,113–115 Acknowledgments. Acknowledgment is made to the donors of the American Chemical Society Petroleum Research Fund for partial support of this research (Grant number 54668-DNI6). AED acknowledges support from the Ralph E. Powe Junior Faculty Enhancement Award from Oak Ridge Associated Universities. Supporting Information. 29



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

B3LYP-optimized geometries and v2RDM-CASSCF energies for the linear polyacene molecules can be found in the Supporting Information. Also provided are the CASSCFoptimized geometries and active-space information for the dinitrene biradical compounds. This information is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Roos, B. O.; Taylor, P. R.; Siegbahn, P. E. Chemical Physics 1980, 48, 157 – 173. (2) Siegbahn, P.; Heiberg, A.; Roos, B.; Levy, B. Phys. Scripta 1980, 21, 323–327. (3) Siegbahn, P. E. M.; Almlöf, J.; Heiberg, A.; Roos, B. O. J. Chem. Phys. 1981, 74, 2384–2396. (4) Roos, B. O. Advances in Chemical Physics; John Wiley & Sons, Inc., 1987; Vol. 68; pp 399–445. (5) White, S. R. Phys. Rev. Lett. 1992, 69, 2863–2866. (6) Wilson, K. G. Rev. Mod. Phys. 1975, 47, 773–840. (7) Schollwöck, U. Rev. Mod. Phys. 2005, 77, 259–315. (8) Yan, S.; Huse, D. A.; White, S. R. Science 2011, 332, 1173–1176. (9) Mitrushenkov, A. O.; Linguerri, R.; Palmieri, P.; Fano, G. J. Chem. Phys. 2003, 119, 4148–4158. (10) Chan, G. K.-L.; Head-Gordon, M. J. Chem. Phys. 2002, 116, 4462–4476. (11) Legeza, O.; Fáth, G. Phys. Rev. B 1996, 53, 14349–14358. (12) Marti, K. H.; Ond´ık, I. M.; Moritz, G.; Reiher, M. J. Chem. Phys. 2008, 128, 014104. (13) Zgid, D.; Nooijen, M. J. Chem. Phys. 2008, 128, 014107. 30


Page 30 of 39

Page 31 of 39

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


(14) Kurashige, Y.; Yanai, T. J. Chem. Phys. 2009, 130, 234114. (15) Ghosh, D.; Hachmann, J.; Yanai, T.; Chan, G. K.-L. J. Chem. Phys. 2008, 128, 144117. (16) Yanai, T.; Kurashige, Y.; Ghosh, D.; Chan, G. K.-L. Int. J. Quantum Chem. 2009, 109, 2178–2190. (17) Yanai, T.; Kurashige, Y.; Neuscamman, E.; Chan, G. K.-L. J. Chem. Phys. 2010, 132, 024105. (18) Wouters, S.; Poelmans, W.; Ayers, P. W.; Neck, D. V. Comput. Phys. Commun. 2014, 185, 1501 – 1514. (19) Garrod, C.; Percus, J. K. J. Math. Phys. 1964, 5, 1756–1776. (20) Garrod, C.; Mihailović, M. V.; Rosina, M. J. Math. Phys. 1975, 16, 868–874. (21) Mihailović, M.; Rosina, M. Nucl. Phys. A 1975, 237, 221 – 228. (22) Rosina, M.; Garrod, C. J. Comput. Phys. 1975, 18, 300 – 310. (23) Erdahl, R. M.; Garrod, C.; Golli, B.; Rosina, M. J. Math. Phys. 1979, 20, 1366–1374. (24) Erdahl, R. Rep. Math. Phys. 1979, 15, 147 – 162. (25) Nakata, M.; Nakatsuji, H.; Ehara, M.; Fukuda, M.; Nakata, K.; Fujisawa, K. J. Chem. Phys. 2001, 114, 8282–8292. (26) Mazziotti, D. A.; Erdahl, R. M. Phys. Rev. A 2001, 63, 042113. (27) Mazziotti, D. A. Phys. Rev. A 2002, 65, 062511. (28) Mazziotti, D. A. Phys. Rev. A 2006, 74, 032501. (29) Zhao, Z.; Braams, B. J.; Fukuda, M.; Overton, M. L.; Percus, J. K. J. Chem. Phys. 2004, 120, 2095–2104. 31



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

(30) Fukuda, M.; Braams, B. J.; Nakata, M.; Overton, M. L.; Percus, J. K.; Yamashita, M.; Zhao, Z. Math. Prog. 2007, 109, 553. (31) Cancès, E.; Stoltz, G.; Lewin, M. J. Chem. Phys. 2006, 125, 064101. (32) Verstichel, B.; van Aggelen, H.; Van Neck, D.; Ayers, P. W.; Bultinck, P. Phys. Rev. A 2009, 80, 032508. (33) Fosso-Tande, J.; Nascimento, D. R.; DePrince III, A. E. Mol. Phys. 2015, 0, 1–8. (34) Verstichel, B.; van Aggelen, H.; Neck, D. V.; Bultinck, P.; Baerdemacker, S. D. Comput. Phys. Commun. 2011, 182, 1235 – 1244. (35) Olivares-Amaya, R.; Hu, W.; Nakatani, N.; Sharma, S.; Yang, J.; Chan, G. K.-L. J. Chem. Phys. 2015, 142, 034102. (36) Colmenero, F.; Valdemoro, C. Phys. Rev. A 1993, 47, 979–985. (37) Nakatsuji, H.; Yasuda, K. Phys. Rev. Lett. 1996, 76, 1039–1042. (38) Mazziotti, D. A. Phys. Rev. A 1998, 57, 4219–4234. (39) Mukherjee, D.; Kutzelnigg, W. J. Chem. Phys. 2001, 114, 2047–2061. (40) Kutzelnigg, W.; Mukherjee, D. J. Chem. Phys. 2004, 120, 7350–7368. (41) Mazziotti, D. A. Phys. Rev. Lett. 2006, 97, 143002. (42) Alcoba, D. R.; Valdemoro, C.; Tel, L. M.; Prez-Romero, E. International Journal of Quantum Chemistry 2009, 109, 3178–3190. (43) Valdemoro, C.; Alcoba, D. R.; Tel, L. M.; Pérez-Romero, E. Int. J. Quantum Chem. 2009, 109, 2622–2638. (44) Mazziotti, D. A. Phys. Rev. Lett. 2012, 108, 263002.

32


Page 32 of 39

Page 33 of 39

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


(45) Mazziotti, D. A. Phys. Rev. A 2012, 85, 062507. (46) Verstichel, B.; van Aggelen, H.; Van Neck, D.; Ayers, P. W.; Bultinck, P. J. Chem. Phys. 2010, 132, 114113. (47) van Aggelen, H.; Verstichel, B.; Bultinck, P.; Neck, D. V.; Ayers, P. W.; Cooper, D. L. J. Chem. Phys. 2011, 134, 054115. (48) Nakata, M.; Anderson, J. S. M. AIP Adv. 2012, 2, 032125. (49) Poelmans, W.; Raemdonck, M. V.; Verstichel, B.; Baerdemacker, S. D.; Torre, A.; Lain, L.; Massaccesi, G. E.; Alcoba, D. R.; Bultinck, P.; Neck, D. V. J. Chem. Theory Comput. 2015, 11, 4064–4076. (50) Burer, S.; Monteiro, R. D. Math. Prog. 2003, 95, 329–357. (51) Mazziotti, D. A. Phys. Rev. Lett. 2004, 93, 213001. (52) Mazziotti, D. A. ESAIM-Math. Model. Num. 2007, 41, 249–259. (53) Povh, J.; Rendl, F.; Wiegele, A. Computing 2006, 78, 277. (54) Malick, J.; Povh, J.; Rendl, F.; Wiegele, A. SIAM J. Optim. 2009, 20, 336. (55) Mazziotti, D. A. Phys. Rev. Lett. 2011, 106, 083001. (56) Gidofalvi, G.; Mazziotti, D. A. J. Chem. Phys. 2008, 129, 134108. (57) Pelzer, K.; Greenman, L.; Gidofalvi, G.; Mazziotti, D. A. J. Phys. Chem. A 2011, 115, 5632–5640. (58) Greenman, L.; Mazziotti, D. A. J. Chem. Phys. 2010, 133, 164110. (59) Sinitskiy, A. V.; Greenman, L.; Mazziotti, D. A. J. Chem. Phys. 2010, 133, 014104.

33



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

(60) Turney, J. M.; Simmonett, A. C.; Parrish, R. M.; Hohenstein, E. G.; Evangelista, F. A.; Fermann, J. T.; Mintz, B. J.; Burns, L. A.; Wilke, J. J.; Abrams, M. L.; Russ, N. J.; Leininger, M. L.; Janssen, C. L.; Seidl, E. T.; Allen, W. D.; Schaefer, H. F.; King, R. A.; Valeev, E. F.; Sherrill, C. D.; Crawford, T. D. WIRES Comput. Mol. Sci. 2012, 2, 556–565. (61) Erdahl, R. M. Int. J. Quantum Chem. 1978, 13, 697–718. (62) Whitten, J. L. J. Chem. Phys. 1973, 58, 4496–4501. (63) Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. J. Chem. Phys. 1979, 71, 3396–3402. (64) Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Chem. Phys. Lett. 1993, 208, 359–363. (65) Vahtras, O.; Almlöf, J.; Feyereisen, M. W. Chem. Phys. Lett. 1993, 213, 514–518. (66) Beebe, N. H. F.; Linderberg, J. Int. J. Quantum Chem. 1977, 12, 683–705. (67) Roeggen, I.; Wisloff-Nilssen, E. Chem. Phys. Lett. 1986, 132, 154–160. (68) Koch, H.; de Meras, A. S.; Pedersen, T. B. J. Chem. Phys. 2003, 118, 9481–9484. (69) Aquilante, F.; Pedersen, T. B.; Lindh, R. J. Chem. Phys. 2007, 126, 194106. (70) Hachmann, J.; Cardoen, W.; Chan, G. K.-L. J. Chem. Phys. 2006, 125, 144101. (71) Mayer, J. E. Phys. Rev. 1955, 100, 1579–1586. (72) Husimi, K. Proceedings of the Physico-Mathematical Society of Japan. 3rd Series 1940, 22, 264–314. (73) Löwdin, P.-O. Phys. Rev. 1955, 97, 1474–1489. (74) Coleman, A. J. Rev. Mod. Phys. 1963, 35, 668–686.

34


Page 34 of 39

Page 35 of 39

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


(75) Braams, B. J.; Percus, J. K.; Zhao, Z. Reduced-Density-Matrix Mechanics: With Application to Many-Electron Atoms and Molecules; John Wiley & Sons, Inc., 2007; pp 93–101. (76) Mazziotti, D. A. Phys. Rev. A 2005, 72, 032510. (77) Gidofalvi, G.; Mazziotti, D. A. Phys. Rev. A 2005, 72, 052505. (78) Pérez-Romero, E.; Tel, L. M.; Valdemoro, C. Int. J. Quantum Chem. 1997, 61, 55. (79) van Aggelen, H.; Verstichel, B.; Bultinck, P.; Neck, D. V.; Ayers., P. W. J. Chem. Phys. 2012, 136, 014110. (80) Shepard, R. In Advances in Chemical Physics; Lawley, K. P., Ed.; John Wiley & Sons, Inc.: New York, 1987; Vol. 69; pp 63–200. (81) Werner, H.-J. In Advances in Chemical Physics; Lawley, K. P., Ed.; John Wiley & Sons, Inc.: New York, 1987; Vol. 69; pp 1–62. (82) McWeeny, R.; Sutcliffe, B. Comput. Phys. Rep. 1985, 2, 219 – 278. (83) Olsen, J.; Yeager, D. L.; Jø rgensen, P. In Advances in Chemical Physics; Prigogine, I., Rice, S. A., Eds.; John Wiley & Sons, Inc.: New York, 1983; Vol. 54; pp 1–176. (84) Dalgaard, E.; Jørgensen, P. J. Chem. Phys. 1978, 69, 3833–3844. (85) Galin, C.; Schmidt, M. W.; Gordon, M. S. Theor. Chem. Acc. 1997, 97, 88–95. (86) Birks, J. B. Photophysics of aromatic molecules; Wily-Interscience: London, 1970. (87) Schiedt, J.; Weinkauf, R. Chemical physics letters 1997, 266, 201–205. (88) Sabbatini, N.; Indelli, M.; Gandolfi, M.; Balzani, V. The Journal of Physical Chemistry 1982, 86, 3585–3591.

35



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

(89) Burgos, J.; Pope, M.; Swenberg, C. E.; Alfano, R. Phys. Status Solidi B 1977, 83, 249–256. (90) Plasser, F.; Paˇsalić, H.; Gerzabek, M. H.; Libisch, F.; Reiter, R.; Burgdörfer, J.; M¨ uller, T.; Shepard, R.; Lischka, H. Angew. Chem. Int. Ed. 2013, 52, 2581–2584. (91) Horn, S.; Lischka, H. J. Chem. Phys. 2015, 142, 054302. (92) Horn, S.; Plasser, F.; M¨ uller, T.; Libisch, F.; Burgdörfer, J.; Lischka, H. Theor. Chem. Acc. 2014, 133, 1511. (93) Dunning Jr, T. H. J. Chem. Phys. 1970, 53, 2823–2833. (94) Dunning Jr, T. H.; Hay, P. J. In Methods in Electronic Structure Theory; Schaefer, H. F., Ed.; Plenum: New York, 1977; Vol. 2. (95) Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. J. Chem. Phys. 1997, 106, 9639–9646. (96) Hajgató, B.; Szieberth, D.; Geerlings, P.; De Proft, F.; Deleuze, M. J. Chem. Phys. 2009, 131, 224321. (97) Wu, C.-S.; Chai, J.-D. J. Chem. Theory Comput. 2015, 11, 20032011. (98) Valiev, M.; Bylaska, E.; Govind, N.; Kowalski, K.; Straatsma, T.; Dam, H. V.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T.; de Jong, W. Comput. Phys. Commun. 2010, 181, 1477 – 1489. (99) Purvis, G. D.; Bartlett, R. J. J. Chem. Phys. 1982, 76, 1910–1918. (100) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479–483. (101) Levchenko, S. V.; Krylov, A. I. J. Chem. Phys. 2004, 120, 175–185.

36


Page 36 of 39

Page 37 of 39

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


(102) Shao, Y.; Gan, Z.; Epifanovsky, E.; Gilbert, A. T. B.; Wormit, M.; Kussmann, J.; Lange, A. W.; Behn, A.; Deng, J.; Feng, X.; Ghosh, D.; Goldey, M.; Horn, P. R.; Jacobson, L. D.; Kaliman, I.; Khaliullin, R. Z.; K´ us, T.; Landau, A.; Liu, J.; Proynov, E. I.; Rhee, Y. M.; Richard, R. M.; Rohrdanz, M. A.; Steele, R. P.; Sundstrom, E. J.; Woodcock III, H. L.; Zimmerman, P. M.; Zuev, D.; Albrecht, B.; Alguire, E.; Austin, B.; Beran, G. J. O.; Bernard, Y. A.; Berquist, E.; Brandhorst, K.; Bravaya, K. B.; Brown, S. T.; Casanova, D.; Chang, C.-M.; Chen, Y.; Chien, S. H.; Closser, K. D.; Crittenden, D. L.; Diedenhofen, M.; DiStasio Jr., R. A.; Dop, H.; Dutoi, A. D.; Edgar, R. G.; Fatehi, S.; Fusti-Molnar, L.; Ghysels, A.; Golubeva-Zadorozhnaya, A.; Gomes, J.; Hanson-Heine, M. W. D.; Harbach, P. H. P.; Hauser, A. W.; Hohenstein, E. G.; Holden, Z. C.; Jagau, T.-C.; Ji, H.; Kaduk, B.; Khistyaev, K.; Kim, J.; Kim, J.; King, R. A.; Klunzinger, P.; Kosenkov, D.; Kowalczyk, T.; Krauter, C. M.; Lao, K. U.; Laurent, A.; Lawler, K. V.; Levchenko, S. V.; Lin, C. Y.; Liu, F.; Livshits, E.; Lochan, R. C.; Luenser, A.; Manohar, P.; Manzer, S. F.; Mao, S.-P.; Mardirossian, N.; Marenich, A. V.; Maurer, S. A.; Mayhall, N. J.; Oana, C. M.; Olivares-Amaya, R.; O’Neill, D. P.; Parkhill, J. A.; Perrine, T. M.; Peverati, R.; Pieniazek, P. A.; Prociuk, A.; Rehn, D. R.; Rosta, E.; Russ, N. J.; Sergueev, N.; Sharada, S. M.; Sharmaa, S.; Small, D. W.; Sodt, A.; Stein, T.; St¨ uck, D.; Su, Y.-C.; Thom, A. J. W.; Tsuchimochi, T.; Vogt, L.; Vydrov, O.; Wang, T.; Watson, M. A.; Wenzel, J.; White, A.; Williams, C. F.; Vanovschi, V.; Yeganeh, S.; Yost, S. R.; You, Z.-Q.; Zhang, I. Y.; Zhang, X.; Zhou, Y.; Brooks, B. R.; Chan, G. K. L.; Chipman, D. M.; Cramer, C. J.; Goddard III, W. A.; Gordon, M. S.; Hehre, W. J.; Klamt, A.; Schaefer III, H. F.; Schmidt, M. W.; Sherrill, C. D.; Truhlar, D. G.; Warshel, A.; Xua, X.; Aspuru-Guzik, A.; Baer, R.; Bell, A. T.; Besley, N. A.; Chai, J.-D.; Dreuw, A.; Dunietz, B. D.; Furlani, T. R.; Gwaltney, S. R.; Hsu, C.P.; Jung, Y.; Kong, J.; Lambrecht, D. S.; Liang, W.; Ochsenfeld, C.; Rassolov, V. A.; Slipchenko, L. V.; Subotnik, J. E.; Van Voorhis, T.; Herbert, J. M.; Krylov, A. I.;

37



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

Gill, P. M. W.; Head-Gordon, M. Mol. Phys. 2015, 113, 184–215. ¨ Fosso-Tande, J.; Chen, P.; White, J. L.; Allen, T. L.; Cherian, J.; (103) G¨ unaydin-Sen, O.; Tokumoto, T.; Lahti, P. M.; McGill, S.; Harrison, R. J.; Musfeldt, J. L. J. Chem. Phys. 2011, 135, 241101. (104) Minato, M.; Lahti, P. M. J. Am. Chem. Soc 1997, 119, 2187–2195. (105) Nimura, S.; Kikuchi, O.; Ohana, T.; Yabe, A.; Kaise, M. Chem. Lett. 1996, 25, 125– 126. (106) Yanai, T.; Chan, G. K.-L. J. Chem. Phys. 2006, 124, 194106. (107) Yanai, T.; Chan, G. K.-L. J. Chem. Phys. 2007, 127, 104107. (108) Neuscamman, E.; Yanai, T.; Chan, G. K.-L. J. Chem. Phys. 2009, 130, 124102. (109) Kurashige, Y.; Yanai, T. J. Chem. Phys. 2011, 135, 094104. (110) Saitow, M.; Kurashige, Y.; Yanai, T. J. Chem. Phys. 2013, 139, 044118. (111) Evangelista, F. A. J. Chem. Phys. 2014, 141, 054109. (112) Li, C.; Evangelista, F. A. J. Chem. Theory Comput. 2015, 11, 2097–2108, PMID: 26574413. (113) Colmenero, F.; Valdemoro, C. Phys. Rev. A 1993, 47, 979–985. (114) Mazziotti, D. A. Chem. Phys. Lett. 1998, 289, 419 – 427. (115) Kutzelnigg, W.; Mukherjee, D. J. Chem. Phys. 1999, 110, 2800–2809.

38


Page 38 of 39

Page 39 of 39


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

Large-Scale Variational Two-Electron Reduced ... - ACS Publications

Recommend Documents