Active Space Selection Based on Natural Orbital Occupation Numbers

Subscriber access provided by UNIV OF WESTERN ONTARIO

Quantum Electronic Structure

Active Space Selection Based on Natural Orbital Occupation Numbers From N-Electron Valence Perturbation Theory Abhishek Khedkar, and Michael Roemelt J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b01293 • Publication Date (Web): 06 May 2019 Downloaded from http://pubs.acs.org on May 7, 2019

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

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

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

Journal of Chemical Theory and Computation

Active Space Selection Based on Natural Orbital Occupation Numbers From N-Electron Valence Perturbation Theory Abhishek Khedkar and Michael Roemelt∗ Lehrstuhl für Theoretische Chemie, Ruhr-Universität Bochum, D-44780 Bochum, Germany and Max-Planck Institut für Kohlenforschung, Kaiser-Wilhelm Platz 1, D-45470 Mülheim an der Ruhr, Germany E-mail: [email protected]

Abstract Efficient and robust approximations to the full configuration interaction (Full-CI) method such as the density matrix renormalization group (DMRG) and the Full-CI quantum Monte-Carlo (FCIQMC) algorithm allow for multi-configurational self-consistent field (MC-SCF) calculations of molecules with many strongly correlated electrons. This opens up the possibility to treat large and complex systems that were previously untractable but at the same time it calls for an efficient and reliable active space selection as the choice of how many electrons and orbitals enter the active space is critical for any multireference calculation. In this work we propose an Active Space Selection based on 1st order perturbation theory (ASS1ST) that follows a ’bottom-up’ strategy and utilizes a set of quasi-natural orbitals together with sensible thresholds for their occupation numbers. The required quasi-natural orbitals are generated by diagonalizing the ∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment

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

virtual and internal part of the one-electron reduced density matrix that is obtained from strongly contracted n-electron valence perturbation theory (SC-NEVPT) on top of a minimal active space calculation. Self consistent results can be obtained when the proposed selection scheme is applied iteratively. Initial applications on four chemically relevant benchmark systems indicate the capabilities of ASS1ST. Eventually, the strengths and limitations are critically discussed.

1

Introduction

For molecular systems with strongly correlated electrons, e.g. many transition metal containing compounds, a reliable electronic structure description can only be obtained by means of multiconfigurational methods. Most commonly used among these methods are the complete active space configuration interaction (CASCI) and complete active space self-consistent field (CASSCF) methods. Since CASSCF and CASCI take into account all possible configurations within the active orbital space they yield a qualitatively correct wavefunction as starting point for further treatment of dynamic electron correlation provided the active space is chosen adequately large. 1,2 Although well-established for years and very successful, a common concern for all multireference methods is their high computational cost, arising from both, the underlying CASSCF and the treatment of dynamic electron correlation. Recently, the latter problem has been addressed using local correlation schemes by the groups of Werner and Neese. 3,4 Nevertheless, as the number of active electrons and orbitals increases the underlying CASSCF or CASCI calculations quickly become unfeasible. 1,2 Consequently, the usage of CAS-based methods has been limited to cases with only few strongly correlated electrons. However, with the emergence of efficient and accurate approximations including the density matrix renormalization group (DMRG), Full-CI quantum Monte-Carlo, selective CI, Heat-Bath CI and others, treatment of large active spaces on the order of 20-50 active electrons and orbitals is possible without greater complication. 5–16 The ability to deal with large active spaces opens up the possibility to treat large and com2


Page 2 of 47

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


plex systems that were previously untractable but at the same time it calls for an efficient and reliable active space selection as the choice of how many electrons and orbitals enter the active space is critical for any multireference calculation. If it is chosen too small, important aspects of the chemical problem at hand might be neglected and not be correctly described while a too large active space causes large unnecessary computational costs and severe convergence problems during the orbital optimization. 17 However, it is by no means trivial to determine a priori how many electrons (and orbitals) are in fact strongly correlated and thus should be included in the active space. For simple systems like many mononuclear transition metal complexes the selection of the active space has traditionally been done on the basis of ’chemical intuition’. In such cases the recently introduced projection on the atomic valence active space (AVAS) by Knizia and coworkers can be used as an aid to identify the desired set of orbitals even for large orbital spaces. 18 Alternatively, the PiOS method by Sayfutyarova and coworkers offers a similar functionality for molecules with π-systems. 19 Nevertheless, while ’chemical intuition’ has been proven successful in many instances it will meet its limits for complicated systems, e.g. molecules with multiple transition metals in complicated coordination environments. Thus it is desirable to define a procedure that utilizes an objective measure to determine the number of strongly correlated electrons and orbitals. An early approach towards this goal uses unrestricted Hartree-Fock natural orbital occupation numbers to select the active space. 20 This procedure has been frequently used and performs well for many organic systems but it exhibits problems for systems in which a strongly occupied orbital has several dominant correlation partners and it is not well suited if excited states are to be calculated. 21 A more robust alternative is provided by an Ansatz based on the ’fractional occupation number weighted electron density’ (FOD) as introduced by Grimme and coworkers. 22,23 However, early studies show that this Ansatz in some cases fails to give fully reliable results for transition metal containing compounds. 23 Similarly, the semi-automated active space selection scheme by Truhlar and coworkers has been designed for excited state calculations of small organic molecules by means of multistate-CASPT2 and 3



