Analytic Gradients for Complete Active Space Pair-Density Functional

Subscriber access provided by READING UNIV

Article

Analytic Gradients for Complete Active Space Pair-Density Functional Theory Andrew M. Sand, Chad E. Hoyer, Kamal Sharkas, Katherine Marie Kidder, Roland Lindh, Donald G. Truhlar, and Laura Gagliardi J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00967 • Publication Date (Web): 06 Dec 2017 Downloaded from http://pubs.acs.org on December 24, 2017

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

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

Page 1 of 42 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

Analytic Gradients for Complete Active Space Pair-Density Functional Theory Andrew M. Sand,† Chad E. Hoyer,† Kamal Sharkas,† Katherine M. Kidder,† Roland Lindh,‡ Donald G. Truhlar,∗,† and Laura Gagliardi∗,† †Department of Chemistry, Chemical Theory Center, and the Minnesota Supercomputing Institute, The University of Minnesota, Minneapolis, MN 55455 USA ‡Department of Chemistry - ˚ Angström, The Theoretical Chemistry Programme, Uppsala University, P.O. Box 518, SE-751 20 Uppsala, Sweden E-mail: [email protected]; [email protected]

Abstract Analytic gradient routines are a desirable feature for quantum mechanical methods, allowing for efficient determination of equilibrium and transition state structures and several other molecular properties. In this work, we present analytical gradients for multiconfiguration pair-density functional theory (MC-PDFT) when used with a state-specific complete active space self-consistent field reference wave function. Our approach constructs a Lagrangian that is variational in all wave function parameters. We find that MC-PDFT locates equilibrium geometries for several small- to mediumsized organic molecules that are similar to those located by complete active space second-order perturbation theory but that are obtained with decreased computational cost.

1

ACS Paragon Plus Environment

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

1

Introduction

Energy gradients play an important role in quantum chemistry. They are needed to determine equilibrium structures, transition state structures, trajectories in reaction dynamics, and many molecular properties such as nuclear magnetic resonance shifts. 1,2 Gradients can be computed numerically (by using finite differences) or analytically. The numerical approach, however, requires a costly number of single-point calculations, and its accuracy must be checked with respect to the finite difference step size. Thus, analytic approaches, 3–14 although more mathematically-involved, are preferred due to lower cost and higher accuracy. Recent advancements in CASPT2 methodologies have enabled analytic gradient calculations for both partially-contracted 15 and fully-contracted 16 formalisms. An accurate modeling of chemical systems requires quantum mechanical methods capable of describing static and dynamic correlation. 17–19 Methods which employ a single Slater determinant often do not give satisfactory results for strongly-correlated systems; the use of multiconfiguration wave function approaches 20 is desirable for such systems. However, these methods usually scale poorly with system size. Multiconfiguration pair-density functional theory 21,22 (MC-PDFT) is an affordable and accurate type of density functional theory that adds dynamic correlation energy to multiconfiguration self-consistent field (MCSCF) wave functions using a density functional called an on-top functional. The total electronic energy is expressed as a sum of the MCSCF kinetic energy and a functional of both the one-electron density and the two-electron on-top pair density of a reference wave function. Whereas the one-electron density is the probability of finding an electron at a point in space, the two-electron on-top pair density is related to the probability of finding two electrons on top of each other at a point in space. 23–27 The advantage MC-PDFT has over other post-MCSCF methods is that MC-PDFT can calculate the dynamic correlation energy correction at a cost that scales as a single iteration of the reference MCSCF calculation. 28 In contrast, post-MCSCF wave function methods such as complete active space second-order perturbation theory (CASPT2), 29 multireference 2


Page 2 of 42

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


Møller-Plesset perturbation theory (MRMP) 30–32 and multireference configuration interaction (MRCI), 33 while often providing a good description of correlation energy for smallersized systems, scale poorly with system size, which limits applications to many medium and large-sized systems. MC-PDFT has been shown to be able to predict energies with an accuracy similar to CASPT2 at a much lower computational cost. 21,22,28,34–39 In this work, we present the formulation and first application of analytic gradients for MCPDFT, using state-specific (SS-) CASSCF wave functions as references (SS-CAS-PDFT). However, SS-CAS-PDFT is not variational; the orbitals and CI expansion coefficients, while variationally optimized to minimize the energy for the reference CASSCF calculation, are not variationally optimized at the SS-CAS-PDFT level. This means the Helmann-Feynman theorem 40–43 does not hold, and coupled-perturbed MCSCF equations 7,44–47 must be solved. However, because the SS-CASSCF reference wave function is variational, it is possible to construct a Lagrangian that is variational. 5,8,48 This allows us to circumvent the use of coupled-perturbed MCSCF. Our paper is organized as follows. First, we provide a brief background on MC-PDFT theory. We next provide the derivation of the equations used to compute MC-PDFT gradients analytically. Then, we demonstrate the use of these equations by performing benchmark geometry optimizations on ammonia and a set of nine small- to medium-sized organic molecules. Finally, we optimize geometries for the highly multiconfigurational cases of ozone and benzyne isomers.

2

Theory

We begin by giving a brief background on MC-PDFT theory, then we discuss the use of Lagrange multipliers to address the non-variationality of the SS-CAS-PDFT energy, and finally we present the analytic form of the MC-PDFT gradient. Throughout this section we use the indices p, q, r, s, t, ... to refer to general molecular orbitals. We also make use of the

3



Page 4 of 42

following excitation operators: Eˆpq = a ˆ†pα a ˆqα + a ˆ†pβ a ˆqβ

(1)

− Eˆpq = Eˆpq − Eˆqp

(2)

eˆpqrs = Eˆpq Eˆrs − δqr Eˆps

(3)

ˆpα are second-quantized creation and annihilation operators, respectively, where a ˆ†pα and a and α and β are spin quantum numbers.

2.1

MC-PDFT theory

The energy in MC-PDFT can be expressed as 21,28

EPDFT = Vnn +

X

hpq Dpq +

pq

1X gpqst Dpq Dst + Eot [ρ, Π, ρ0 , Π0 ] 2 pqst

(4)

where hpq and gpqst are the one- and two-electron integrals: Z hpq = Z Z gpqst =

φ∗p (r)h(r)φq (r)dr

φ∗s (r1 )φt (r1 )

1 φ∗ (r2 )φq (r2 )dr1 dr2 |r1 − r2 | p

(5)

(6)

where the one-electron operator h(r) accounts for both electronic kinetic energy and electronnuclear potential energy: h(r) =

−∇2 X ZA − , 2 |r − rA | A

(7)

Molecular orbitals are indicated by φi (we assume real orbitals to keep the notation simple), while ρ and Π are the electronic density and the on-top pair density, respectively, and ρ0 and Π0 are their derivatives. These densities can be expressed in terms of the orbitals, the

4


Page 5 of 42 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-body density matrix D, and the two-body density matrix d:

ρ(r) =

X

φp (r)φq (r)Dpq

(8)

φp (r)φq (r)φs (r)φt (r)dpqst

(9)

φ0p (r)φq (r) + φp (r)φ0q (r) Dpq

(10)

pq

Π(r) =

X pqst

ρ0 (r) =

X pq

Π0 (r) =

X

φ0p (r)φq (r)φs (r)φt (r) + φp (r)φ0q (r)φs (r)φt (r)

pqst

+ φp (r)φq (r)φ0s (r)φt (r) + φp (r)φq (r)φs (r)φ0t (r) dpqst (11) The current generation of on-top energy functionals Eot [ρ, Π, ρ0 , Π0 ] is formed by translating existing local density approximation and generalized gradient approximation KohnSham (KS) density functionals. 21,49 The parent KS density functionals depend on the spin-up and spin-down electron densities ρα and ρβ as well as their derivatives ρ0α and ρ0β . We express these KS functionals as Exc [ρα , ρβ , ρ0α , ρ0β ]. Our original translation scheme 21 defines the on-top energy functional (with no dependence on Π0 ) as