multiconfigurational pair-density functional theory (MC-PDFT) only. 24 In contrast, Wouters and coworkers as well as Reiher and coworkers proposed procedures that have been shown to produce meaningful results for transition metal containing compounds as well. 25,26 Their common starting point is a set of orbitals obtained from a preceding SCF calculation. These orbitals are then used in a DMRG calculation with a small bond dimension (number of renormalized states) that features a large active space on the order of 50-100 orbitals around the HOMO-LUMO gap. In this way it is ensured that all chemically relevant orbitals are taken into account. From the converged DMRG calculation one can choose a set of strongly correlated orbitals on the basis of either the von Neumann entropy (Reiher) or the natural orbital occupation number (NOON, Wouters). On account of this sequence of two consecutive active space calculations one might label both procedures as ’top-down’- approaches. Although rather recently proposed the entropy based method by Reiher and coworkers has been successfully applied to a variety of systems ranging from organic compounds to transition metal complexes. Notably, it produced satisfactory results even for excited states and over the course of reaction profiles. 27 However, for large systems where hundreds of orbitals fall within a relative narrow energy region around the HOMO-LUMO gap all ’top-down’ approaches will inevitably face some problems. There is no way to include all possibly relevant orbitals in the first, approximate DMRG calculation without manual preselection which in turn leads to the original problem. In a recent work, this problem has been addressed by a scheme thats splits the full valence space of large molecules into smaller sets manageable size. 28 In this work we propose a ’bottom-up’ approach, the Active Space Selection based on 1st order perturbation theory (ASS1ST), as a feasible alternative for such cases. Its starting point is a CASSCF or DMRGSCF calculation with a minimal and thus tractable active space. The corresponding ground state wavefunction serves as unperturbed wavefunction Ψ(0) for a subsequent strongly contracted second order n-electron valence perturbation theory (SC-NEVPT2) calculation. 29,30 The resulting first order perturbed wavefunction then 4


Page 4 of 47

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


gives rise to non-vanishing elements in the internal/internal and external/external blocks of the reduced density matrix D. Separate diagonalization of these two blocks yields a set of quasi-natural internal and external orbitals. On the basis of their corresponding occupation numbers and a sensible threshold a new and more suitable active space is chosen. In that sense ASS1ST follows the spirit of the extremely successful pair-natural orbital approaches applied in many forms throughout quantum chemistry 3,4,31–36 that are ultimately based on Löwdin’s finding that natural orbitals provide the basis for efficient CI expansions. 37 Therefore, ASS1ST can be regarded as a generalization of the basic concept of using natural orbitals from second order Møller Plesset perturbation theory (MP2) 2 for genuine multiconfigurational cases where MP2 does not yield satisfactory results.A similar strategy as the presented one was previously used by Nooijen and coworkers to identify a reduced virtual orbital set in multireference equation-of-motion coupled cluster calculations. 38 In principle, the proposed ’bottom-up’ Ansatz is flawed by the fact that the obtained densities are only approximate as the underlying wavefunction is restricted to the first order perturbed wavefunction which might in some cases miss important configurations of the exact wavefunction. Nevertheless, we assume that through repeated application of ASS1ST until self-consistency is achieved all important configurations will be captured as observed in the Configuration Interaction by Perturbation with multiconfigurational zeroth-order wave functions Selected by Iterative process (CIPSI) and more recent variants thereof. 14,15,39,40 A set of test applications of ASS1ST on a set of four chemically relevant systems, the Hieber anion, Fe(II) porphyrin, a mixed-valence Mn dimer and deca-1,3,5,7,9-pentaene support this assumption as for all four compounds active spaces were obtained that constitute a reasonable trade-off between computational cost and accuracy, which is the ultimate goal of the procedure. To this end we would like to emphasize that ASS1ST does not aim to make the decision for an active space a ’black-box’ procedure but rather aims at providing a tool for theoretical chemists to put their chemically motivated choice on solid grounds and, where needed, refine it. In this regard, ASS1ST appears to be in the middle between completely unguided 5



Page 6 of 47

methods like Reiher’s entanglement entropy analysis and guided options like Knizia’s AVAS scheme. Both approaches to the active space selection problem, unguided and guided, have distinct advantages and disadvantages in different application scenarios. While the former set of methods allows for a full automation of the process and hence the rigorous definition of a ’model chemistry’ the latter set of techniques provides the possibility to adapt the choice to specific user-defined requirements, e.g. spin state energetics of weakly correlated complexes. The results obtained from ASS1ST can to a certain degree be influenced by a the user through the choice of the initial space and the selection thresholds. Of course, the final active space suggested by ASS1ST can furthermore be supplemented by any weakly correlated orbital if deemed necessary. In principle, one might even employ a combination of guided methods and ASS1ST to arrive at an active space that suits specific chemical requirements and is self-consistent in the sense of ASS1ST. However, such a combination is not explored within this work and remains to be examined in the future.

2

Theory

This section is structured as follows. First, N-Electron Valence Perturbation Theory and its strongly contracted variant are briefly introduced before we outline how the first order perturbed wavefunction is evaluated and used to construct quasi-natural orbitals. Eventually, the workflow for the ASS1ST scheme and a discussion of its computational cost is given.

2.1

Strongly Contracted NEVPT2

NEVPT2, introduced by Angeli and coworkers in 2001, is a second order perturbation theory 29,30 in which the zeroth order wavefunction is multireference in nature:

Ψ(0) m

=

CAS X

CI |Ii

I

6


(1)

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


(2)

(0) 0 ˆ CAS Ψ(0) PCAS HP m = E m Ψm (0)

where, Ψm is the m’th root of the underlying CASSCF problem, Em denotes the energy (0)

eigenvalue corresponding to Ψm , and PCAS is is the projector onto the CAS space. The space of perturber functions in NEVPT2 is defined as the first order interacting space (FOIS) which (0)

is generated by the action of the electronic Hamiltonian on the reference function Ψm . In the (k)

context of NEVPT2 the FOIS is divided into different classes of subspaces Sl

with k being

the number of electrons promoted to or from the active space (−2 < k < 2), and l denotes the fixed occupation pattern of the inactive orbitals. In the majority of implementations (including the present one) Dyall’s Hamiltonian is used as an approximation to the zero’th order Hamiltonian: 41 (3)

ˆ 0Dyall = H ˆc + H ˆv + C H X X ˆc = H i Eii + a Eaa

(4)

a

i

1X (tu|vw){Eut Ewv − δuv Ewt } 2 tuvw tu X X X C=2 hii + {2(ii|jj) + (ij|ij)} − 2 i .

ˆv = H

X

f u hef tu Et +

i

ij

(5) (6)

i