Eot [ρ, Π, ρ0 ] = Exc [˜ ρα , ρ˜β , ρ˜0α , ρ˜0β ]

5


(12)


Page 6 of 42

where the translated densities are defined as

ρ˜α (r) =

ρ˜β (r) =

ρ˜0α (r) =

ρ˜0β (r) =

    ρ(r) (1 + ζt (r)) R(r) ≤ 1 2    ρ(r) R(r) > 1 2     ρ(r) (1 − ζt (r)) R(r) ≤ 1 2    ρ(r) R(r) > 1 2   0   ρ (r) (1 + ζt (r)) R(r) ≤ 1 2  0   ρ (r) R(r) > 1 2   0   ρ (r) (1 − ζt (r)) R(r) ≤ 1 2  0   ρ (r)

(13)

(14)

(15)

(16)

R(r) > 1

2

where ζt (r) =

p

R(r) =

1 − R(r)

Π(r) [ρ(r)/2]2

(17) (18)

and

ρ(r) = ρα (r) + ρβ (r)

(19)

ρ0 (r) = ρ0α (r) + ρ0β (r).

(20)

This scheme is denoted by “t”. We have also developed a fully-translated (“ft”) scheme 49 in which the on-top functional depends on Π0 . For simplicity, the main text will only treat the “t” functional case, but we present an appendix that gives the changes for the “ft” functional case.

6


Page 7 of 42 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.2

SS-CAS-PDFT energy Lagrangian

The Hellmann-Feynman theorem 40–43 shows that if a wave function is variationally optimized in all of its parameters, the first-order energy response to a perturbation depends only on the expectation value of the derivative of the Hamiltonian operator with respect to the perturbation. ˆ d ˆ var i = hΨvar | dH |Ψvar i hΨvar |H|Ψ dλ dλ

(21)

In this work, we will only use SS-CASSCF wave functions as references for MC-PDFT calculations. In a SS-CASSCF wave function, the energy is parameterized in the following way: 50 ˆ

ˆ

ˆ −ˆκ e−P |0i E0 = h0|eP eκˆ He

(22)

where κ ˆ represents the orbital rotation operators

κ ˆ=

X

− κpq Eˆpq

(23)

p 1 2  h h i i  2Π(r) 1   1 − ζ (r) − φ (r)φ (r) R(r) ≤ 1 t p q 2 ζt (r)[ρ(r)]2 ∂ ρ˜β (r) = h i  ∂Dpq   1 φp (r)φq (r) R(r) > 1 2  h ih i h i   2ρ0 (r)Π(r) 1 0 0  0 1 + ζ (r) φ (r)φ (r) + φ (r)φ (r) + φ (r)φ (r) t q p p q p q ∂ ρ˜α (r) 2 ζt (r)[ρ(r)]3 = h i  ∂Dpq   1 φ0p (r)φq (r) + φp (r)φ0q (r) 2  h ih i h i   2ρ0 (r)Π(r) 1 0 0  0 1 − ζ (r) φ (r)φ (r) + φ (r)φ (r) − φ (r)φ (r) t q p p q ∂ ρ˜β (r) p q 2 ζt (r)[ρ(r)]3 = h i  ∂Dpq   1 φ0p (r)φq (r) + φp (r)φ0q (r) 2

10


(40)

(41)

R(r) ≤ 1 (42) R(r) > 1 R(r) ≤ 1 (43) R(r) > 1

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


Likewise, we define the two-electron on-top potential as

vpqst =

∂Eot [ρ, Π, ρ0 ] . ∂dpqst

(44)

Similarly, by using the translation scheme in eqs. (13-16), we can write the two-electron on-top potential in terms of the derivatives of the KS density functional:

vpqst =

∂Exc [˜ ρα , ρ˜β , ρ˜0α , ρ˜0β ] ∂ ρ˜α ∂Exc [˜ ρα , ρ˜β , ρ˜0α , ρ˜0β ] ∂ ρ˜β + ∂ ρ˜α ∂dpqst ∂ ρ˜β ∂dpqst 0 0 0 ∂Exc [˜ ρα , ρ˜β , ρ˜0α , ρ˜0β ] ∂ ρ˜0β ∂Exc [˜ ρα , ρ˜β , ρ˜α , ρ˜β ] ∂ ρ˜α + (45) + ∂ ρ˜0α ∂dpqst ∂ ρ˜0β ∂dpqst

The derivatives of the translated densities and gradients with respect to the two-body density matrix are obtained by taking derivatives of eqs. (13-16) over a grid of points:    

−1 ρ(r)ζt (r)

h

i φp (r)φq (r)φs (r)φt (r) R(r) ≤ 1

∂ ρ˜α (r) =  ∂dpqst  0 R(r) > 1  h i  1   φ (r)φ (r)φ (r)φ (r) R(r) ≤ 1 p q s t ρ(r)ζt (r) ∂ ρ˜β (r) =  ∂dpqst  0 R(r) > 1  h i  −ρ0 (r)   φ (r)φ (r)φ (r)φ (r) R(r) ≤ 1 0 p q s t [ρ(r)]2 ζt (r) ∂ ρ˜α (r) =  ∂dpqst  0 R(r) > 1  h i  ρ0 (r)   0 φ (r)φ (r)φ (r)φ (r) R(r) ≤ 1 p q s t ∂ ρ˜β (r) [ρ(r)]2 ζt (r) =  ∂dpqst  0 R(r) > 1

(46)

(47)

(48)

(49)

By combining these results, eq. 34 can be written succinctly as ∂EPDFT = 2(Fxy − Fyx ) ∂κxy

11


(50)


Page 12 of 42

where Fxy are elements of the generalized SS-CAS-PDFT Fock matrix:

Fxy =

X

(hpy + Vpy )Dpx +

p

X

(gpqry Dpq Drx + 2vpqry dpqrx ).

(51)

pqs

Next, we consider the CI state transfer response: X ∂ 1X ∂EPDFT 0 = Vnn + hpq Dpq + gpqst Dpq Dst + Eot [ρ, Π, ρ ] . ∂PI ∂PI 2 pqst pq

(52)

The derivative of the one-body and two-body density matrix with respect to a CI coefficient are given by ∂Dpq = hi|Eˆpq |0i + h0|Eˆpq |ii ∂Pi

(53)

∂dpqst = hi|ˆ epqst |0i + h0|ˆ epqst |ii ∂Pi

(54)

where |0i is the SS-CASSCF reference wavefunction and |ii is a CSF. The derivative of the on-top energy functional with respect to a CI coefficient is given as ∂Eot [ρ, Π, ρ0 ] X ∂Eot [ρ, Π, ρ0 ] ∂Dpq X ∂Eot [ρ, Π, ρ0 ] ∂dpqst = + ∂Pi ∂D ∂P ∂dpqst ∂Pi pq i pq pqst X ∂dpqst ∂Dpq X vpqst Vpq + = ∂Pi ∂Pi pqst pq

(55) (56)

These results allow us to write eq. 52 as ∂EPDFT = 2hi|Fˆ1 + Fˆ2 − E|0i ∂Pi

(57)

where Fˆ1 =

X

hpq + Vpq +

pq

X

gpqst Dst Eˆpq

(58)

st

Fˆ2 =

X

vpqst eˆpqst

pqst

12


(59)

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


E = h0|Fˆ1 + Fˆ2 |0i.

2.4

(60)

Response equations

The Lagrange multipliers are obtained by solving the linear system of equations given in eq. 29. In lieu of the explicit construction and direct diagonalization of the SS-CASSCF Hessian matrix, we solve eq. 29 via the iterative preconditioned conjugate gradient (PCG) algorithm, 52 during which Hessian and trial vector multiplications are performed on-the-fly. In our approach, the standard CASSCF preconditioner is used. 48 Further details about the PCG procedure are available in refs. 8,48

2.5

SS-CAS-PDFT nuclear gradients

The gradient of the SS-CAS-PDFT energy, upon determination of the Lagrange multipliers, is given by X X X X ∂ X λ ∂ X λ dEPDFT (λ) λ λ = EPDFT + zxy ( hpq Dpq + gpqst dpqst )+ zI ( hpq Dpq + gpqst dpqst ) dλ ∂κ ∂P xy I xy pq pqst pq pqst I (61) λ where hλpq and gpqst are derivative integrals given by

dhpq dλ i ∂hpq 1 X h ∂Spw ∂Swq = − hwq + hpw ∂λ 2 w ∂λ ∂λ

hλpq =

(62) (63)

dgpqst dλ i ∂gpqst 1 X h ∂Spw ∂Sqw ∂Ssw ∂Stw = − gwqst + gpwst + gpqwt + gpqsw ∂λ 2 w ∂λ ∂λ ∂λ ∂λ

λ gpqst =

(64) (65)

where Sxy is an element of the orbital overlap matrix. The terms involving derivatives of the overlap matrix arise due to the response of the basis set to the perturbation. This is often (λ)

referred to as the ‘connection’ or ‘renormalization’ contribution. EPDFT is the SS-CAS-PDFT

13



Page 14 of 42

energy expression evaluated with the derivatives of energy operators and functionals: (λ)

EPDFT =

dVnn X λ 1X λ dEot [ρ, Π, ρ0 ] + hpq Dpq + gpqst Dpq Dst + dλ 2 pqst dλ pq

(66)

Like the one- and two-electron integral derivatives, the derivative of the on-top functional also contributes a ‘renormalization’ contribution: i ∂E [ρ, Π, ρ0 ] X dEot [ρ, Π, ρ0 ] X λ h X ot = Sxy Vxp Dyp + vxpqs dypqs + . dλ ∂λ xy p pqs

(67)

The final term in eq. 67 is evaluated using standard DFT techniques, 53 using the translated densities in the evaluation of the derivative of the KS-DFT functional. This treatment explicitly includes the derivatives of grid weights in the evaluation of the functional gradient. The remaining derivatives in eq. 61 can be evaluated using eqs. 35-36 and eqs. 53-54, and upon rearrangement we obtain dEPDFT dVnn 1 X ∂Spq eff X ∂hpq eff X ∂gpqst eff ∂Eot [ρ, Π, ρ0 ] = − Fpq + Dpq + dpqst + dλ dλ 2 pq ∂λ ∂λ ∂λ ∂λ pq pqst

(68)

where we have introduced an effective one-body density matrix

eff ˘ pq + D ¯ pq Dpq = Dpq + D X ˘ pq = D (Dsq zps − Dps zsq )

(69) (70)

s

¯ pq = D

X