Here, we follow the convention that orbital labels i, j, k, l denote internal orbitals, t, u, v, w correspond to active orbitals and a, b, c, d are virtual or external orbital labels. The constant ˆ 0Dyall C is chosen to ensure that the zeroth order wavefunction is still an eigenfunction of H (0)

with eigenvalue Em . The operators Epq = a†qα aqα + a†qβ aqβ

(7)

are spin traced excitation operators for general orbital indices p, q while (ai|tu) are twoelectron integrals in Mulliken notation for spatial orbitals:

7



∞

Z

φa (r1 )φi (r1 )

(ai|tu) = −∞

Page 8 of 47

1 φt (r2 )φu (r2 )dr1 dr2 . r12

(8)

f ef f Eventually, a , i are orbital energies. The effective one-electron integrals hef pq and hpq

0

incorporate the interaction between active and internal orbitals and are defined as f hef pq = hpq +

X

(9)

2(ii|pq) + (pi|iq)

i 0

f f hef = hef pq pq −

X

(10)

(pt|tq).

t

(k)

In the strongly contracted variant of NEVPT2 (SC-NEVPT2) each subspace Sl

is repre-

sented by a single function (k)

Ψl

(k)

(11)

(0) = PS (k) HΨ(0) m = V l Ψm . l

According to their different excitation patterns, eight distinct classes of perturber operators Vl

(k)

can be defined: (0)

Vi,a =

X

(12a)

a {Eia Eut (ai|tu) + Eit Eua (au|ti)} + heff ai Ei

tu (0)

Vij,ab =γab γij {Eia Ejb (ia|jb) + Eib Eja (ib|ja)} X X a Va(−1) = Eua Evt (au|tv) + heff’ at Et tuv (−1)

(12b) (12c)

t

X

Vi,ab =γab

i≤j , a≤b

{Eia Etb (ai|bt) + Eib Eta (bi|at)}

a≤b

(12d)

t (+1)

Vi

=

X

Eiu Evt (ui|tv) +

tuv (+1) Vij,a

=γij

X

(12e)

t heff ti Ei

t

X

{Eja Eit (it|ja) + Eia Ejt (ia|jt)}

i≤j

(12f)

{Eiu Ejt (ui|tj)}

i≤j

(12g)

t (2)

Vij =γij

X tu

8


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


(−2)

Vab

=γab

X

{Eua Etb (au|bt)}

(12h)

a≤b

tu

where γij = 1 − 12 δij . It should be noted that the one and two-electron integrals act as weights in the summation over the active indices in order to include each contributing function weighted by their specific importance. The perturber functions corresponding to the perturber operators in equation (12) are in general not normalized with their squared norm (k)

Nl

being (k)

Nl

(k) (k) ˆ (k) † ˆ (k0 ) |Ψ(0) i . = hΨl |Ψl i = hΨ(0) m |(Vl ) Vl0 m

(13)

Using this norm it is straightforward to define normalized perturber functions as (k)

˜ (k) = qΨl . Ψ l (k) Nl

2.2

(14)

First Order Perturbed Wavefunction (k)

The strongly contracted perturber function Ψl

(k)

belonging to any Sl

interacts with the

(0)

zeroth order wavefunction Ψm only through Vlk . Therefore, the first order correction to the (k)

total wavefunction from the subspace Sl

becomes

D E E (k) (k) (0) (k) ˆ (0) ˜ ˜ E E Ψl Vl Ψm Ψl H Ψm ˜ (k) (k) (1) ˜ = |Ψm,kl i = Ψl Ψl (k) (k) ∆l ∆l D E q E (k) (k) (k) ˜ E E Ψ(k) Ψl Ψl N l l ˜ (k) ˜ (k) = = = Ψl Ψl (k) (k) (k) ∆l ∆l ∆l D

(k)

where ∆l

(15)

(k)

0 = Em − El . The complete first order perturbed wavefunction can then be given

as (1) (0) Ψm = Ψm +

E (k) Ψ X l (k)

k,l

9

∆l


(16)


Page 10 of 47

together with its norm X N (k) (1)

l N (1) = Ψ(1) . Ψ = 1 + m m (k) 2 (∆ ) l k,l

(17)

The reader is referred to the Supporting Information for the derivation of the latter.

2.3

Quasi-Natural Orbitals

In second quantization formalism the elements of the one-electron reduced density matrix for a normalized N -electron wavefunction Ψ, D are defined as

Dpq = Ψ Eqp Ψ

(18)

for general orbital indices p and q. D is a N orb ×N orb positive semi-definite hermitian matrix, where N orb is the orbital dimension. As the matrix D is hermitian it is possible to bring it to a diagonal form by means of a unitary transformation according to D = UDdiag U† . The columns of the transformation matrix U are composed of the eigenvectors {ur } of the density matrix. These eigenvectors constitute the natural orbitals and satisfy

Dur = nr ur .

(19)

The corresponding eigenvalues are usually interpreted as the natural orbitals occupation P numbers (NOONs) with nr ≥ nr+1 and T r(D) = T r(Ddiag ) = r nr = Nel . Löwdin has shown in 1955 that the natural orbitals originating from the diagonalization of the exact density are optimal to describe electron correlation in the sense that a configuration interaction procedure with excited determinants constructed from these orbitals leads to the most rapid convergence of the correlation energy possible. 37 In this work we make use of this finding by constructing natural orbitals from the first order perturbed density of SC-NEVPT2. Their occupation numbers are then used as a measure for strong correlation and as criterion for 10


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


the inclusion of natural orbitals in the active space. However, instead of diagonalizing the complete density matrix the internal/internal and external/external blocks (Dint and Dext ) of D are set up and diagonalized separately as indicated in Figure 1, leading to two sets of quasi-natural orbitals. It should be noted that since the active/active block of D is not calculated with the perturbed wavefunction, the sum over all diagonal elements does not match the number of electrons any more. Within the chosen approach the external/external block of D does not feature any contri-

Internal-Internal

External-External

Figure 1: The important blocks of the 1-electron reduced density matrix D (0)

butions from Ψm and is thus given by ext Dab =