zI (hI|Eˆpq |0i + h0|Eˆpq |Ii,

I

14


(71)

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


an effective two-body density matrix

˘ ¯ deff pqst = dpqst + dpqst + dpqst X d˘pqst = (duqst zpu − dpust zqu + dpqut zsu − dpqsu ztu )

(72) (73)

u

d¯pqst =

X

zI (hI|ˆ epqst |0i + h0|ˆ epqst |Ii,

(74)

I

and an effective Fock matrix

eff Fpq =

X

(hpt +vpt )Dqt +

t

X

X X X ˘ qt +D ¯ qt )+ gprst Dqr Dst + vprst dqrst + hpt (D gprst (d˘qrst +d¯qrst ).

rst

rst

t

rst

(75) Equations 69-74 have been previously reported in the derivation of other Lagrangian-based gradient formulations. 8,15

3

Methods

NH3

pyrrole

HCN

CH2O

HCCH

oxirane

acrolein

butadiene

pyridine

maleic anhydride

Figure 1: Systems from the SE47 database included in this study.

A set of ten molecules from the SE47 database of structures 54 was utilized in the first 15



4 3

3

5 2

6 1

2

o-benzyne

ozone

4

5

6

5

6

3

2

4

1

m-benzyne

1

p-benzyne

Figure 2: Additional systems considered in this study. Table 1: Active space selection [denoted (electrons, orbitals)] and symmetry employed for each system. The index is used in the labeling of Figs. 3 and 4. Index Active Space Symmetry Constraints NH3 1 (6,6) Cs HCN 2 (8,8) C2v CH2 O 3 (12,9) C2v HCCH 4 (10,10) D2h oxirane 5 (10,10) C2v pyrrole 6 (6,5) C2v acrolein 7 (4,4) Cs butadiene 8 (4,4) C2h pyridine 9 (6,6) C2v maleic anhydride 10 (8,7) C2v ozone 11 (12,8) C2v o-benzyne 12 (8,8) C2v m-benzyne 13 (8,8) C2v p-benzyne 14 (8,8) D2h 16


Page 16 of 42

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


stage of this study (Fig. 1) and in the second stage, we considered four additional cases (Fig. 2): ozone and the three isomers of benzyne (ortho-, meta-, and para-). All CASSCF, CASPT2, Kohn-Sham density functional theory (KS-DFT), and MC-PDFT calculations were performed with a locally-modified version of Molcas 8.1 . 55 The same CASSCF reference was used for all multireference calculations on a given system, and the active space selection and symmetry employed are listed in Table 1. All tPBE on-top density functional was used for MC-PDFT calculations on all systems, and additional MC-PDFT calculations on the benzyne isomers used the tBLYP on-top density functional. 21 The fine integration grid was employed. An MC-PDFT numerical gradient implementation based on finite differences was also developed as part of this project. All analytically-determined structures were verified to agree with the numerically-determined structures within the precision reported for bond lengths (0.001 ˚ A) and bond angles (0.1 degrees) in the supporting information. We employed the cc-pVDZ and cc-pVTZ basis sets 56,57 for all calculations. The M diagnostic 58 was used to evaluate the degree of multireference character. The M diagnostic is given by nX SONO 1 2 − n(MCDONO) + n(MCUNO) + |n(SON Oj ) − 1| M= 2 j=1

(76)

where n(MCDONO), n(MCUNO), and n(SONOj ) are natural orbital occupation numbers for the most-correlated doubly-occupied natural orbital, the most-correlated unoccupied natural orbital, and singly-occupied natural orbitals, respectively, while nSONO is the number of singly-occupied orbitals. The distinctions of doubly-occupied, singly-occupied and unoccupied are made by examining the dominant configuration in the CASSCF wave function.

4

Results

A total of 14 different molecules were considered, and the resultant geometric parameters are available in Table S1. M diagnostic results are given in Table 2. Generally, systems exhibit17



Page 18 of 42

Table 2: M diagnostics 58 and orbital occupation numbers for the CASSCF equilibrium geometries. M values greater than 0.1 indicate significant multireference character. M

HONO occupation

LUNO occupation

cc-pVDZ cc-pVTZ cc-pVDZ cc-pVTZ cc-pVDZ cc-pVTZ NH3 0.023 0.022 1.977 1.978 0.022 0.023 HCN 0.066 0.064 1.934 1.936 0.066 0.064 CH2 O 0.065 0.063 1.936 1.939 0.067 0.064 HCCH 0.068 0.065 1.932 1.935 0.068 0.066 oxirane 0.036 0.034 1.965 1.966 0.036 0.035 pyrrole 0.076 0.074 1.925 1.926 0.077 0.0753 acrolein 0.107 0.103 1.897 1.900 0.110 0.106 butadiene 0.118 0.113 1.884 1.889 0.120 0.116 pyridine 0.103 0.101 1.898 1.900 0.105 0.102 maleic anhydride 0.104 0.100 1.904 1.907 0.111 0.107 ozone 0.256 0.240 1.764 1.779 0.275 0.259 o-benzyne 0.177 0.168 1.824 1.832 0.177 0.168 m-benzyne 0.346 0.331 1.654 1.669 0.346 0.331 p-benzyne 0.772 0.772 1.228 1.228 0.772 0.772 ing M values greater than 0.1 are considered to have significant multireference character. Acrolein, butadiene, pyridine, and maleic anhydride all were found to have M diagnostic values slightly larger than 0.1, and ozone and the three isomers of benzyne were found to have M diagnostic values significantly higher than 0.1. In order to evaluate the performance of each method, each converged structure was best-fit (via rigid rotation and translation) to the corresponding reference structure from the SE47 database by minimizing the root mean square (RMS) distance between the atomic centers. 59 We subdivide our results into two sections: systems with predominantly single-reference character (1-6) and systems with strong multireference character (7-14).

4.1

Systems with predominantly single-reference character

Although CAS-PDFT is a method designed for multireference calculations, it is important to test whether tPBE can determine accurate geometries for systems with predominantly single-reference character. Molecules 1—6 were found to have M diagnostic values smaller

18


Page 19 of 42

0 .0 3 5 tP B C A C A P B

R M S d i s p l a c e m e n t ( Å)

0 .0 3 0

E S S C F S P T 2 E

0 .0 2 5 0 .0 2 0 0 .0 1 5 0 .0 1 0 0 .0 0 5 0 .0 0 0 1

2

3

4

5

6

M o le c u le

Figure 3: Root mean square (RMS) displacements of the atomic centers of the converged structure relative to the SE47 reference structure. 54 The cc-pVDZ basis set was used.

0 .0 3 5 tP B C A C A P B

0 .0 3 0

RMS displacement ( Å)

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


0 .0 2 5

E S S C F S P T 2 E

0 .0 2 0 0 .0 1 5 0 .0 1 0 0 .0 0 5 0 .0 0 0 1

2

3

4

5

6

M o le c u le

Figure 4: RMS displacements of the atomic centers of the converged structure relative to the SE47 reference structure. 54 The cc-pVTZ basis set was used.

19



than 0.1, indicating low amounts of multireference character. Results using the cc-pVDZ basis set are given in Figure 3, and results using the cc-pVTZ basis set are given in Figure 4. Increasing the basis set size lowers the RMS displacement for tPBE, CASPT2, and PBE while the CASSCF RMS displacements become larger for several molecules. A total of 27 geometric parameters (16 bond lengths and 11 bond angles) were optimized in molecules 1—6. If we compare the tPBE, CASSCF, and CASPT2 results (which all employ the same active space), we see that the tPBE results are more similar to the CASPT2 results than to the CASSCF results for 22 of these data when using the cc-pVDZ basis set and 23 of these data when using the cc-pVTZ basis set. Additionally, if we view the tPBE and CASPT2 methods as “correcting” the CASSCF result, tPBE and CASPT2 agree on the direction of the change in the geometric parameter for 24 of the geometrical variables when using the cc-pVDZ and 25 of them when using the cc-pVTZ basis set. A comparison between the tPBE and PBE results shows that both methods predict very similar geometries for these species.

4.2

Systems with strong multireference character

0 .0 3 5 tP B C A C A P B

0 .0 3 0


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

E S S C F S P T 2 E

0 .0 2 5 0 .0 2 0 0 .0 1 5 0 .0 1 0 0 .0 0 5 0 .0 0 0 7

8

9

1 0

M o le c u le

Figure 5: Root mean square (RMS) displacements of the atomic centers of the converged structure relative to the SE47 reference structure. 54 The cc-pVDZ basis set was used. 20


Page 20 of 42

Page 21 of 42

0 .0 3 5 tP B C A C A P B

0 .0 3 0

RMS displacement ( Å)

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


E S S C F S P T 2 E

0 .0 2 5 0 .0 2 0 0 .0 1 5 0 .0 1 0 0 .0 0 5 0 .0 0 0 7

8

9

1 0

M o le c u le

Figure 6: RMS displacements of the atomic centers of the converged structure relative to the SE47 reference structure. 54 The cc-pVTZ basis set was used. Molecules 7—10 from the SE47 database all had M diagnostic values larger than 0.1 and are thus considered to have significant multireference character. Results using the cc-pVDZ basis set are given in Figure 5, and results using the cc-pVTZ basis set are given in Figure 6. Considerably larger discrepancies are seen between PBE and tPBE are seen in these cases, and tPBE shows a noticeable improvement over PBE. The selection of active spaces can possibly affect the quality of results for certain bonds and angles. For example, C-H bonding and antibonding orbitals are often not included in the active space of the largest systems due to active space size limitations. This is reflected in the resultant tPBE bond lengths–the C–H and N–H bonds tend to have larger errors than for the C–C, N–C, and C–O bonds. In the case of pyridine, a (6,6) full-π active space was selected. The errors in the tPBE bond lengths (cc-pVTZ) are about 0.001 ˚ A for the C–C and C–N bonds, but the C–H bond errors are near 0.010 ˚ A. A similar trend can be seen in the CASSCF results. CASPT2 appears less susceptible to active space limitations as its errors are quite similar across all bonds. Because CAS-PDFT does not improve or alter the reference wave function (it only supplies an energy correction), the reference wave function likely has the potential to exhibit greater influence over the quality of the CAS-PDFT result 21



Page 22 of 42

than in a method such as CASPT2, which computes a first-order correction to the reference wave function. A total of 47 geometric parameters (23 bond lengths and 24 bond angles) were optimized for molecules 7—10. The results obtained with tPBE are more similar to the CASPT2 results than to the CASSCF results for 37 of these data when using the cc-pVDZ basis set and 34 of these data when using the cc-pVTZ basis set. Additionally, tPBE and CASPT2 agree on the direction of the change in the geometric parameter for 39 of the geometrical variables when using the cc-pVDZ and 40 of them when using the cc-pVTZ basis set. In contrast with the predominantly single-reference results, tPBE gives noticeably more accurate results than PBE. Table 3: Mean signed error (MSE) and mean unsigned error (MUE) over the set of 39 bond lengths and 35 bond angles. Errors are calculated relative to the SE47 database. 54 cc-pVDZ

cc-pVTZ

tPBE CASSCF CASPT2

PBE

tPBE CASSCF CASPT2

PBE

MSE bonds lengths bond angles

0.014 -0.1

0.005 -0.1

0.013 -0.1

0.016 -0.3

0.006 0.0

-0.003 0.0

0.001 -0.3

0.008 0.0

MUE bonds lengths bond angles

0.014 0.4

0.008 0.6

0.013 0.6

0.016 0.7

0.007 0.3

0.008 0.6

0.003 0.5

0.009 0.4

The mean signed error (MSE) and mean unsigned error (MUE) over the set of 39 bond lengths and 35 bond angles from the SE47 database (molecules 1—10) for each method are given in Table 3. The errors for all methods are small. The MC-PDFT results obtained with tPBE show similar MSEs and MUEs for both basis sets. Notably, the tPBE results show an improvement upon the KS-DFT results obtained with PBE. Further, the tPBE results parallel the CASPT2 and PBE results in terms of performance in going from the smaller cc-pVDZ to the larger cc-pVTZ basis set, obtaining greater accuracy with the larger basis set. In contrast, the MUE for the CASSCF method does not show an improvement upon going to a larger basis set. 22


Page 23 of 42

Table 4: Calculated equilibrium geometries of ozone cc-pVDZ

cc-pVTZ

tPBE CASSCF CASPT2 PBE tPBE CASSCF CASPT2 PBE r(O-O) 1.277 1.290 1.290 1.279 1.274 1.282 1.282 1.276 a(O-O-O) 116.9 116.4 116.4 118.0 117.2 116.8 116.7 118.1 RMSD 0.001 0.008 0.008 0.008 0.003 0.003 0.003 0.008

Experiment 60 1.278 116.8

Ozone has long provided an important test for electronic structure methods as it requires the inclusion of both static and dynamic correlation effects in order to correctly predict its molecular properties. 61,62 The M diagnostics given in Table 2 reflect the very high multiconfigurational nature of ozone. Optimized structures are provided in Table 4. The geometries calculated by CAS-PDFT with the tPBE functional are found to closely match the experimentally-derived geometry 60 with improved accuracy over both CASSCF and PBE. In comparison to CASPT2, tPBE exhibits smaller errors at the cc-pVDZ level and similarly-small errors at the cc-pVTZ level.

0 .1 4 tP B E tB L Y P C A S S C F C A S P T 2 P B E B L Y P

0 .1 2 0 .1 0


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


0 .0 8 0 .0 6 0 .0 4 0 .0 2 0 .0 0

o -b e n z y n e

m -b e n z y n e

p -b e n z y n e

Figure 7: RMS displacements of the atomic centers of the converged benzyne structure relative to the reference structures. The cc-pVTZ basis set was used. The isomers of benzyne have been extensively studied with a variety of computational methods. 63–73 Previous theoretical studies 63–65 and the M diagnostics given in Table ?? indicate the following trend in the degree of biradical character: ortho- < meta- < para-. 23



0 .1 4

tP B E tB L Y P C A S S C F C A S P T 2 P B E B L Y P

0 .1 2 0 .1 0


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

0 .0 8 0 .0 6 0 .0 4 0 .0 2 0 .0 0

o -b e n z y n e

m -b e n z y n e

p -b e n z y n e

Figure 8: RMS displacements of the atomic centers of the converged benzyne structure relative to the reference structures. The cc-pVTZ basis set was used. The geometry optimization of m-benzyne is particularly challenging, with several methods predicting a bicyclic structure (indicated in Fig. 2 with a dashed line) with much less biradical character than the monocyclic structure. We have compared our computed structures with an experimentally-determined structure for o-benzyne 74 and theoretically-determined [RMR-CCSD(T)/6-311G(dp)] structures for m-benzyne and p-benzyne. 70 RMS displacements for the cc-PVDZ basis set and cc-PVTZ basis sets are given in Figs. 7 and 8, respectively, and complete structural information (coordinates, bond distances, and bond angles) are given in the supporting information. For the o-benzyne and p-benzyne structures, both CAS-PDFT functionals employed show good agreement with the reference structures, with RMS displacements very similar to the CASPT2 values. For p-benzyne, the most multiconfigurational molecule considered, the CAS-PDFT structures show significant improvement over the structures determined via KS-DFT with the corresponding functional. For the m-benzyne structure, CAS-PDFT with the tBLYP on-top functional predict a monocyclic structure in good agreement with the reference structure, with a deviation similar in magnitude to the predicted CASPT2 structure. The C1-C3 bond distance for tBLYP/cc-pVTZ was calculated to be 2.004 ˚ A. In contrast, CAS-PDFT with the tPBE 24


Page 24 of 42

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


on-top functional predicts a bicyclic m-benzyne structure, with a C1-C3 bond distance of 1.573 ˚ A. The CASPT2 C1-C3 distance of 2.063 ˚ A agrees best with the reference value of 2.108 ˚ A. The performance of each CAS-PDFT functional mirrors the performance of the corresponding KS-DFT functional; tPBE and PBE predict bicyclic structures and tBLYP and BLYP predict monocyclic structures. Despite the large deviation in the C1-C3 distance, the tPBE energy of the CASPT2/cc-pVTZ minimum-energy structure is only 3.12 kcal/mol higher in energy than the energy of the tPBE optimized structure. Thus, for p-benzyne, the most inherently multiconfigurational of the benzynes, both tPBE and tBLYP perform quite well, greatly improving over their KS-DFT counterparts.

5

Conclusions

In this work we have derived the working equations necessary to calculate analytical gradients for the CAS-PDFT method. We have used these routines to perform geometry optimizations on a set of ten predominantly single-reference molecules, including four molecules with significant (non-negligible) multireference character. We have found that CAS-PDFT with the tPBE functional performs comparably to CASPT2 on these predominantly-singlereferenced systems. For the four multiconfigruational systems considered, CAS-PDFT was similarly found to perform comparably to CASPT2, with an exception for tPBE applied to m-benzyne. However, good agreement was obtained with the use of the tBLYP functional. Throughout all systems considered, the CAS-PDFT results paralleled (but improved) the KS-DFT results with the corresponding untranslated functional; further improvement of CAS-PDFT geometries is also possible through the implementation of new on-top pair density functionals as the current generation of on-top functionals are translated from KS-DFT functionals which are often not the best functional for geometry optimizations. 75 In future work we will apply the analytical gradient treatment of CAS-PDFT to additional multireference systems, including transition states and excited states.

25



The current implementation for analytical CAS-PDFT gradients is restricted to statespecific CASSCF reference wave functions. In future work we plan to develop state-averaged CAS-PDFT analytic gradient routines.

6

Associated Content

6.1

Supporting Information

The supporting information is available free of charge on the ACS Publications website at DOI: 10.1021/ The Supporting information contains optimized Cartesian coordinates and some examples showing results for ftPBE.

7

Acknowledgment

This work was supported in part by the National Science Foundation by grants CHE-1464536. K.M.K. acknowledges support from summer research fellowship CHE-1359181. R.L. acknowledges support by the Swedish Research Council (Grant 2016-03398).

A

Appendix: Analytic gradients for fully-translated functionals

In the fully-translated (“ft”) functional scheme, the on-top functional depends on the derivative of the on-top pair density Π0 . The functional transformation becomes Eot [ρ, Π, ρ0 , Π0 ] =

26


Page 26 of 42

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


Exc [˜ ρα , ρ˜β , ρ˜0α , ρ˜0β ], where the ft-densities are written as 49

ρ˜α (r) =

   ρ(r)  (1 + ζt (r))  2   

R(r) < R0

ρ(r) (1 + ζf t (r)) R0 ≤ R(r) ≤ R1 2        ρ(r) R(r) > R1 2    ρ(r)  (1 − ζt (r)) R(r) ≤ R0    2 

ρ(r) (1 − ζf t (r)) R0 ≤ R(r) ≤ R1 2        ρ(r) R(r) > R1 2    ρ0 (r)  (1 + ζt (r)) + ρ(r) ζ 0 (r) R(r) ≤ R0  2 2 t    ρ˜0α (r) = ρ0 (r) (1 + ζf t (r)) + ρ(r) ζf0 t (r) R0 ≤ R(r) ≤ R1 2 2      0   ρ (r) R(r) > R1 2    ρ0 (r)  (1 − ζt (r)) − ρ(r) ζ 0 (r) R(r) ≤ R0  2 2 t    ρ˜0β (r) = ρ0 (r) (1 − ζf t (r)) − ρ(r) ζf0 t (r) R0 ≤ R(r) ≤ R1 2 2      0   ρ (r) R(r) > R1 2

ρ˜β (r) =

(77)

(78)

(79)

(80)

where ζt (r), ζf t (r),ζt0 (r) and ζf0 t (r) are defined as

ζt (r) =

p

1 − R(r)

(81)

ζf t (r) = A(R(r) − R1 )5 + B(R(r) − R1 )4 + C(R(r) − R1 )3 ζt0 (r) = −

1 R0 (r) 2 ζt (r)

(82) (83)

ζf0 t (r) = R0 (r)[5A(R(r) − R1 )4 + 4B(R(r) − R1 )3 + 3C(R(r) − R1 )2 ]

27


(84)


Page 28 of 42

with the following parameters:

R0 = 0.9

(85)

R1 = 1.15

(86)

A = −475.60656009

(87)

B = −379.47331922

(88)

C = −85.38149682

(89)

The gradient R0 (r) is written as

R0 (r) =

Π0 (r) Π(r)ρ0 (r) − [ρ(r)/2]2 [ρ(r)/2]3

(90)

The evaluation of the one-electron and two-electron on-top potentials Vpq (Eq. 38) and vpqrs (Eq. 44) require the derivatives of Eqs. 77-80 with respect to the one- and two-body density

28


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


matrices. Derivatives with respect to the one-body density matrix are given by    (1+ζt (r)) ∂ρ(r)  +  2 ∂Dpq   

ρ(r) ∂ζt (r)) 2 ∂Dpq

∂ ρ˜α (r) = (1+ζf t (r)) ∂ρ(r) + 2 ∂Dpq  ∂Dpq       1 ∂ρ(r)

R(r) < R0

ρ(r) ∂ζf t (r)) 2 ∂Dpq

R(r) > R1

2 ∂Dpq

   (1−ζt (r)) ∂ρ(r)  −  2 ∂Dpq   

ρ(r) ∂ζt (r)) 2 ∂Dpq

∂ ρ˜β (r) = (1−ζf t (r)) ∂ρ(r) − 2 ∂Dpq  ∂Dpq       1 ∂ρ(r)

R(r) < R0

ρ(r) ∂ζf t (r)) 2 ∂Dpq

R(r) > R1 ρ0 (r) ∂ζt (r) 2 ∂Dpq

∂ ρ˜0α (r) = (1+ζf t (r)) ∂ρ0 (r) + 2 ∂Dpq  ∂Dpq     0   1 ∂ρ (r)

+

ρ0 (r) ∂ζf t (r) 2 ∂Dpq

ρ(r) ∂ζt0 (r) 2 ∂Dpq

+

ζt0 (r) ∂ρ(r) 2 ∂Dpq

0

+

ρ(r) ∂ζf t (r) 2 ∂Dpq

+

ζf0 t (r) ∂ρ(r) 2 ∂Dpq

∂ ρ˜0β (r) ∂Dpq

=

(1−ζf t (r)) ∂ρ0 (r) 2 ∂Dpq 

ρ0 (r) ∂ζt (r) 2 ∂Dpq