=

1 (1) a (1) Ψm Eb Ψm N (1) D E (k) a (k0 ) E Ψ Ψ 1 X X l b l0 N (1)

(k)

(k0 )

∆l ∆l0 D E (0) ˆ (k) † a ˆ (k0 ) (0) Ψ ( V ) E V Ψ X X 0 m m b l l 1 = (1) . (k) (k0 ) N k,l k0 ,l0 ∆l ∆l0

(20a) (20b)

k,l k0 ,l0

11


(20c)


Page 12 of 47

(0)

In contrast, Ψm contributes to the internal-internal part of D but all terms of the form (0)

(k)

hΨm | Eji |Ψl i vanish. Hence,

(1) i (1) 1 2δ + Ψm Ej Ψm ij N (1) D E  (0) ˆ (k) † i ˆ (k0 ) (0) Ψ ( V ) E V Ψ X X 0 m m j l l 1 . = (1) 2δij + (k) (k0 ) N ∆l ∆l0 k,l k0 ,l0

int Dij =

(21a) (21b)

Importantly, since the number of holes and particles in the external and internal space is not altered by Eba and Eji only terms that include two perturber functions of the same class are non-vanishing in equations (20) and (21). For the evaluation of Dint and Dext according to equations (20) and (21) the energy denominators and perturber function norms from the original formulation of SC-NEVPT2 can be readily used. 30 The numerators can be derived and evaluated in a similar fashion. (−1)

Here, one example for the contributions from subspace Si,ab to Dext is presented to illustrate the procedure. The expression for Vˆi,ab from equation (12d) is reproduced here for clarity. (−1)

The perturber function Ψi,ab is given by (−1)

(−1)

|Ψi,ab i =Vˆi,ab |Ψ(0) ,a ≤ b m i X =γab Eia Etb (ia|tb) + Eib Eta (ib|ta) |Ψ(0) m i ,a ≤ b

(22)

t

where

γab = 1 −

12

1 · δab 2


(23)

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


It is advantageous to lift the restrictions on the two external labels for the evaluation of contributions to the external density. D (i,ab),ext

Def

=

XX

(0) Ψm

E ˆ (−1) † e ˆ (−1) (0) (Vi,ab ) Ef Vi0 ,a0 b0 Ψm (−1)

(−1)

∆i,ab ∆i0 ,a0 b0 XX

2(ei|ua)(f i|ta) (0) u (0) = Ψ m Et Ψ m ∆i,ea ∆i,f a i,a t,u i,ab i0 ,a0 b0

(i,ab),ext

Def

(24)

2(f t|ia)(eu|ia) (f i|ta)(eu|ia) (ei|ua)(f t|ia) − − + ∆i,ea ∆i,f a ∆i,ea ∆i,f a ∆i,ea ∆i,f a

where (−1)

(25)

∆i,ab = a + b − i − i,ab

In order to produce efficient computer code, the above terms have been reformulated in terms of matrices T of the large virtual orbital dimension: (i,ab)ext

D

=

XX i

=

2·

(0) Dtu

˜ iu T†it + T ˜ ui T†ti T

Tiu T†it

+2·

Tui T†ti

t,u

XX i

(0) Dtu

−

Tiu T†ti

−

Tui T†it

(26)

t,u

where ˜ it =2 · Tit − Tti T

(27)

(ia|tb) a + b − i + i,ab

(28)

u (0)

(0) Dtu = Ψ(0) m Et Ψ m .

(29)

Tit,ab = and

13



Page 14 of 47

The effective ionization potential

i,ab =

1 (−1)

Ni,ab

D

h i E ˆ (−1)† ˆ ˆ (−1) (0) Ψ(0) V H , V Ψ m v m i,ab i,ab

(30)

in the denominator of equation (28) is evaluated as described in the original work of Angeli and coworkers. 30 Analogous formulas for the other seven perturber classes are given in the Supporting Information.

2.4

Choosing Orbitals

The ASS1STed choice of an active space as proposed in this work includes several steps. First, the two sets of quasi-natural orbitals obtained by diagonalizing the two blocks of density matrices, DInt and DExt , are sorted according to their NOONs. Then the user defines a suitable NOON threshold TN OON (2 − TN OON for internal orbitals) to select the orbitals to be included in a larger, more adequate active space as indicated in Figures 2. This means each orbital with an occupation np that satisfies 2 − TN OON ≥ np ≥ TN OON is considered as strongly correlated and therefore suggested to be included in the updated active space. In principle, the entire procedure including the steps i) CASSCF or DMRGSCF calculation ii) SC-NEVPT2 density generation iii) natural orbital construction and iv) choice of orbitals can be repeated until self-consistency is achieved (see Figure 3 and Results section). The code required to conduct steps i) through iii) is implemented in the MOLBLOCK program which is being developed by our group and will soon be made available to the public. 17 In the first presented test case, a rather conservative threshold (TN OON ≈ 0.05) was used to exclusively include the orbitals expected to be responsible for strong or static correlation. A separate threshold for internals and externals can, in principle, be used if desired. In general the thresholds can and should be adjusted on grounds of a careful analysis of the investigated system and the sets of quasi-natural orbitals as demonstrated for the othertest cases. If the ASS1ST scheme is applied multiple times the intermediate CASSCF or DMRGSCF 14


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


Figure 2: Schematic representation of the distribution of NOONs together with the selection thresholds 2 − TN OON and TN OON . calculations can be performed with relaxed convergence criteria (∆E < 10−5 Eh and |g| < 0.001) and are expected to converge within a few (≈ 8 − 10) cycles as the quasi-natural orbitals would prove to be a near-ideal input. Importantly, as will be demonstrated below, the occupation numbers within the active space of these intermediate CASSCF calculations might suggest to reassign active orbitals to the internal or external space. In this way, the active space is not only allowed to grow but also to shrink. Towards the end of the iterative procedure, the comparison between different active spaces should give some indication of the stability of the results. Since both growing and shrinking of the active space in the course of the ASS1ST scheme is allowed one might encounter the situation where the procedure is stuck in an infinite loop between two or more active space choices. In such cases the user must either opt for one of the suggested active spaces or manually select a compromise based on ’chemical intuition’. A critical aspect of the presented procedure concerns the choice of an initial active space 15