−

R(r) < R0 R0 ≤ R(r) ≤ R1

(93)

R(r) > R1

2 ∂Dpq

   (1−ζt (r)) ∂ρ0 (r)  −  2 ∂Dpq   

(92)

R0 ≤ R(r) ≤ R1

2 ∂Dpq

   (1+ζt (r)) ∂ρ0 (r)  +  2 ∂Dpq   

(91)

R0 ≤ R(r) ≤ R1

−

ρ0 (r) ∂ζf t (r) 2 ∂Dpq

ρ(r) ∂ζt0 (r) 2 ∂Dpq

−

−

0 ρ(r) ∂ζf t (r) 2 ∂Dpq

ζt0 (r) ∂ρ(r) 2 ∂Dpq

−

ζf0 t (r) ∂ρ(r) 2 ∂Dpq

    0   1 ∂ρ (r)

R(r) < R0 R0 ≤ R(r) ≤ R1 R(r) > R1

2 ∂Dpq

29


(94)


Page 30 of 42

where we employ the following intermediate derivatives: ∂ρ(r) ∂Dpq ∂ρ0 (r) ∂Dpq ∂ζt (r) ∂Dpq ∂R(r) ∂Dpq ∂ζt0 (r) ∂Dpq ∂R0 (r) ∂Dpq ∂ζf t (r) ∂Dpq ∂ζf0 t (r) ∂Dpq

= φp (r)φq (r)

(95)

= φ0p (r)φq (r) + φp (r)φ0q (r)

(96)

1 ∂R(r) 2ζt (r) ∂Dpq 8Π(r) =− φp (r)φq (r) [ρ(r)]3 R0 (r) ∂ζt (r) 1 ∂R0 (r) = − 2[ζt (r)]2 ∂Dpq 2ζt (r) ∂Dpq 0 8Π(r) ∂ρ0 (r) 4Π (r) 24Π(r) ∂ρ(r) + = − [ρ(r)]3 [ρ(r)]4 ∂Dpq [ρ(r)]3 ∂Dpq =−

(97) (98) (99) (100)

∂R(r) (101) ∂Dpq ∂R(r) = 20A(R(r) − R1 )3 + 12B(R(r) − R1 )2 + 6C(R(r) − R1 ) R0 (r) ∂Dpq 0 ∂R (r) + [5A(R(r) − R1 )4 + 4B(R(r) − R1 )3 + 3C(R(r) − R1 )2 ] (102) ∂Dpq = [5A(R(r) − R1 )4 + 4B(R(r) − R1 )3 + 3C(R(r) − R1 )2 ]

30


Page 31 of 42 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 derivatives with respect to the two-body density matrix are given by    ρ(r) ∂ζt (r)   2 ∂dpqst   

∂ ρ˜α (r) = ρ(r) ∂ζf t (r) 2 ∂dpqst  ∂dpqst      0

R(r) < R0 R0 ≤ R(r) ≤ R1 R(r) > R1

   ∂ζt (r)  − ρ(r)  2 ∂dpqst   

R(r) < R0

∂ ρ˜β (r) = − ρ(r) ∂ζf t (r) 2 ∂dpqst  ∂dpqst      0    ρ(r) ∂ζt0 (r)  +  2 ∂dpqst   

∂ ρ˜0α (r) 0 = ρ(r) ∂ζf t (r) + 2 ∂dpqst  ∂dpqst      0

∂ ρ˜0β (r) ∂dpqst

=

R0 ≤ R(r) ≤ R1

∂ζf0 t (r) ∂dpqst

(104)

R(r) > R1 ρ0 (r) ∂ζt (r) 2 ∂dpqst ρ0 (r) ∂ζf t (r) 2 ∂dpqst

R(r) < R0 R0 ≤ R(r) ≤ R1

(105)

R(r) > R1

  ∂ζt0 (r)   − ρ(r) −  2 ∂dpqst    − ρ(r) 2       0

(103)

−

ρ0 (r) ∂ζt (r) 2 ∂dpqst ρ0 (r) ∂ζf t (r) 2 ∂dpqst

R(r) < R0 R0 ≤ R(r) ≤ R1 R(r) > R1

31


(106)


Page 32 of 42

where we have used the intermediates ∂ζt (r) ∂dpqst ∂R(r) ∂dpqst ∂Π(r) ∂dpqst ∂ζt0 (r) ∂dpqst ∂R0 (r) ∂dpqst ∂Π0 (r) ∂dpqst

1 ∂R(r) 2ζt (r) ∂dpqst 4 ∂Π(r) = [ρ(r)]2 ∂dpqst =−

(107) (108)

= φp (r)φq (r)φs (r)φt (r) 1 −R0 (r) ∂ζt (r) 1 ∂R0 (r) =− + 2 [ζt (r)]2 ∂dpqst ζt (r) ∂dpqst 0 0 4 ∂Π (r) 8ρ (r) ∂Π(r) = − 2 [ρ(r)] ∂dpqst [ρ(r)]3 ∂dpqst

(109) (110) (111)

= φ0p (r)φq (r)φs (r)φt (r) + φp (r)φ0q (r)φs (r)φt (r) + φp (r)φq (r)φ0s (r)φt (r) + φp (r)φq (r)φs (r)φ0t (r)

(112)

∂R(r) ∂ζf t (r) = [5A(R(r) − R1 )4 + 4B(R(r) − R1 )3 + 3C(R(r) − R1 )2 ] (113) ∂dpqst ∂dpqst ∂ζf0 t (r) ∂R(r) = 20A(R(r) − R1 )3 + 12B(R(r) − R1 )2 + 6C(R(r) − R1 ) R0 (r) ∂dpqst ∂dpqst 0 ∂R (r) + [5A(R(r) − R1 )4 + 4B(R(r) − R1 )3 + 3C(R(r) − R1 )2 ] (114) ∂dpqst (115)

References (1) Pulay, P. Ab initio calculation of force constants and equilibrium geometries in polyatomic molecules. Mol. Phys. 1969, 17, 197–204. (2) Helgaker, T.; Jørgensen, P. Analytical Calculation of Geometrical Derivatives in Molecular Electronic Structure Theory. Adv. Quantum Chem. 1988, 19, 183–245. (3) Kato, S.; Morokuma, K. Energy gradient in a multi-configurational SCF formalism and its application to geometry optimization of trimethylene diradicals. Chem. Phys. Lett. 1979, 65, 19–25.

32


Page 33 of 42 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) Rice, J. E.; Amos, R. D.; Handy, N. C.; Lee, T. J.; Schaefer, H. F. The analytic configuration interaction gradient method: Application to the cyclic and open isomers of the S3 molecule. J. Chem. Phys. 1986, 85, 963–968. (5) Helgaker, T.; Jørgensen, P. Configuration-interaction energy derivatives in a fully variational formulation. Theoret. Chim. Acta 1989, 75, 111–127. (6) Shepard, R.; Lischka, H.; Szalay, P. G.; Kovar, T.; Ernzerhof, M. A general multireference configuration interaction gradient program. J. Chem. Phys. 1992, 96, 2085–2098. (7) Yamaguchi, Y.; Goddard, J. D.; Osamura, Y.; Schaefer, H. F. A New Dimension to Quantum Chemistry: Analytic Derivative Methods in Ab Initio Molecular Electronic Structure Theory; Oxford University Press: New York, 1994. (8) St˚ alring, J.; Bernhardsson, A.; Lindh, R. Analytical gradients of a state average MCSCF state and a state average diagnostic. Mol. Phys. 2000, 99, 103–114. (9) Jagau, T.-C.; Prochnow, E.; Evangelista, F. A.; Gauss, J. Analytic gradients for Mukherjees multireference coupled-cluster method using two-configurational selfconsistent-field orbitals. J. Chem. Phys. 2010, 132, 144110. (10) Shiozaki, T.; Gyrffy, W.; Celani, P.; Werner, H.-J. Communication: Extended multistate complete active space second-order perturbation theory: Energy and nuclear gradients. J. Chem. Phys. 2011, 135, 081106. (11) Liu, F.; Kurashige, Y.; Yanai, T.; Morokuma, K. Multireference Ab Initio Density Matrix Renormalization Group (DMRG)-CASSCF and DMRG-CASPT2 Study on the Photochromic Ring Opening of Spiropyran. J. Chem. Theory Comput. 2013, 9, 4462– 4469. (12) Schutski, R.; Jimnez-Hoyos, C. A.; Scuseria, G. E. Analytic energy gradient for the projected HartreeFock method. J. Chem. Phys. 2014, 140, 204101. 33



(13) 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. (14) Park, J. W.; Shiozaki, T. Analytical Derivative Coupling for Multistate CASPT2 Theory. J. Chem. Theory Comput. 2017, 13, 2561–2570. (15) P.Celani,; Werner, H.-J. Analytical energy gradients for internally contracted secondorder multireference perturbation theory. J. Chem. Phys. 2003, 119, 5044–5057. (16) MacLeod, M. K.; Shiozaki, T. Automatic code generation enables nuclear gradient computations for fully internally contracted multireference theory. J. Chem. Phys. 2015, 142, 051103. (17) Roos, B. O. In Applications of Computational Chemistry: The First Forty Years; Dykstra, C. E., Frenking, G., Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005; pp 725–764. (18) Nakano, H.; Tsuneda, T.; Hirao, K. In Applications of Computational Chemistry: The First Forty Years; Dykstra, C. E., Frenking, G., Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005; pp 507–557. (19) Gordon, M. S.; Schmidt, M. W. In Applications of Computational Chemistry: The First Forty Years; Dykstra, C. E., Frenking, G., Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005; pp 1167–1189. (20) Szalay, P. G.; M¨ uller, T.; Gidofalvi, G.; Lischka, H.; Shepard, R. Multiconfiguration Self-Consistent Field and Multireference Configuration Interaction Methods and Applications. Chem. Rev. 2012, 112, 108–181. (21) Li Manni, G.; Carlson, R. K.; Luo, S.; Ma, D.; Olsen, J.; Truhlar, D. G.; Gagliardi, L. Multiconfiguration Pair-Density Functional Theory. J. Chem. Theory Comput. 2014, 10, 3669–3680. 34


Page 34 of 42

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


(22) Gagliardi, L.; Truhlar, D. G.; Li Manni, G.; Carlson, R. K.; Hoyer, C. E.; Bao, J. L. Multiconfiguration Pair-Density Functional Theory: A New Way To Treat Strongly Correlated Systems. Acc. Chem. Res. 2017, 50, 66–73. (23) Carlson, R.; Truhlar, D. G.; Gagliardi, L. On-Top Pair Density as a Measure of LeftRight Correlation in Bond Breaking. J. Phys. Chem. A 2017, 121, 5540–5547. (24) Becke, A. D.; Savin, A.; Stoll, H. Extension of the local-spin-density exchangecorrelation approximation to multiplet states. Theoret. Chim. Acta 1995, 91, 147–156. (25) Moscardó, F.; San-Fabián, E. Density-functional formalism and the two-body problem. Phys. Rev. A 1991, 44, 1549–1553. (26) Perdew, J. P.; Savin, A.; Burke, K. Escaping the symmetry dilemma through a pairdensity interpretation of spin-density functional theory. Phys. Rev. A 1995, 51, 4531– 4541. (27) Gusarov, S.; Malmqvist, P.-˚ A.; Lindh, R. Using on-top pair density for construction of correlation functionals for multideterminant wave functions. Mol. Phys. 2004, 102, 2207–2216. (28) Sand, A. M.; Truhlar, D. G.; Gagliardi, L. Efficient algorithm for multiconfiguration pair-density functional theory with application to the heterolytic dissociation energy of ferrocene. J. Chem. Phys. 2017, 146, 034101. (29) Andersson, K.; Malmqvist, P. ˚ A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. Secondorder perturbation theory with a CASSCF reference function. J. Phys. Chem. 1990, 94, 5483–5488. (30) Hirao, K. Multireference Møller-Plesset method. Chem. Phys. Lett. 1992, 190, 374–380. (31) Hirao, K. Multireference Møller-Plesset perturbation theory for high-spin open-shell systems. Chem. Phys. Lett. 1992, 196, 397–403. 35



(32) Hirao, K. State-specific multireference Møller-Plesset perturbation treatment for singlet and triplet excited states, ionized states and electron attached states of H2O. Chem. Phys. Lett. 1993, 201, 59–66. (33) Lischka, H.; Shepard, R.; Brown, F. B.; Shavitt, I. New implementation of the graphical unitary group approach for multireference direct configuration interaction calculations. Int. J. Quantum Chem. 1981, 20, 91–100. (34) Ghosh, S.; Sonnenberger, A. L.; Hoyer, C. E.; Truhlar, D. G.; Gagliardi, L. Multiconfiguration Pair-Density Functional Theory Outperforms KohnSham Density Functional Theory and Multireference Perturbation Theory for Ground-State and Excited-State Charge Transfer. 2015, 11, 3643–3649. (35) Hoyer, C. E.; Gagliardi, L.; Truhlar, D. G. Multiconfiguration Pair-Density Functional Theory Spectral Calculations Are Stable to Adding Diffuse Basis Functions. 2015, 6, 4184–4188. (36) Hoyer, C. E.; Ghosh, S.; Truhlar, D. G.; Gagliardi, L. Multiconfiguration Pair-Density Functional Theory Is as Accurate as CASPT2 for Electronic Excitation. 2016, 7, 586– 591. (37) Odoh, S. O.; Manni, G. L.; Carlson, R. K.; Truhlar, D. G.; Gagliardi, L. Separated-pair approximation and separated-pair pair-density functional theory. Chem. Sci. 2016, 7, 2399–2413. (38) Ghosh, S.; Cramer, C. J.; Truhlar, D. G.; Gagliardi, L. Generalized-active-space pairdensity functional theory: an efficient method to study large, strongly correlated, conjugated systems. Chem. Sci. 2017, 8, 2741–2750. (39) Wilbraham, L.; Verma, P.; Truhlar, D. G.; Gagliardi, L.; Ciofini, I. Multiconfiguration Pair-Density Functional Theory Predicts Spin-State Ordering in Iron Complexes with

36


Page 36 of 42