Figure 3: Schematic representation of the steps during the ASS1ST procedure. since the outcome of the ASS1ST scheme is somewhat reliant upon the suitability of the initial active space used during the generation of the reference wavefunction. As a general rule of thumb, we advise potential users to employ a small but ’chemically reasonable’ initial active space. In cases with open shells this might correspond to only the singly occupied orbitals while for first-row transition metals one might include all orbitals with significant 3d character (see below). For organic compounds with large π systems a small number of occupied and unoccupied π-orbitals should provide a satisfactory starting point. In this regard, we would like to emphasize that the ASS1ST scheme is designed as a tool to aid or refine the active space choice of a potential user rather than making the decision a ’blackbox’ procedure. Hence, strict rules for the initial active space that can be applied regardless of the chemical system at hand cannot be made but the choice should be made on the 16


Page 16 of 47

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


basis of some a priori insight into the electronic structure under investigation, i.e. ’chemical intuition’. Although it has not been tested explicitly in this work, it might be a useful idea to utilize other, cheaper methods like Knizia’s AVAS scheme or Pulay’s NOON based criterion to identify a suitable initial active space. Importantly, if combined with conservative thresholds, different suitable choices for the initial active space will converge to the same active space as is demonstrated below.

2.5

Computational Costs

An important concern for virtually any multirefence method is the high computational cost associated with it. Therefore, a detailed discussion of the computational bottlenecks of the ASS1ST scheme is imperative. The first aspect to be discussed is the formation and transformation of two-electron integrals. Like the underlying NEVPT2 method, ASS1ST requires two-electron integrals with up to two virtual indices. Accordingly, with increasing size of the orbital space the formation and transformation of these orbitals becomes computationally expensive. A full four-index driven algorithm scales as O(N 5 ) and quickly becomes unfeasible. However, the resolution-of-identity approximation implemented in our MOLBLOCK code drastically reduces the computational cost for the integral evaluation as only threeindex quantities are generated and transformed. 42–45 With the help of the highly optimized LibInt library of Valeev and coworkers a highly efficient code is obtained that can readily handle at least 2000 basis functions. 46 Of course, the formation of one-electron integrals comes at much smaller cost than the two-electron integrals and can thus be neglected in this discussion. The second computational bottleneck arises from the generation of active space density matrices. Most notably the four-electron terms required the computation of energy denom(1)

inators of Vi

(−1)

and Va

become extremely expensive with increasing active space size as

8 they formally scale as O(Nactive ). 29,30 A discussion of the computational costs of all active

space density matrices is given in reference 47 that describes their implementation for DMRG 17



based wavefunctions in the BLOCK code. Without going into detail we simply state that when active spaces on the order of 15-20 orbitals were used for the two presented cases in this work, the generation of the aforementioned four-electron active space densities dominated the computational cost by taking between 75-80 % of the CPU time. Owing to this steep increase ASS1ST in its current form is limited to active space sizes on the order of 20-30 orbitals putting many interesting systems out of reach that are readily tractable by the DMRG and FCIQMC. 47 However, we envision a replacement of the exact energy denominators by a mere sum of orbital energies as it occurs in CASPT2 which would relief ASS1ST of the need to evaluate the four-electron densities altogether. 48,49 Work along these lines is currently being conducted in our group. With the required integrals and active space densities at hand the calculation of contribu(0)

(−1)

(−2)

tions from Vij,ab , Vi,ab and Vab

to Dext are most expensive since they all involve summa-

tions over three external indices (see Supporting Information). More precisely, the evalua2 3 2 3 tion of the corresponding terms scale as O(Ninternal Nexternal ), O(Ninternal Nactive Nexternal ) and 4 3 O(Nactive Nexternal ), respectively. Accordingly, the relative cost of these terms depends on the

relative size of the internal and active orbital spaces. To this end it must be noted that the evaluation of density matrix contributions is by construction more expensive than the energy evaluation within regular SC-NEVPT2. As such, it significantly adds to the computational costs of any computational investigation in contrast to other methods, e.g. Pulay’s approach based on unrestricted natural orbitals which comes at almost no cost at all. 20,21 Nevertheless, ASS1ST is feasible for medium sized problems with many strongly correlated electrons as is demonstrated below. Moreover, as a correct description of static electron correlation effects usually does not require large basis sets one could resort to using small basis sets for large systems without losing any essential information.

18


Page 18 of 47

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


2.6

Computational Details

All presented CASSCF, DMRGSCF and SC-NEVPT2 calculations, including the generation of active space density matrices and quasi-natural orbitals were performed using the MOLBLOCK program which at this stage is an in-house code but will be made available to the public in 2019. 17 Importantly, the MOLBLOCK code utilizes the extremely efficient BLOCK DMRG-CI code by Chan and coworkers. 5,47,50,51 As no dedicated CAS-CI code is yet implemented in MOLBLOCK the CASSCF calculations reported in this work in fact correspond to DMRGSCF calculations with a bond dimension that renders the DMRG-CI step exact. All reported multireference calculations employed the def2-TZVP basis set. The resolution-of-identity approximation was applied together with the def2-TZVP/C basis set to reduce the computational cost and storage requirements during the evaluation of the required two-electron integrals. 42–45,52 If not stated otherwise,all geometries underlying the aforementioned correlated calculations were optimized on the DFT level of theory using the ORCA program package in its version 4.0.1. 53,54 The hybrid density functional B3LYP was employed together with the def2TZVP basis set for the geometry optimizations. 55–57 In the following the conventional notation (n, m) will be adopted to denote an active space with n electrons in m orbitals, e.g. CASSCF(n, m).

3 3.1

Applications Hieber Anion

The first test case for the ASS1ST scheme presented in this work is the complex anion [Fe(CO)3 (NO)]− (1, Figure 4) also referred to as Hieber anion which exhibits catalytic activity in various organic reactions and has been extensively studied both experimentally and theoretically. 58–61 Complex 1 is isoelectronic with the anion of Collman’s reagent [Fe(CO)4 ]2−

19



and it is known to play a crucial role in highly-unusual nitrosyl-ligand based oxidation. 60 Furthermore, analysis of its singlet ground-state wavefunction indicates a profound multiconfigurational character calling for the use of multiference methods to describe its electronic structure. A previous study by Plietker and coworkers applied an active space of 14 electrons in 9 orbitals that is supposed to include all Fe 3d orbitals as well as NO π and π∗ orbitals. In contrast, a more recent study by Knizia and coworkers suggested to use either a (10,8) or (14,10) active space including orbitals of Fe 3d and optionally N 2p character which is further testimony to the ambiguity in the active space choice. The following describes how ASS1ST may be used to determine a sensible active space for this intriguing example. All presented calculations were performed using a geometry that was optimized for the singlet ground state at the aforementioned level of DFT.

Figure 4: The Hieber anion denoted complex 1 in this work. The starting point of the procedure is an initial CASSCF(8,5) calculation inspired by the notion that 1 features an Fe center with a formal 3d8 configuration. As depicted in Figure 5 the subsequent ASS1ST(8,5) calculation together with a threshold of TN OON = 0.05 suggests to incorporate four additional internal and external orbitals, respectively, resulting in a (16,13) active space. Before engaging in a further run of ASS1ST this choice is reviewed by performing a CASSCF(16,13) calculation, which reveals that three active orbitals now feature occupation numbers close to 2.0 (see upper right side of Figure 6). Accordingly, they were reassigned to the internal space leading to a (10,10) active space. Visual inspection together with an analysis of orbital compositions reveals that the three omitted natural 20


Page 20 of 47

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


Figure 5: Quasi-natural orbitals and their occupation numbers from a ASS1ST(8,5) calculation. A threshold of TN OON = 0.05 (dashed line) is employed resulting in the inclusion of four additional internal (top) and four additional external (bottom) orbitals. orbitals have major contributions from Fe 3p orbitals thus underlining the importance of a critical assessment of the suggested active orbitals by the user. A following ASS1ST(10,10) calculation rendered no further orbitals suitable for inclusion in the active space (see Figure 7). The same is found true by the ASS1ST(16,13) calculation (see Figure 6). In the following the final choice obtained from the described sequence of steps will be denoted (10,10)β . A second set of calculations with a CASSCF(10,8) starting point was carried out to assess the dependence of the outcome of the procedure on the initial choice of active space. As illustrated in Figure S12 the ASS1ST(10,8) calculation together with the same threshold as above (TN OON = 0.05) suggests to include two internal and two external quasi-natural

21



Figure 6: Quasi-natural orbitals and their occupation numbers from a ASS1ST(16,13) calculation. When a threshold of TN OON = 0.05 (dashed line) is employed no further orbitals will be included in the active space.

22


Page 22 of 47

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


Figure 7: Quasi-natural orbitals and their occupation numbers from a ASS1ST(10,10)α calculation. When a threshold of TN OON = 0.05 (dashed line) is employed no further orbitals will be included in the active space.

23



orbitals into the new active space thus leading to a (14,12) choice. While ASS1ST(14,12) does not recommend to include any further internal or external orbitals, the natural orbital occupancies within the active space (see Figure S13) suggest to discard two almost doubly occupied active orbitals. The resulting (10,10)α active space is virtually indistinguishable from the previously generated (10,10)β one upon visual inspection of the orbitals, composition of orbitals and their natural occupancies (Figure S1 and Figure S2). Analogous procedures with starting points (8,5) and (10,8) were carried out for the first triplet state of complex 1. In both cases the initial choice evolved to an active space of (14,12).1 As observed above two orbitals in that (14,12) space exhibit occupation numbers close to 2.0 (Figure S15). Reassignment of these to the internal space again resulted in a (10,10)γ active space (Figure S14). In order to evaluate the adequacy of the different active space choices, the singlet-triplet gap was calculated on the CASSCF and CASSCF + SC-NEVPT2 level of theory. Unfortunately there is no experimental reference data available for this system but vertical excitation energies on the highly accurate MRCI + Q level of theory were reported previously. 61 Table 1 summarizes the calculated vertical transition energies in comparison to the reference value. Unsurprisingly, the closest agreement with the reference value is obtained with the CASSCF + SC-NEVPT2(14,12) calculation (∆E = 0.01 eV) as it features the largest active space. However, almost equally excellent agreement is achieved from the smaller (10,10) active space which lies within 0.11 eV of the reference value. Interestingly, the marginally smaller (10,8) active space already gives a significantly larger deviation (1.1 eV) on the CASSCF + SC-NEVPT2 level of theory while the mere CASSCF(10,8) value is even closer to the reference value the CASSCF(10,10) value.As expected, the results from the considerably smaller (8,5) active space also deviate significantly from the reference value on the order of 0.5 eV and 0.9 eV with and without SC-NEVPT2, respectively. Although the accuracy of the calculated singlet triplet gap is certainly not an ideal measure for the adequacy of the 1

It must be noted that for ASS1ST(10,8) in the triplet case one orbital was included with an occupation number of 1.96.

24


Page 24 of 47

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


active space to describe chemical properties we take this result as a confirmation that the suggested (10,10) active space yields a suitable basis to describe the electronic structure of complex 1. Potential differences between the different choices for other physical or chemical properties will not be investigated as it is not the main subject of this work. Table 1: Calculated vertical excitation energies from the singlet ground state to the lowest triplet state in eV. Active Space CASSCF CASSCF+SC-NEVPT2 (8, 5) 1.905 2.924 (10, 8) 2.694 3.485 (10, 10)α 2.863 2.435 (14, 12) 2.953 2.330 MRCI + Q 2.320

3.2

Fe(II) porphyrin