Page 37 of 42 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 Same Accuracy as Complete Active Space Second-Order Perturbation Theory at a Significantly Reduced Computational Cost. J. Phys. Chem. Lett. 2017, 8, 2026–2030. (40) Hellmann, H. Einf¨ uhrung in die Quantenchemie; Franz Deuticke: Leipzig, 1937. (41) Feynman, R. P. Forces in Molecules. Phys. Rev. 1939, 56, 340–343. (42) Kern, C. W.; Karplus, M. Analysis of Charge Distributions: Hydrogen Fluoride. J. Chem. Phys. 1964, 40, 1374–1389. (43) Nakatsuji, H. Common nature of the electron cloud of a system undergoing change in nuclear configuration. J. Am. Chem. Soc. 1974, 96, 24–30. (44) Hoffmann, M. R.; Fox, D. J.; Gaw, J. F.; Osamura, Y.; Yamaguchi, Y.; Grev, R. S.; Fitzgerald, G.; Schaefer, H. F.; Knowles, P. J.; Handy, N. C. Analytic energy second derivatives for general MCSCF wave functions. J. Chem. Phys. 1984, 80, 2660–2668. (45) Helgaker, T. U.; Almlöf, J.; Jensen, H. J. ˚ A.; Jørgensen, P. Molecular Hessians for largescale MCSCF wave functions. J. Chem. Phys. 1986, 84, 6266–6279. (46) Page, M.; Saxe, P.; Adams, G. F.; Lengsfield, B. H. Multireference CI gradients and MCSCF second derivatives. J. Chem. Phys. 1984, 81, 434–439. (47) Dudley, T. J.; Olson, R. M.; Schmidt, M. W.; Gordon, M. S. Parallel coupled perturbed CASSCF equations and analytic CASSCF second derivatives. J. Comput. Chem. 2006, 27, 352–362. (48) Bernhardsson, A.; Lindh, R.; Olsen, J.; F¨ ulscher, M. A direct implementation of the second-order derivatives of multiconfigurational SCF energies and an analysis of the preconditioning in the associated response equation. Mol. Phys. 1998, 96, 617–628. (49) Carlson, R. K.; Truhlar, D. G.; Gagliardi, L. Multiconfiguration Pair-Density Functional Theory: A Fully Translated Gradient Approximation and Its Performance for 37



Transition Metal Dimers and the Spectroscopy of Re2Cl82. J. Chem. Theory Comput. 2015, 11, 4077–4085. (50) Helgaker, T.; Jørgensen, P.; Olsen, J. Molecular Electronic-Structure Theory; Wiley: West Sussex, England, 2008. (51) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (52) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran 77: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, 1992. (53) Johnson, B. G.; Frisch, M. J. Analytic second derivatives of the gradient-corrected density functional energy. Effect of quadrature weight derivatives. Chem. Phys. Lett. 1993, 216, 133–140. (54) Piccardo, M.; Penocchio, E.; Puzzarini, C.; Biczysko, M.; Barone, V. SemiExperimental Equilibrium Structure Determinations by Employing B3LYP/SNSD Anharmonic Force Fields: Validation and Application to Semirigid Organic Molecules. J. Phys. Chem. A 2015, 119, 2058–2082. (55) Aquilante, F.; Autschbach, J.; Carlson, R. K.; Chibotaru, L. F.; Delcey, M. G.; Vico, L. D.; Galvn, I. F.; Ferr, N.; Frutos, L. M.; Gagliardi, L.; Garavelli, M.; Giussani, A.; Hoyer, C. E.; Manni, G. L.; Lischka, H.; Ma, D.; Malmqvist, P.; M¨ uller, T.; Nenov, A.; Olivucci, M.; Pedersen, T. B.; Peng, D.; Plasser, F.; Pritchard, B.; Reiher, M.; Rivalta, I.; Schapiro, I.; Segarra-Mart, J.; Stenrup, M.; Truhlar, D. G.; Ungur, L.; Valentini, A.; Vancoillie, S.; Veryazov, V.; Vysotskiy, V. P.; Weingart, O.; Zapata, F.; Lindh, R. Molcas 8: New capabilities for multiconfigurational quantum chemical calculations across the periodic table. J. Comput. Chem. 2016, 37, 506–541.

38


Page 38 of 42

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


(56) Dunning Jr., 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. (57) Kendall, R. A.; Dunning Jr., T. H.; Harrison, R. J. Electron affinities of the firstrow atoms revisited. Systematic basis sets and wave functions. J. Chem. Phys. 1992, 96, 6796–6806. (58) Tishchenko, O.; Zheng, J.; Truhlar, D. G. Multireference Model Chemistries for Thermochemical Kinetics. J. Chem. Theory Comput. 2008, 4, 1208–1219. (59) Calculate Root-mean-square deviation (RMSD) of Two Molecules Using Rotation, GitHub, http://github.com/charnley/rmsd, 1.2.5-19-g258e15a. (60) Herzberg, G. Molecular Spectra and Molecular Structure. III. Electronic Spectra and Electronic Structure of Polyatomic Molecules; Van Nostrand Reinhold: Princeton, New Jersey, 1966. (61) Hino, O.; Kinoshita, T.; Chan, G. K.-L.; Bartlett, R. J. Tailored coupled cluster singles and doubles method applied to calculations on molecular structure and harmonic vibrational frequencies of ozone. J. Chem. Phys. 2006, 124, 114311. (62) 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. (63) Cramer, C. J.; Nash, J. J.; Squires, R. R. A reinvestigation of singlet benzyne thermochemistry predicted by CASPT2, coupled-cluster and density functional calculations. Chem. Phys. Lett. 1997, 277, 311. (64) Lindh, R.; Bernhardsson, A.; Sch¨ ulz, M. Benzyne Thermochemistry: A Benchmark ab Initio Study. J. Phys. Chem. A 1999, 9913, 9913.

39



(65) de Visser, S. P.; Filatov, M.; Shaik, S. REKS calculations on ortho-, meta- and parabenzyne. Phys. Chem. Chem. Phys. 2000, 2, 5046–5048. (66) Crawford, T. D.; Kraka, E.; Stanton, J. F.; Cremer, D. Problematic p-benzyne: Orbital instabilities, biradical character, and broken symmetry. J. Chem. Phys. 2001, 114, 10638. (67) Slipchenko, L. V.; Krylov, A. I. Singlet-triplet gaps in diradicals by the spin-flip approach: A benchmark study. J. Chem. Phys. 2002, 117, 4694. (68) Smith, C. E.; Crawford, T. D. The structures of m-benzyne and tetrafluoro-m-benzyne. J. Chem. Phys. 2005, 122, 174309. (69) Wang, E. B.; Parish, C. A.; Lischka, H. An extended multireference study of the electronic states of para-benzyne. J. Chem. Phys. 2008, 129, 044306. (70) Li, X.; Paldus, J. Force field of para- and meta-benzyne diradicals: A multireference coupled-cluster study. J. Chem. Phys. 2010, 132, 114103. (71) Jagu, T.-C.; Prochnow, E.; Evangelista, F. A.; Gauss, J. Analytic gradients for Mukherjees multireference coupled-cluster method using two-configurational self-consistentfield orbitals. J. Chem. Phys. 2010, 132, 144110. (72) Schutski, R.; Jiméenez-Hoyos, C. A.; Scuseria, G. E. Analytic energy gradient for the projected HartreeFock method. J. Chem. Phys. 2010, 140, 114103. (73) Ray, S. S.; Ghosh, A.; Chattopadhyay, S. Taming the Electronic Structure of Diradicals through the Window of Computationally Cost Effective Multireference Perturbation Theory. J. Phys. Chem. A 2016, 120, 5897–5916. (74) Groner, P.; Kukolich, S. G. Equilibrium structure of gas phase o-benzyne. J. Molec. Struct. 2006, 780, 178–181.

40


Page 40 of 42

Page 41 of 42 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) Yu, H. S.; He, X.; Truhlar, D. G. MN15-L: A New Local Exchange-Correlation Functional for KohnSham Density Functional Theory with Broad Accuracy for Atoms, Molecules, and Solids. J. Chem. Theory Comput. 2016, 12, 1280–1293.

41



Graphical TOC Entry

42


Page 42 of 42

Analytic Gradients for Complete Active Space Pair-Density Functional

Recommend Documents