Figure 8: Fe(II) porphyrin denoted complex 2 in this work. The well studied Fe(II) porphyrin complex 2 (Figure 8) has been chosen as second test system for ASS1ST. Owing to its unique electronic structure featuring multiple low-lying electronic states, complex 2 has been the subject of a number of experimental and theoretical investigations targeting the nature of its electronic ground state. 8,62–73 While Experimental studies on multiple substituted derivatives such as tetraphenylporphine (FeTPP) and Fe octaethylporphine (FeOEP) agree on a triplet ground state, quantum chemical calculations of complex 2 obtained triplet and quintet ground states depending on the methods used. 25



However, the most recent and accurate calculations seem to agree on a triplet ground state as well but with extremely close-lying triplet and quintet excited states. 72,73 The presented calculations on complex 2 in this work are based on geometries that were optimized on the DFT level for each spin state separately. We first focused on selecting a suitable active space for the high-spin (HS) quintet spin state. The straightforward minimal active space that includes six electrons in five orbitals showed severe convergence issues during the orbital optimization and therefore a starting point of ten electrons in nine orbitals was chosen (see Figure S3). This proved to be the smallest possible option with stable convergence behavior for both spin states under consideration.2 An ASS1ST(10,9) calculation with the same thresholds as above does not recommend to select any further orbitals into the active space. However, inspection of the natural orbitals and their occupation number distribution in Figure 9 reveals several orbitals with ligand π character that are strongly correlated but do not quite cross the threshold. Therefore, the threshold was slightly relaxed to a value of TN OON = 0.03 which in turn leads to selection of six additional internal and two additional external orbitals and hence a (22,17) active space (see Figure S6 and Figure S7). The fact that the corresponding DMRGSCF calculation with 17 orbitals converged within three cycles confirms that the quasi-natural orbitals produced in the ASS1ST run are indeed near-ideal input for subsequent CASSCF or DMRGSCF calculations. Reassignment of four active orbitals with an occupation close to 2.0 (see Figure 10) eventually yields the (14,13) active space visualized in Figure 11. As depicted in Figure 10 the corresponding ASS1ST(14,13) calculation does not recommend to select any further orbitals for the active space. At this point it should be noted that the ASS1ST(14,13) calculations described here and below were not conducted using a CASSCF reference but with a DMRGSCF reference with m = 750. An analogous sequence of calculations was conducted for the triplet state. The initial ASS1ST(10,9) calculation with a threshold of TN OON = 0.03 suggested to include seven 2 Although one active orbital of the (10,9) space features an occupation of 1.99 it was not reassigned to the internal space to retain comparability with the procedure for the triplet state.

26


Page 26 of 47

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


Figure 9: Quasi-natural orbitals and their occupation numbers from a ASS1ST(10,9) calculation. A threshold of TN OON = 0.05 leads to no further orbitals to be included while a threshold of TN OON = 0.03 leads to selection of six additional internal (top) and two additional external (bottom) orbitals

Figure 10: Quasi-natural orbitals and their occupation numbers from a ASS1ST(14,13) calculation with m = 750. No further orbitals are suggested for the active space with a threshold of TN OON = 0.04.

27



Figure 11: The final CASSCF(14,13) active natural orbitals and their occupancies for the quintet state internal and two external orbitals (see Figure S16). Augmenting the active space with these orbitals lead to a (24,18) space which after removal of five natural orbitals with occupations close to 2.0 lead to a (14,13) active space. Like for the quintet state ASS1ST(14,13) did not render any further orbitals suitable for inclusion into the active space (see Figure S17). The fact that the ASS1ST procedure arrives at the same active space choice for both spin states is fortuitous since it allows to make a meaningful comparison of their energies at the CASSCF and CASSCF+SC-NEVPT2 level of theory. Before proceeding to compare calculated energy gaps between the triplet and quintet states we would like to briefly discuss the character of some of the chosen active orbitals. Already in 1961, Gouterman identified four ligand based π orbitals as central unit to achieve a qualitative understanding of the observed absorption and emission spectra of typical D4h symmetrical porphyrins. 74 The recent study of Li Manni and coworkers confirms the importance of the ’Gouterman orbitals’ also for quantitative calculations. However, the ’Gouterman orbitals’ could only be identified in the active space of large-scale FCIQMC-SCF calculations when an all-valence active space of 32 electrons in 34 orbitals was employed. In contrast, the much

28


Page 28 of 47

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


smaller (14,13) active space obtained by ASS1ST for both spin states already incorporates four orbitals with occupation numbers significantly off 0.0 and 2.0, respectively (0.12 and 1.89 in Figure 11 and Figure S5). To this end we will briefly discuss the size of calculated triplet-quintet gaps using different active spaces. Table 2 summarizes the calculated triplet-quintet gap for the (10,9) and (14,13) active spaces on the CASSCF and CASSCF+SC-NEVPT2 level of theory in kcal mol−1 and compares it to recent calculations by Pierloot and coworkers and Li Manni and coworkers. 72,73,75 Note, a negative number in Table 2 corresponds to a quintet ground state. The obtained results show that, as expected, the larger active space corresponds to a smaller gap between the two states. Furthermore, the results reiterate the critical role of dynamic electron correlation for a realistic description of the triplet-quintet gap. However, even when SC-NEVPT2 is invoked to incorporate dynamic electron correlation, the gap between the different spin states appears to be slightly overestimated. Whether the best strategy to obtain more accurate results corresponds to enlarging the active space further or to employ more sophisticated methods to describe dynamic electron correlation on top the here suggested active spaces remains open to debate. Regardless of this discussion we believe that the active space choice provided by ASS1ST encompasses at least the majority of strongly correlated orbitals including the ’Gouterman orbitals’. Table 2: Calculated triplet-quintet gaps form complex 2 from various levels of theory in kcal mol−1 Active Space (10, 9) (14, 13) CASPT2(8,11) 75 FCIQMC(32,34) 72 CASPT2/CC(8,11) 73

CASSCF CASSCF+SC-NEVPT2 -30.44 -10.94 -34.88 -9.98 -4.9 3.1 0.2

29



Figure 12: The mixed valence Mn(III/IV) dimer denoted complex 3 in this work.

3.3

Mixed-Valence Mn Dimer

While the two test cases presented above feature a single transition metal center, complex 3 shown in Figure 12 is a mixed-valence bis-µ-oxo/acetato Mn(III/IV) dimer with a 1.4.7trimethyl-1,4,7-triazacyclononane and acetates as ligand. First reported by Wieghardt and coworkers in 1992, it is relevant to Mn based single molecular magnets and at the same time mimics some of the structural and electronic features of the oxygen evolving complex (OEC) of photosystem II (PSII). 76–78 Owing to its formal d4 configuration on one of the Mn centers, it is subject to strong axial Jahn-Teller distortion leading to an approximate squarepyramidal coordination geometry (see Figure 12). In the present context we are mostly interested in the magnetic exchange coupling between the local spins associated with the two metal centers since a physically sound description of the multiple involved spin states requires multireference methods and numerically accurate results are extremely difficult to obtain. 79 A recent study by Pantazis and coworkers used 3 as a test case to demonstrate the capabilities of DMRG based methods to quantitatively predict magnetic exchange coupling constants in model systems that feature structural motifs from the OEC of PSII. 17 In the course of that study multiple sets of chemically related orbitals from the core motif (e.g. Mn 3d’s, O 2p’s etc) and various combinations thereof have been employed in the active space to establish an order of importance for the calculation of the magnetic exchange coupling 30


Page 30 of 47

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


constant J. Importantly, it has been found that in order to achieve numerical agreement with the experimental result (J = −90 cm−1 ) the metal 3d orbitals and the 2p orbitals of the two oxo-bridges need to be included in the active space while it is sufficient to account for interactions with and between the remaining orbital space by means of second order perturbation theory. Here we apply ASS1ST to investigate its capabilities to capture the subtle interactions required to accurately describe magnetic exchange coupling in 3 and comparable model complexes. All test calculations reported in the following use a structure reported earlier which was optimized on the TPSS/def2-TZVP level of theory. 17,57,80,81 We chose the spin doublet ground state together with the minimal active space that features

Figure 13: Quasi-natural orbitals and their occupation numbers from a ASS1ST(7,7) calculation. Two internal and two virtual orbitals appear to be subject to strong correlation. seven unpaired electrons in seven Mn d orbitals as starting point for the ASS1ST scheme in the case of 3. The corresponding occupation number distribution depicted in Figure 13 indicates two internal and two virtual orbitals being strongly correlated whose inclusion leads to an (11,11) active space. A subsequent ASS1ST(11,11) calculation renders this choice self-consistent (see Figure S18). Interestingly, the pairs of additional internal and virtual orbitals comprise the bonding and antibonding combinations of Mn(III) d orbitals 31



and O based p orbitals (see Figure S10). Although this result appears chemically sensible and somewhat agrees with the pronounced importance of the oxo-bridges established by Pantazis and coworkers the present active space choice deviates considerably from their proposition for the calculation of the exchange coupling constant. 17 It is thus unsurprising that the relative spin state energies calculated at the state-averaged CASSCF(11,11) and SC-NEVPT2(11,11) levels of theory (see Table 3) exhibit the wrong order, i.e. the high spin state is most stable while the low-spin state is highest in energy. More precisely, the spin state energies in Table 3 correspond to almost identical, positive magnetic exchange coupling constants of J = 4.8 cm−1 (CASSCF) and J = 4.5 cm−1 (SC-NEVPT2). In summary, while the active space suggested by ASS1ST for 3 is from our point of view chemically sensible as it covers all singly occupied orbitals and all orbitals involved in the Mn(III)-oxo bond, it is not sufficient as a basis to obtain accurate exchange coupling constants. Table 3: Calculated energies of the different spin states of 3 relative to the octet ground obtained from state-averaged CASSCF(11,11) and CASSCF(11,11) + SC-NEVPT2. All results are given in cm−1 . Spin Multiplicity 6 4 2

3.4

CASSCF CASSCF+SC-NEVPT2 32.5 30.1 56.4 52.5 71.1 66.5

Deca-1,3,5,7,9-pentaene

Figure 14: The deca-1,3,5,7,9-penatene molecule denoted 4 in this work. Eventually, the performance of ASS1ST for extended conjugated π systems is tested through application to the deca-1,3,5,7,9-pentaene molecule 4 shown in Figure 14. Previous experimental and theoretical works by Hudson and coworkers reported and explained 32


Page 32 of 47

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


a dark electronic state lying at lower energies than the optically allowed HOMO-LUMO transition. 82,83 The theoretical work by Karplus and coworkers demonstrated that a multiconfigurational treatment was necessary to properly account for this unusual low energy state. 84 Interestingly, the electronic structure of this and analogous low lying states in structurally similar compounds lies at the heart of energy transport in systems ranging from conjugated organic semiconductors to the biological centers of light harvesting and vision. 85 Here, ASS1ST is applied to determine a suitable active space for 4 based on the density

Figure 15: Quasi-natural orbitals and their occupation numbers from a ASS1ST(6,6) calculation. Two internal and two virtual orbitals appear to be subject to strong correlation. of only the electronic ground state. An extension of the ASS1ST scheme towards simultaneously taking into account multiple electronic states is envisaged for the near future and will be tackled soon in our group. Following the recommendation from the previous section an initial active space with six electrons in six π orbitals was chosen. The corresponding quasi-natural orbital occupation number distribution in Figure 15 reveals the existence of two pairs of rather strongly correlated internal and external orbitals. As expected, these orbitals coincide with the remaining π orbitals and their inclusion leads to the all-π (10,10) active space previously employed by Hirao and coworkers (see Figure S11). 86 Subsequent 33



application of ASS1ST shows that no further alteration of the active space is recommended (see Figure S19). If the active space is nevertheless extended by two internal and two virtual orbitals the corresponding CASSCF(14,14) occupation numbers of the additional orbitals are >1,98 and

Active Space Selection Based on Natural Orbital Occupation Numbers

Recommend Documents