A Nonorthogonal State-Interaction Approach for Matrix Product State

The key element is the transformation of the MPS wave functions of different ... Ab Initio Study of Circular Dichroism and Circularly Polarized Lumine...
0 downloads 0 Views 908KB Size
Subscriber access provided by University of Otago Library

Article

A nonorthogonal state-interaction approach for matrix product state wave functions Stefan Knecht, Sebastian Keller, Jochen Autschbach, and Markus Reiher J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00889 • Publication Date (Web): 07 Nov 2016 Downloaded from http://pubs.acs.org on November 8, 2016

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

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

Page 1 of 37

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

Journal of Chemical Theory and Computation

A nonorthogonal state-interaction approach for matrix product state wave functions Stefan Knecht,1, a) Sebastian Keller,1 Jochen Autschbach,2 and Markus Reiher1 1)

ETH Z¨ urich, Laboratorium f¨ ur Physikalische Chemie, Vladimir-Prelog-Weg 2,

8093 Z¨ urich, Switzerland 2)

University at Buffalo, State University of New York, Department of Chemistry,

New York 14260-3000, United States of America We present a state-interaction approach for matrix product state (MPS) wave functions in a nonorthogonal molecular orbital basis. Our approach allows us to calculate for example transition and spin-orbit coupling matrix elements between arbitrary electronic states provided that they share the same one-electron basis functions and active orbital space, respectively. The key element is the transformation of the MPS wave functions of different states from a nonorthogonal to a biorthonormal molecular orbital basis representation exploiting a sequence of non-unitary transformations following a proposal by Malmqvist (Int. J. Quantum Chem. 30, 479 (1986)). This is well-known for traditional wave-function parametrizations but has not yet been exploited for MPS wave functions.

a)

Electronic mail: [email protected]

ACS Paragon 1 Plus Environment

Journal of Chemical Theory and Computation

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

I.

INTRODUCTION The calculation of electronic and vibronic transition matrix elements between electronic

states of the same or different spin and/or spatial symmetry is a ubiquitous task in the modeling of photochemical and photophysical processes. Prime examples include the modeling of non-adiabatic dynamics processes1,2 as well as light-induced excited spin-state trapping phenomena3,4 . The theoretical description of these processes builds upon the calculation of intersystem crossing rates5 which, besides electronic and vibronic coupling elements, requires the evaluation of spin-orbit (SO) coupling (SOC) matrix elements. Similarly, calculating magnetic properties6 such as molecular g-factors and electron-nucleus hyperfine coupling, which are central parameters in electron paramagnetic resonance (EPR) spectroscopy, requires spin-orbit coupled wave functions7,8 . To this end, correlated two- and four-component ab initio wave function9–12 , and density functional theory approaches13–15 . In the present study we focus on EPR g-tensors for testing purposes, but it should be noted that the underlying novel method of gaining access to wave functions that include the effects from SOC has a vast range of applications. While it is possible to treat SOC variationally, a considerable number of two-step correlated wave function approaches for the calculation of molecular g-factors were developed over the past decades (see, for example, Refs. 9,16–23 and literature citations in these works and in Refs. 8–11,13–15). In these schemes, the calculation of a number of non- or scalar-relativistic many-particle spin-free states that are eigenfunctions of the spin-squared operator S 2 , is decoupled from a subsequent perturbative or variational mixing of the latter through the SO coupling operator to obtain SO coupled many-electron wave functions (e.g. by diagonalization “state-interaction”). It is straightforward to calculate properties such as g-factors in the basis of the eigenstates of the SO operator subsequently. Appealing features of the two-step approaches are that valuable insight into contributions of each (ground or excited) spin-state to the g-tensor is gained, and that the underlying wave function basis is spin-adapted. The price to pay is the need to calculate a sufficient number of spin-free states to interact, which can amount up to several hundred states to achieve convergence for (heavy-element containing) molecules where electron correlation and spin-orbit coupling contributions can be of similar order of magnitude. Open shell electronic structures are often governed by strong electron correlation effects. ACS Paragon 2 Plus Environment

Page 2 of 37

Page 3 of 37

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

Journal of Chemical Theory and Computation

In this context, multiconfigurational methods are the preferred methods of choice24,25 which typically split electron correlation into a static and a dynamic contribution. However, such a separation requires careful attention26 . A well-established approach to handle static correlation is the complete active space self-consistent field (CASSCF) ansatz27 which requires to select a tailored number of (partially occupied) active orbitals. The selection of active orbitals is a tedious procedure but can be automatized26,28 . Since the computational cost of traditional CASSCF scales exponentially with the number of active orbitals and electrons, tractable active orbital spaces are presently limited to about 18 electrons in 18 orbitals29 . These limitations can be overcome by resorting to the density matrix renormalization group (DMRG) approach30–33 in quantum chemistry34–44 which, in combination with a self-consistent-field orbital optimization ansatz (DMRG-SCF)45,46 , is capable of approximating CASSCF wave functions to chemical accuracy with merely a polynomial scaling. DMRG-SCF therefore allows to handle much larger active orbital spaces that can boldly surpass the CASSCF limit. To account in addition for spin-orbit coupling in a DMRG framework, variational SO approaches47 as well as two-step approaches22,23 based on spinfree DMRG wave functions have been reported recently. In this work, we present a generalized state-interaction approach for nonorthogonal spinfree matrix product state (MPS) wave functions which enables the evaluation of arbitrary one- and two-particle transition matrix elements as well as SO coupling matrix elements. Diagonalization of the SO Hamiltonian matrix, for instance, yields spin-orbit coupled wave functions as linear combinations of the uncoupled, spin-pure MPS states. The latter can (but do not have to) be obtained as results from one or several DMRG-SCF orbital optimization calculations. This allows for utmost flexibility in the individual DMRG-SCF steps as each state-specific or spin-specific state-averaged orbital optimization is given the possibility to reflect potential relative differences in open-shell occupancies. For example, transition metal as well as lanthanide and actinide complexes often exhibit different s- and d-occupations (transition metals) and s-, d- and f -occupations (lanthanides/actinides) in ground- and electronically excited electronic states of various spin symmetries. However, a set of wave functions that were optimized individually generally implies mutual nonorthogonality of the respective MO bases. Moreover, the MO bases may neither be orthogonal to each other nor non-interacting which can, e.g., strongly affect the calculation of transition moments between such electronic states48 . As first shown in a landmark paper ACS Paragon 3 Plus Environment

Journal of Chemical Theory and Computation

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

Page 4 of 37

by Malmqvist49 , an elegant approach for the calculation of matrix elements and transition density matrices is the transformation to a biorthonormal basis for the bra and ket orbital basis of the respective wave functions. A change of the MO basis to a biorthonormal orbital basis necessitates, however, not only to transform all one- and two-electron integrals but to also “counter-rotate” the configuration basis of the wave function. For configurationinteraction (CI)-type expansions, the rotations and counter-rotations can be achieved by a sequence of single-orbital transformations49,50 that require only one-electron operations. Subsequently, standard second-quantization algebra can be exploited for the evaluation of overlap matrix elements as well as arbitray one- and two-particle matrix elements between the states in the biorthonormal basis. In a recent work, Olsen51 exploited the potential of the biorthonormal approach further to devise an efficient algorithm for CI and orbital optimization schemes based on nonorthogonal orbitals. In contrast to previous nonorthogonal CI approaches (cf. Ref. 52), this newly proposed algorithm51 requires only the calculation of one- and two-particle reduced density matrices. In this paper, we derive the working equations to calculate one- and two-particle matrix elements between MPS wave functions that may be originally expressed in different, mutually nonorthogonal molecular orbital bases. Following the work of Malmqvist49 , the central element of our algorithm is the transformation of the bra and ket MPS wave functions to a biorthonormal basis representation. It is important to stress that the latter transformation is not needed if the MPS wave functions that are considered for state interaction share a common MO basis (cf. Refs. 22 and 23). We emphasize that our approach is applicable to the general case with MPS wave functions built from mutually nonorthogonal molecular orbital bases. It therefore provides the desired flexibility to find the best individual molecular orbital basis to represent wave functions of different spin and/or spatial symmetry. After solving a generalized eigenvalue equation of the form Hc = ESc ,

(1)

with the Hamiltonian matrix H expressed in the basis of the DMRG-SCF MPS wave functions and the wave function overlap matrix S we obtain a set of fully orthogonal and noninteracting states as linear combinations of the DMRG-SCF MPS wave functions with the expansion coefficients given by c. Since Malmqvist’s approach49 assumes either a full CI expansion or, in general terms, a ACS Paragon 4 Plus Environment

Page 5 of 37

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

Journal of Chemical Theory and Computation

wave function expansion that is closed under de-excitation53–55 we probe the closedness of our MPS wave function transformation by systematically increasing its numerical accuracy for a given active orbital space. In Section II we briefly discuss the theoretical framework for a second quantization formalism based on nonorthogonal orbitals. In Section III we introduce an algorithm to calculate expectation values for nonorthogonal wave functions in an MPS and matrix-product operator (MPO) representation of the wave function and operators56,57 , respectively, based on a nonunitary orbital transformation and demonstrate in Section IV how the latter can be exploited for a nonorthogonal MPS state-interaction (MPS-SI) ansatz. Numerical examples for the calculation of g-factors for f 1 - and f 2 -type actinide complexes are presented in Section V.

II.

SECOND QUANTIZATION FOR NONORTHOGONAL ORBITALS In this section, we briefly summarize the second quantization formalism for nonorthogonal

orbitals50,51,58,59 . In what follows, all quantities expressed in an orthonormal orbital basis are denoted by a tilde, those in the biorthonormal basis will be labeled by a bar whereas quantities with no extra label refer to a nonorthogonal orbital basis. Assuming an arbitrary set of N (linearly independent) MOs ϕ = {ϕp } with the general overlap matrix S = {Spq } Spq = hϕp | ϕq i =

Z

dr ϕ∗p (r)ϕq (r) ,

(2)

we can define a new set of orbitals ϕ ˜ = {ϕ˜p } that form an orthonormal basis by applying a L¨owdin symmetric orthogonalization60 ϕ ˜ = ϕ S−1/2 .

(3)

For an element of the overlap matrix in the new basis then holds S˜pq = S−1/2† SS−1/2



pq

= δpq ,

(4)

and consequently the new orbital set {ϕ} ˜ is orthonormal. The corresponding set of creation (annihilation) operators {˜ a†pσ } ({˜ apτ }) for spin orbitals, ϕ˜p (r)σ(ms ), in the new basis is ACS Paragon 5 Plus Environment

Journal of Chemical Theory and Computation

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

Page 6 of 37

related to the set of creation (annihilation) operators {a†pσ } ({apτ }) for spin orbitals of the original (nonorthogonal) basis59 , a ˜†pσ =

X

−1/2 a†qσ Sqp ,

(5)

−1/2 aqσ Spq .

(6)

q

a ˜pσ =

X q

Note that the creation {˜ a†pσ } and annihilation {˜ apσ } operators satisfy the well-known anticommutation rules. Expressing the creation (annhihilation) operators of the nonorthogonal basis in terms of the operators defined in the orthonormal basis yields59 a†pσ =

X

1/2 a ˜†qσ Sqp ,

(7)

1/2 a ˜qσ Spq .

(8)

q

apσ =

X q

Inserting Eq. (7) into the definition of the anticommutator, it is easy to verify that two creation operators in the nonorthogonal basis anticommute as it is the case in the orthonormal basis, {a†pσ , a†qτ } =

X

1/2 † 1/2 {˜ a†rσ Srp ,a ˜sτ Ssq } =

X rs στ

rs στ

1/2 1/2 Srp Ssq {˜ a† , a ˜† } = 0 , | rσ{z sτ }

(9)

=0

where we exploited in the second step the anticommutation of two creation operators for orthonormal spin orbitals. The same result holds for the anticommutation relation of two annihilation operators which is shown by Hermitian conjugation of Eq. (9). The anticommutator between a creation and annihilation operator reads {a†pσ , aqτ } =

X rs στ

1/2 1/2 {˜ a†rσ Srp ,a ˜sτ Sqs } =

X rs στ

1/2 1/2 Srp Sqs {˜ a† , a ˜† } = | rσ{z sτ } =δrs δστ

X

1/2 1/2 Sqr Srp δστ = Sqp δστ .

r στ

(10)

Hence, whereas the anticommutator for a general pair a ˜†pσ , a ˜†qτ of creation- and annihilation operators reduces in the orthonormal orbital basis to {˜ a†pσ , a ˜qτ } = δpq δστ ,

(11)

it depends on the (in general non-vanishing) overlap matrix element Sqp in the nonorthogonal case. ACS Paragon 6 Plus Environment

Page 7 of 37

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

Journal of Chemical Theory and Computation

As a consequence, the action of an annihilator aqσ on a given occupation number vector (ONV) |ki |ki =

L Y

a†pσ

p=1

with kpσ leads to a sum of ONVs59 aqσ |ki =

kpσ

|vaci = |k1 k2 . . . kL i ,

  1 if ϕ (r)σ(m ) occupied p s =  0 if ϕp (r)σ(ms ) unoccupied L X

(−1)

"

p−1 P j=1

kjσ

#

Sqp kpσ |k1 k2 . . . 0p . . . kL i ,

(12)

(13)

(14)

p

rather than to a single ONV as it would be the case for an orthonormal orbital basis. For the sake of completeness, we provide the proof for Eq. (14) in the Appendix which exploits the anticommutator relation in Eq. (10). In the last decades, considerable efforts were made to make nonorthogonal approaches such as valence-bond theory52,61–65 and nonorthogonal CI51,66,67 efficient. The latter suffer from an increasing computational complexity compared to a standard orthonormal formalism68–76 because of the nonorthogonal molecular orbital (MO) basis. In the framework of second quantization, the additional complexity can be attributed to the non-vanishing anticommutator in Eq. (10). For example, evaluating efficiently a matrix element of the type hΨ| O |Φi in a nonorthogonal approach, where O could be an arbitrary one-electron (two-electron) operator and |Ψi and |Φi are many-particle wave functions optimized for two different sets of MOs, respectively, requires a generalization68,69,72 of the Slater-Condon rules77,78 . In order to arrive at a formalism for evaluating matrix elements that closely resembles an orthonormal approach, the anticommutator in Eq. (11) must vanish. To this end, it is useful to define a new orbital basis ϕ ¯ = {ϕ¯p } through the transformation50,58,79 , ϕ ¯ = ϕ S−1 ,

(15)

where {ϕ¯p } is referred to as the dual of the nonorthogonal basis {ϕp }58,79 . Moreover, {ϕ¯p } and {ϕp } are said to form a biorthonormal system50 since we have hϕ¯p | ϕq i = δpq , ACS Paragon 7 Plus Environment

(16)

Journal of Chemical Theory and Computation

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

Page 8 of 37

and further, if both bases span the same space, they are jointly called a biorthonormal orbital basis. The definition of the biorthonormal creation operators follows from Eq. (15) as a ¯†pσ =

X

−1 a†rσ Srp .

(17)

r

To illustrate the biorthonormality, we consider the anticommutator for a biorthonormal creation operator with an annihilation operator in the original basis, {¯ a†pσ , aqτ } =

X

−1 {a†rσ Srp , aqτ } =

r

X r

−1 Srp {a†rσ , aqτ } = δpq δστ , | {z }

(18)

=Sqr δστ

which yields the standard anticommutation relation (cf. Eq. (11)).

III.

MATRIX PRODUCT STATES AND MATRIX PRODUCT

OPERATORS A.

Concepts We briefly introduce the concepts of expressing a quantum state as an MPS and a (Her-

mitian) operator as an MPO. Our notation follows the presentation of Ref. 56. Consider an arbitrary state |Ψi in a Hilbert space spanned by L spatial orbitals which we express as a linear superposition of ONVs |ki with the CI coefficients ck1 ...kL as expansion coefficients |Ψi =

X k

ck |ki =

X

ck1 ...kL |k1 . . . kL i ,

(19)

k1 ,...,kL

where each local space is of dimension four corresponding to the basis states kl = |↑↓i , |↑i , |↓i , |0i of a spatial orbital. In an MPS representation of |Ψi, we encode the CI coefficients ck1 ...kL l as a product of ml−1 × ml -dimensional matrices M kl = {Makl−1 al }

|Ψi =

X

X

k1 ,...,kL a1 ,...,aL−1

k1 L M1a Mak12a2 · · · MakL−1 1 |k1 . . . kL i = 1

X

M k1 M k2 · · · M kL |ki ,

(20)

k

where the last equality follows from collapsing the summation over the al indices (sometimes referred to as virtual indices or bonds) as matrix-matrix multiplications. Since the final contraction of the matrices M kl must yield the scalar coefficient ck1 ...kL , the first and the last matrices are in practice 1 × m1 -dimensional row and mL−1 × 1-dimenisional column vectors, ACS Paragon 8 Plus Environment

Page 9 of 37

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

Journal of Chemical Theory and Computation

respectively. Allowing for the introduction of some maximum dimension m for the matrices M kl , with m commonly referred to as number of renormalized block states 31 , is the central idea that facilitates a reduction of the exponentially scaling full CI ansatz in Eq. (19) to a polynomial-scaling MPS wave function ansatz. For further details on the actual variational search algorithm for ground- and excited states in an MPS framework, we refer the reader to the review by Schollw¨ock33 and, for its formulation in a quantum chemical program package, for instance to our recent works56,57 . c in MPO form We may express an operator W c= W

=

X

X

X

k k′

k k′

′ b ,...,b k1 ,...,kL k1′ ,...,kL 1 L−1

X



k k′

L L ′ ′ W1b11 1 Wb12b22 · · · WbL−1 1 |k1 . . . kL i hk1 . . . kL |





W k1 k1 W k2 k2 · · · W kL kL |ki hk′ | ≡

X

wkk′ |ki hk′ | ,

(21)

kk′

kk′

with the incoming and outgoing physical states kl and kl′ and the virtual indices bl−1 and bl . In analogy to Eq. (20), we may recognize the summation over pairwise matching indices bl as matrix-matrix multiplications which leads to the second line on the right-hand side of Eq. (21). For practical applications, we rearrange the summations in Eq. (21) and contract first over the local site indices kl kl′ cb b = W l−1 l

which then yields c= W

X

X

k k′

l l ′ Wbl−1 bl |kl i hkl | ,

(22)

1 W1b · · · Wbll−1 bl · · · WbLL−1 1 . 1

(23)

kl kl′

b1 ,...,bL−1

The local, operator-valued matrix representation introduced in Eq. (22) is a central element for an efficient MPO-based implementation of the quantum-chemical DMRG approach56,57 that offers the same polynomial scaling as a “traditional” (non-MPO) DMRG implementation.

B.

MPO expectation values with nonorthogonal orbitals The calculation of overlap matrix elements and expectation values for N -electron oper-

ators in an MPO framework based on an orthonormal basis common for the bra and ket states was outlined, for instance, in Refs. 33 and 56, respectively. Here, we illustrate the ACS Paragon 9 Plus Environment

Journal of Chemical Theory and Computation

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

Page 10 of 37

complexity that results when the bra and ket states are expressed in two different MO bases which are in general not orthonormal. Following the notation by Malmqvist and Roos48 , orbital basis superscripts X and Y Y denote the original (orthonormal) MO bases {ϕX p } and {ϕp }, respectively, whereas the B superscripts A and B refer to different MO bases {ϕA p } and {ϕp }, respectively, which have

been manipulated in some way and are therefore in general not orthonormal. We assume that the MPS wave functions are optimized with the same active orbital spaces and a common atomic orbital basis set. The latter restrictions can, however, be lifted50,80 .

1.

General considerations Let |Ψi and |Φi denote two MPS wave functions based on the definition in Eq. (20), |Ψi =

X k

|Φi =

A

X

X

a1 ,...,aL−1

X

A A A kL kA , M1a11 Mak12a2 · · · MaL−1 1 k kB

kB

kB

N1a1′ Na′2a′ · · · Na′L 1

1 2

kB a′1 ,...,a′L−1

L−1 1

(24)

B k ,

(25)

B which are constructed from MO sets {ϕA p } and {ϕp }, respectively.

We start by considering the overlap hΦ| Ψ i which can be written as33,56,59   A X X  kB † A kL k1 k1B † k2B † k2A L M1a1 Ma1 a2 · · · MaL−1 1 kB kA i N1a′ · · · Na′ a′ Na′ 1 hΦ| Ψ i = kA kB a1′ ,...,aL−1 a1 ,...,a′L−1

=

X

X

kA kB a1′ ,...,aL−1 a1 ,...,a′L−1



kB † N1aL′ L−1

1

2 1

L−1

kB † kB † · · · Na′2a′ Na′11 1 2 1

×δNkB NkA det Sk

B A

k

 A  A kL k1 k2A M1a1 Ma1 a2 · · · MaL−1 1

(26)

,

with NkB and NkA being the number of electrons comprised in the ONVs and det Sk determinant of the overlap matrix Sk

B A

k

B A

k

the

of all overlap integrals between occupied orbitals of

the bra and ket ONVs. Compared to the standard expression of the inner product of ONVs B 59 that are built from a common orthonormal basis {ϕ}, ˜ implying for example {ϕA p } ≡ {ϕp },



k



B

A

k i=

L Y

δkpB kpA ,

(27)

p=1

the inner product in Eq. (26) leads to a rather complex expression with a large number of non-zero terms. The reason for the latter is that in the general case of nonorthogonal ACS Paragon10 Plus Environment

Page 11 of 37

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

Journal of Chemical Theory and Computation

orbitals the action of an annihilator on a given ONV, as illustrated in Eq. (14), creates a sum of ONVs rather than a single ONV.

ˆ |Ψi of an N -electron operator O ˆ given Turning next to transition matrix elements hΦ| O in MPO form (cf. Eq (21)), the expression in the case of nonorthogonal orbitals reads

ˆ |Ψi = hΦ| O

X

X

kA kB a1′ ,...,aL−1 a1 ,...,a′L−1





kB † N1aL′ L−1

X

× kB  ×



X

kB † kB † · · · Na′2a′ Na′11 1 2 1

X

kB′ kB′′ kA′ b1 ...bL−1

A kA M1a11 Mak12a2

A kL · · · MaL−1 1





 

kB′ kA′ kB′ kA′ kB′ kA′ O 1 1 O 2 2 · · · O L L kB′ kB′′ kA′ kA′  1b1

b1 b2

 k A .

bL−1 1

(28)



Note that, similarly to Eq. (26), the inner product kB′′ kA′ of the ONVs will not sim-

ply reduce to a product of delta functions since the incoming and outgoing basis states of the operator are nonorthogonal. Moreover, the fermionic anticommutation which is encoded explicitly in the matrix representation of the creation and annihilation operators (see, in particular, Section IV in Ref. 56 for a detailed discussion) by using their WignerJordan-transformed form81 will introduce additional overlap terms as becomes clear from

the anticommutation rules in Eqs. (9) and (10).

If, however, we assume that |Ψi and |Φi are expressed in a common, orthonormal basis ACS Paragon11 Plus Environment

Journal of Chemical Theory and Computation

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

Page 12 of 37

where Eq. (27) holds, Eq. (28) will reduce to the well-known expression33,56 ,  X X  kB † kB † kB † ˆ |Ψi = N1aL′ · · · Na′2a′ Na′11 hΦ| O kA kB a1′ ,...,aL−1 a1 ,...,a′L−1



X

× kB 

X

1

2 1

L−1

X

kB′ kB′′ kA′ b1 ...bL−1



  B′ B′′ A′ A′  B′ k A′ kL kB′ kA′ kB′ kA′ L k k k O1b11 1 Ob12b2 2 · · · ObL−1 k  1 | {z } =δ

B′′ A′ =I

k k  A A k k1A L k A × M1a1 Mak12a2 · · · MaL−1 1    B′ A′ B′ A′ X X X X  kB † B′ k A′ kL k k k k kB † kB † L O1b11 1 Ob12b2 2 · · · ObL−1 = N1aL′ · · · Na′2a′ Na′11 1



kA′ kA kB′ kB b1 ...bL−1 a1′ ,...,aL−1 a1 ,...,a′L−1

2 1

L−1

1

  A A A kL k B B′ ih kA′ | kA i × M1a11 Mak12a2 · · · MaL−1 1 h k |k

=

X

B kA kL L aL−1 a′L−1 bL−1



=δkB kB′

=δkA′ kA

   X X B kA kB † kL kA kB † kB kA kB † kB kA L  Na′11 O1b11 1 M1a11  N1aL′ ObL−1 Na′2a′ Ob12b22  · · · 1 1 2 1 L−1  B A B A k2 k 2 a1 a′1 b1



k 1 k1

 A  kLA × Mak12a2 · · ·  MaL−1 1 . 

(29)

Note that the last equality follows from the original expression by regrouping the summations to minimize the operational cost33,56 .

2.

Transformation of orbitals and MPS wave functions To bring the expressions for the overlap and transition matrix elements, Eqs. (26) and

(28), into a form that is comparable to the case for a common orthonormal basis requires the fulfillment of the biorthonormality condition58,79

B ϕ¯A p ϕq i = δpq .

(30)

The price to pay is, however, that a transformation to a biorthonormal MO basis entails an additional transformation step of the wave function expansion parameters, in the present case the MPS tensors. This is most easily seen by inspection of Eq. (17) which defines a creation operator in the biorthonormal basis as a suitable linear combination of creation operators in the original basis. Because of the latter basis change of the creation (annihilation) operators ACS Paragon12 Plus Environment

Page 13 of 37

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

Journal of Chemical Theory and Computation

the wave function expansion parameters need to be transformed as well. Moreover, as already pointed out by Malmqvist and Roos48 , the objective is to find efficient wave function Y A B transformations from the original basis {ϕX p } ({ϕp }) to the biorthogonal basis {ϕp } ({ϕp })

while explicitly avoiding any additional wave function optimization. Such an algorithm was proposed by Malmqvist for CI-type wave functions49 . Following up on Malmqvist’s original idea, we outline in the following how such a transformation can be achieved for MPS-type wave function expansions by a sequence of single-orbital replacements. Similarly to CI-type expansions, each step in the sequence is approximately twice as expensive to compute as the matrix-vector product y = Tb x .

for a one-electron-operator Tb . a.

(31)

Orbital transformation Considering the biorthonormality condition (cf. Eq. (30)),

Malmqvist49 and Olsen et al.50 showed that an LU-partitioning of the inverse of the orbital −1 defined according to Eq. (2) with the bra (ket) orbitals in {ϕX overlap matrix SXY p}

({ϕY p }) yields

SXY

−1

= CYB CXA

†

,

(32)

where CXA and CYB are the desired orbital transformation matrices. The LU-partitioning −1 has an arbitrary degree of freedom and the of the inverse of the overlap matrix SXY only requirement is that Eq. (32) is fullfilled in each step of a row-wise LU-partitioning −1 of SXY 49 . Additionally, no permutations of rows and columns are allowed to avoid a

recoupling of orbitals. We strictly employ the algorithm described in detail in Section V of Ref. 49 which proceeds in each step of row-wise LU-partitioning such that the final CXA and CYB matrices are numerically best suited for a subsequent LU-partitioning (vide infra). The non-unitary orbital transformation and subsequent counterrotation of the MPS wave function presented in this work will not change the physical content of our wave function. In addition, having qualitatively similar MOs in the active space but appearing in different order is of no concern for our approach. The only requirement is that the size of the active spaces, namely number of MOs and number of electrons, for the left- and righthand side MPS match. As shown on the left-hand side of Fig. 1, the MO overlap matrix between the optimized sets of active MOs for the triplet and singlet wave functions of the PuO2+ 2 molecule (see Section V B for more details) contains significant departures from a ACS Paragon13 Plus Environment

Journal of Chemical Theory and Computation

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

Page 14 of 37

unit matrix which, however, does not pose a particular problem for the subsequent orbital and MPS transformation. In practice, active orbital spaces describing the same molecule but in electronic states of different spin and/or spatial symmetry will be qualitatively the same if the active orbital space can be chosen to be sufficiently large. In this regard, having the possibility to treat large active spaces as it is the case for the DMRG-optimized MPS parametrization will be beneficial for the numerical stability of the transformation to a biorthonormal basis representation and of the subsequent state-interaction approach outlined in Section IV. Transforming the wave function to a new orbital basis can introduce changes in orbital entanglement measures. To this end, we analyzed the single-orbital entropy si (1)82 for the ith orbital and the mutual information Iij 83,84 for the i-th and j-th orbital in the active space. As an example, the right-hand side of Fig. 1 illustrates the mutual information difference matrix ∆(I˜ij − I¯ij ) obtained for the lowest-lying triplet state of PuO2+ 2 (see Section V B for computational details) in orthonormal and biorthonormal orbital bases, respectively. Hence, as expected, a non-unitary orbital transformation gives rise to changes in the mutual information. By contrast, the single-orbital entropy si (1), which can be employed as a descriptor to determine in an automated fashion the most suitable active orbital space26,28 , remains constant up to a deviation by 10−4 . The latter in turn indicates that the algorithm outlined in Ref. 28 would lead to the proposition of the same active orbital space for PuO2+ 2 irrespective of the molecular orbital basis. Hence, the proposed non-unitary orbital transformation is not likely to lead to an artificial creation of entanglement that would hamper the fidelity of our MPS transformation. Having found suitable transformation matrices CXA and CYB according to Eq. (32) and the algorithm outlined in Ref. 49, allows us to express the biorthonormal bases {ϕA p } and X Y 49,50 {ϕB p } in terms of the original orbital bases {ϕp } and {ϕp }, respectively,

ϕA = ϕX CXA ,

(33)

ϕB = ϕY CYB .

(34)

Expressing the orbital transformation in Eq. (33) as a sequence of single orbital transforma49,50 tions the substitution of orbital ϕX 1 reads A X X ϕX 1 = ϕ1 t11 + ϕ2 t21 + ϕ3 t31 + . . . ,

ACS Paragon14 Plus Environment

(35)

Page 15 of 37

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 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 16 of 37

The same considerations hold for CYB and its factorization into upper and lower triangular matrices. b. MPS wave function transformation With the elements of the matrix t at hand (vide supra), the wave function expansion parameters can be transformed by a sequence of singleorbital transformations49,50 . A detailed account of this transformation approach can be found for CI-type wave functions in Refs. 49 and 50. Following closely their ansatz, we outline below the essential steps required to transform an MPS-type wave function representation A from an expansion in the original basis {ϕX p } to an expansion in a biorthonormal basis {ϕp }.

As discussed in Subsection III B, the latter will allow us to calculate N -particle transition density matrices in a framework of standard second-quantization algebra. In our pilot implementation presented in this work, we chose for simplicity to carry out the biorthonormal transformation for a non-spin adapted MPS with a given initial spin S and spin projection MS . This MPS is obtained after a transformation of the MPS from a spinadapted (SU(2)) to a non-spin adapted representation which generates an MPS for each of the possible MS values ranging from –S, −S +1, . . . , S −1, S. All transition densities are then evaluated for an MPS of a given MS value. Exploiting the Wigner-Eckart theorem5,85 allows us to generate a so-called Wigner-Eckart reduced transition density matrix (see Section IV for more details) from which we proceed to calculate matrix elements of different operators D E such as the spin-orbit mean-field operator. By calculating the expectation value Sˆ2 of the squared spin operator Sˆ2 for a given transformed MPS we confirmed that the transformation into the biorthonormal basis preserved the SU(2) symmetry of the initial spin-adapted MPS with a given spin S up to a numerical accuracy of 10−4 a.u. in all cases presented in this work. We start by transforming the wave function |Ψi |Ψi =

X kX

X X X M k1 M k 2 · · · M kL k X ,

(41)

where we have introduced above an additional superscript X to emphasize that the MPS tensors of |Ψi refer to the original orbital basis {ϕX }. With the parameters tL and tU calculated from an LU-factorization of CXA according to Eq. (40), the MPS wave function transformation proceeds differently for inactive and active orbitals. No special action is required for the secondary (virtual) orbital space. ACS Paragon16 Plus Environment

Page 17 of 37

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

Journal of Chemical Theory and Computation

Transformation with respect to the inactive orbital space As shown in Ref. 49, this transformation step reduces for CAS-type wave functions to a simple scaling of |Ψi by a factor α, |Ψi ≡ α · |Ψi =

X



α· M

kX

k1X

M

k2X

···M

X kL

 k X

(42)

with α=

nI Y

(tii )2 ,

(43)

i=1

where the inactive orbital space comprises nI orbitals.

Transformation with respect to the active orbital space Assuming an active orbital space comprising L active orbitals, we set the orbital counter to j = 1 and proceed as follows: X

1. Scale the j-th MPS tensor M kj of |Ψi consisting of a set of matrices (one for each basis state occupation kjX ) with respect to the occupation number of the j-th orbital,

X

M kj

and X

 X kj,|↑↓i 2  t · M  jj   X   t · M kj,|↑i jj ≡ X   tjj · M kj,|↓i    X  M kj,|0i X

M ki ≡ M k i

for kjX = |↑↓i for kjX = |↑i

,

(44)

∀ i 6= j ∧ i = 1, . . . , L .

(45)

for kjX = |↓i for kjX = |0i

c for the one-electron operator Tb , 2. Construct a new MPO W Tb =

L X tmj † a ajσ , tjj mσ m6=j σ

(46)

where the coefficients tmj are given by the matrix elements of tL and tU , respectively. Note that σ denotes here the eigenvalue of the spin eigenfunction. 3. Carry out a sequence of MPO–MPS operations33 , starting with the exact product ACS Paragon17 Plus Environment

Journal of Chemical Theory and Computation

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

Page 18 of 37

(analogous to Section 5.1, pages 140–141 of Ref. 33) (1) c |Ψi Ψ =W   X′  X  X X′ X X′ X X′ X′ X′ = W k1 k1 W k 2 k 2 · · · W k L kL M k 1 M k2 · · · M kL k X kX kX′

=

X X

kX k1X′

W1b11

kX kX′ ab

=

X X

kX kX′ kX′ W1b11 1 M1a11

kX kX′ ab

=

XX kX

=

X k

X

ab

kX kX′

kX kX′

L L Wb12b2 2 · · · WbL−1 1

   X′ X′ X′ kL k k X M1a11 Mak12a2 · · · MaL−1 1

   X X′   X X′ X′ kL k L kL k2 k 2 k2X′ Wb1 b2 Ma1 a2 · · · WbL−1 1 MaL−1 1 kX

X X kL k1X k2X k N(11)(b N · · · N (b a )(11) a ) (b a )(b a ) 1 1 1 1 2 2 L−1 L−1

X X X N k1 N k2 · · · N k L k X ,

(47)

where the compact notation

kX

N(bii−1 ai−1 )(bi ai ) =

X

kX kX′

kX′

i i i Wbi−1 bi Mai−1 ai ,

(48)

kiX′

has been introduced in the second to last step on the right-hand side of Eq. (47). 4. Compress Ψ(1) by means of a singular-value decomposition (SVD) with a truncation threshold ǫ of (at most) ǫ = 10−8 . This is followed by

(2) c Ψ(1) Ψ =W  X X  X  X X′ X X′ X X′ X = W k1 k 1 W k 2 k2 · · · W kL kL N k1 N k2 · · · N k L k X kX kX′

=

XX kX

=

X k

X

ab

X X kL k2X k1X · · · O O O(11)(b (bL−1 aL−1 )(11) k (b1 a1 )(b2 a2 ) 1 a1 )

X X X O k1 O k2 · · · O k L k X ,

(49)

in complete analogy to Eq. (47).

5. Set |Ψi ≡ |Ψi + Ψ(1) + 12 Ψ(2) and perform an SVD compression of |Ψi.

6. Increase the orbital counter j by one. If j ≤ L repeat steps (1)-(5), otherwise exit the transformation algorithm. Note that the final MPS tensors of |Ψi now refer to an expansion in the biorthonormal basis {ϕ¯A }, for example, A

X

M k i ≡ M ki

∀ i = 1, . . . , L ,

ACS Paragon18 Plus Environment

(50)

Page 19 of 37

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

Journal of Chemical Theory and Computation

such that we can express |Ψi as |Ψi =

X k

A

A A A M k 1 M k 2 · · · M kL k A

(51)

Trial studies revealed that a rather conservative cutoff value of ǫ = 10−8 for the SVD compression required in steps 4 and 5 is sufficient to provide numerical stability in the transformation step to biorthonormal basis. Moreover, the conservative compression threshold ensures that the numerical accuracy of the MPS wave function is kept after the transformation to the biorthonormal basis. As a numerical test, we calculated for the lowest triplet state of PuO2+ discussed in Section V B the one-electron part of the Hamiltonian in the 2 original and biorthonormal bases from the MPS and corresponding CI wave functions. The energies obtained from the MPS and CI wave functions were identical up to 10−8 Hartree before and after the wave function transformation. Having found a representation of |Ψi in the biorthonormal basis {ϕ¯A p }, we repeat the transformation steps above for the inactive and active orbital spaces of the MPS wave function |Φi with tL and tU calculated from an LU-factorization of CYB .

IV.

THE STATE-INTERACTION APPROACH FOR MPS WAVE

FUNCTIONS With the overlap and transition matrix elements at hand, calculated from spin-free MPS wave functions (i.e., SU(2) invariant MPS wave functions |Ψ(S)i with a well-defined total spin quantum number S as, for instance implemented in the QCMaquis program57 ) in a biorthonormal basis, we are now able not only to calculate Hamiltonian but also matrix elements for arbitrary one- and two-particle property operators such as transition dipole moments (from which we can calculate oscillator strengths), angular momentum eigenvalues and magnetic transition dipoles, or electric field gradients. Based on the discussion in the previous section, Fig. 2 illustrates the typical workflow of our MPS-SI approach where the MPSs are obtained from different DMRG-SCF optimizations. We implemented this scheme in a development version of the Molcas 8.086 software which, to a large extent, exploits the existing framework of the CASSI/RASSI implementation for CI-type wave functions by Malmqvist and co-workers48,80 . Further details on the SI approach common to both CI-type ACS Paragon19 Plus Environment

Journal of Chemical Theory and Computation

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

Page 20 of 37

and MPS wave functions can be found in Refs. 48 and 80 and we here provide only a brief summary of the most important steps as outlined in Fig. 2. Considering only the lower triangular matrix of dimension N (N + 1) originating from N converged spin-free (non- or scalar-relativistic) MPS wave functions of arbitrary spatial symmetry and total spin S, we first test whether the i-th MPS shares the same MO basis as the j-th state. If the outcome of the test is negative, both MPSs are transformed to the biorthonormal basis according to Section III B 2 b, before proceeding with the calculation of the wave function overlap and transition density matrix elements. The latter are then combined with the corresponding AO-integrals to evaluate property and Hamiltonian matrix elements. In the last step we diagonalize the Hamiltonian matrix and calculate property matrix elements over the resulting N spin-free eigenstates. Moreover, the calculation of magnetic properties such as g-tensors and hyperfine coupling tensors as well as the prediction of intersystem crossing rates requires access to SO coupled wave functions. To this end, we diagonalize the full Hamiltonian ˆ=H ˆ el,sf + H ˆ SO , H

(52)

represented in the basis of the multiplet of states |Ψ(S, M )i that can be obtained from a spin-free calculation of |Ψ(S)i. Here, the state label (S, M ) indicates the pair of total spin S and spin magnetic quantum number M the latter which can take values in the range from −S ˆ el,sf is typically a scalar-relativistic, to S with unit increments. In our current approach, H spin-free electronic Hamiltonian, for example the Douglas-Kroll Hess Hamiltonian87–89 or the exact two-component Hamiltonian (X2C)90–93 with the energy expectation value EA ˆ SO we consider an effective one-electron for the eigenfunction ΨA (S, M ) of state A. For H mean-field SO Hamiltonian94 in which the two-electron terms essentially serve as screening

corrections of the (dominating) one-electron terms. Without significant loss of accuracy18,95

we apply an atomic-mean field approximation of the SO integrals as implemented in the AMFI program of B. Schimmelpfennig96 where all multicenter integrals are neglected. In the AMFI representation the effective one-electron mean-field SO Hamiltonian reads80 ˆ SO = H

X

z ˆz ˆx + V y T ˆy Vpqx T pq pq pq + Vpq Tpq

pq



,

(53)

where Vpq are the AO spin-orbit interaction integrals and Tpq the triplet operators in Cartesian representation. The latter are related to the more common triplet operators in spinACS Paragon20 Plus Environment

Page 21 of 37

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 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

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

Page 22 of 37

ˆ S,M is a spin-tensor operator whose eigenfunctions are a tensor state with spin Here, T eigenvalues S and M . It is one of the 2S + 1 spin-tensor operators comprising the spinˆ S of half-integer or integer rank S. To calculate the matrix representation tensor operator T of the SO operator given in Eq. (53) we need to evaluate matrix elements between state D E A ′ ˆ S,M ΨB (S ′′ , M ′′ ) . These Ψ (S , M ′ ) and state ΨB (S ′′ , M ′′ ) of the form ΨA (S ′ , M ′ ) T can be calculated very efficiently by employing the Wigner-Eckart (WE) theorem5,85 which

states that the (2S ′ + 1) × (2S + 1) × (2S ′′ + 1) matrix elements above can be obtained as a S E D ˆ B ′′ product from a single (WE-reduced) matrix element ΨA (S ′ ) T Ψ (S ) and a Wigner 3j-symbol (calculated from Clebsch-Gordon coefficients), D E ′′ ′ ˆ S,M B ′′ ′′ ΨA (S ′ , M ′ ) T Ψ (S , M ) = (−1)S +M −S   S ′′ ′ D E S S S ˆ ΨB (S ′′ ) . (58)  ΨA (S ′ ) T × M ′′ M M ′

Taking advantage of the WE theorem, we calculate WE-reduced one-particle (spin) transition density matrices80

D E A ′ ˆ S B ′′ ΓWE = Ψ (S ) T Ψ (S ) , pq pq

which can now be employed to calculate a WE-reduced matrix element80 X VAB = ΓWE pq Vpq .

(59)

(60)

pq

Since the matrix elements over spin states will only be non-zero if |S ′′ − S ′ | is an integer ≤ 1, while also |M ′′ − M ′ | is an integer ≤ 1, we are left with a total of nine non-zero SO Hamiltonian matrix elements, for example, p D E  ˆ SO ΨB (S + 1, M ± 1) = − (S±M + 1) (S±M + 2) ±V AB,x + ıV AB,y ΨA (S , M ) H 2 (61)

and

D

E ˆ SO B Ψ (S , M ) H Ψ (S , M ) = M V AB,z , A

(62)

where the remaining six cases can be found in Ref. 80.

With the matrix elements of the SO operator at hand, we next diagonalize the total Hamiltonian matrix H of the Hamiltonian H given in Eq. (52) with the matrix elements given by D E D E ˆ B ′′ ˆ SO B ′′ ΨA (S ′ , M ′ ) H Ψ (S , M ′′ ) = SAB δS ′ S ′′ δM ′ M ′′ EA + ΨA (S ′ , M ′ ) H Ψ (S , M ′′ ) ,

(63)

ACS Paragon22 Plus Environment

Page 23 of 37

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

Journal of Chemical Theory and Computation

where SAB is the overlap between states A and B. Note that it is also possible, prior to the diagonalziation of H, to correct the state energies for dynamic electron correlation contributions from, e.g., a DMRG-CASPT297–99 or DMRG-NEVPT2100–104 calculation by shifting the elements of the total Hamiltonian by80 ∆HAB = 0.5 (∆EA + ∆EB ) SAB ,

(64)

where ∆EA is the dynamical electron correlation contribution shift for state A and similarly for state B. Diagonalization of H yields a set of SO coupled eigenstates which are linear combinations of all MPS wave functions with different spin and magnetic quantum numbers and the corresponding eigenvalues. In the final step (lower box in Fig. 2), we can now calculate the desired property matrix elements such as g-tensor matrix elements (see for example Ref. 9 for further details) in the basis of the SO eigenstates.

V.

NUMERICAL EXAMPLES

A.

Closedness of MPS wave function expansion The transformation algorithm outlined in Section III B 2 b is based on Malmqvist’s

approach49 which can only be applied to wave function expansions that are closed under de-excitation53–55 . The latter implies that in each sequence of the wave function transformation no configurations outside of the original wave function expansion are generated. This prerequisite is fulfilled by full CI or CAS-CI wave functions and was also shown to be the case for restricted active space wave functions53 . Here, we probe the closedness and (approximate) orbital rotation invariance of our MPS wave function transformation by systematically increasing the numerical accuracy of the MPS with a varying number of renormalized block states m for a given active orbital space since for finite m, the MPS ansatz will only approximately represent the corresponding CAS-CI wave function. As an example, we choose to study the zero-field splitting (ZFS) of the spin-free ground 105 state 3 Σ− g of the tellurium dimer Te2 for which CASSCF-SO reference data are available 1 + for both the ZFS and the SO matrix element V between the spin-free 3 Σ− g and Σg states.

In accordance with the reference work of Rota et al.105 , we employ an all-electron ANORCC basis set of TZP quality106,107 for Te along with the 2nd-order Douglas-Kroll-Hess ACS Paragon23 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

(DKH2) scalar-relativistic Hamiltonian108 . The internuclear Te-Te distance of 2.557 ˚ A is taken from experiment109 . In the chalcogenide dimers, the low-lying electronic spectrum originates from the valence (π ∗ )2 configuration which gives rise to the X 0+ g electronic ground and X2 1g , a 2g and b 0g + excited states. We consider a CAS(8,6) with 8 electrons allowed to freely occupy the six spatial orbitals (σ, σ ∗ , π and π ∗ ) resulting from combinations of the atomic Te 5p valence orbitals. For this active orbital space, we carried out individual state-averaged spin-free DMRG-SCF calculations for five triplet and six singlet states which are subsequently allowed to mix through SO interaction in our MPS-SI ansatz denoted as DMRG-SCF-SO. The latter requires a transformation to a biorthonormal MO basis and a corresponding rotation of the MPS wave functions to this basis. We therefore expect that any significant violation of the closedness and/or orbital rotation invariance of the MPS ansatz for finite values of m would lead to considerable differences in the calculated spectrum and/or the SO matrix element V compared to the reference CASSCF data. All results are compiled in Table I. We first note that for the small CAS(8,6) space all our DMRG-SCF-SO data for the ZFS and the SO matrix element V , irrespective of the number of renormalized block states m ranging from m = 64 to m = 2048, is in excellent agreement with the reference CASSCF-SO data (obtained from spin-orbit state-interaction calculations for the spin-free CASSCF triplet and singlet states; third-to-last row in Table I) calculated in this work. This result may not be surprising at first glance because of the size of the considered CAS but it nevertheless indicates that the closedness assumption holds for full CI-type wave functions represented by MPS expansions with finite m values. Moreover, orbital rotation invariance which is in principle only fullfilled exactly for MPS wave function expansions with m → ∞ appear not to have an impact on our MPS transformation approach to the biorthonormal basis. Compared to the CASSCF-SO data by Rota et al., our DMRG-SCF-SO/CASSCF-SO data deviate by several hundred wave numbers for the excitation energies of the spin-free and SO coupled excited states as well as by 25 cm−1 for the SO matrix element V . This is partly due to the fact that we took into account only the lowest five (six) triplet (singlet) states in the wave function optimization. Considering all spin-free states up to 70000 cm−1 (second-to-last row in Table I) partly reduces the observed deviations, in particular for V . In summary, we emphasize that our CASSCF-SO and DMRG-SCF-SO data are consistent. ACS Paragon24 Plus Environment

Page 24 of 37

Page 25 of 37

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

Journal of Chemical Theory and Computation

TABLE I. Relative first excited state energies of Te2 obtained from spin-free (SF) DMRG-SCF and spin-orbit coupled (SOC) DMRG-SCF state interaction calculations denoted as DMRG-SCFSO. For comparison, corresponding spin-orbit coupled CASSCF state interaction calculations (CASSCF-SO) were carried out. The number of renormalized block states m in the DMRGSCF calculations was varied from m = 64 to m = 2048. V denotes the spin-orbit coupling matrix 1 + −1 element between the 3 Σ− g and Σg states. All data are given in cm .

SF

B.

method

CAS

m

1∆

DMRG-SCF-SO

(8,6)

64

3712 6588

1759 5471 10107

3832

DMRG-SCF-SO

(8,6)

128

3711 6587

1759 5471 10106

3832

DMRG-SCF-SO

(8,6)

256

3711 6587

1759 5470 10105

3832

DMRG-SCF-SO

(8,6)

512

3711 6586

1759 5470 10105

3832

DMRG-SCF-SO

(8,6)

1024

3711 6587

1759 5470 10105

3832

DMRG-SCF-SO

(8,6)

2048

3711 6586

1759 5470 10105

3832

CASSCF-SOa

(8,6)

-

3711 6586

1759 5470 10105

3832

CASSCF-SOb

(8,6)

-

3507 6298

1657 4917 9317

3808

CASSCF-SOc

(8,6)

-

2716 6009

1700 4180 9153

3807

g

1 Σ+ g

SOC X2 1g a2g b 0+ g

V

a

This work.

b

This work. All spin-free states up to 70,000 cm−1 were considered.

c

Ref. 105; CAS(8,6)SCF-SO/TZP data.

Magnetic resonance properties of sample f 1 and f 2 actinide complexes To illustrate the capabilities of our MPS-SI approach, we calculate g-factors of pro-

totypical f 1 - and f 2 -type actinide complexes, namely NpO2+ and PuO2+ 2 2 , and compare to CASSCF/CASPT2 data from Gendron et al.110 . To facilitate a comparison, we follow the same computational protocol as outlined in Ref. 110 which we briefly summarize here. We employed all-electron ANO-RCC basis sets of TZP quality106,107 along with the scalarrelativistic DKH2 Hamiltonian108 and a Cholesky-decomposition (CD) of the two-electron repulsion integrals as implemented in Molcas86 (keyword Cholesky in the integral module SEWARD). Dynamical electron correlation is included through CD-DMRG-NEVPT2 ACS Paragon25 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

(second order n-electron valence perturbation theory) calculations43,104 with a state-averaged DMRG-SCF reference wave function. Active orbital spaces for the latter comprise 7 (8) elec2+ trons in 10 orbitals for NpO2+ 2 (PuO2 ) denoted in the following as CAS(7,10) (CAS(8,10)). 2+ For NpO2+ 2 , we simultaneously optimized six doublet states while for PuO2 we carried out

individual state-averaged DMRG-SCF calculations for 20 triplet and 20 singlet states, respectively. In all DMRG-SCF calculations, the number renormalized block states m was set to 1024 which is sufficient to yield results of CASSCF quality for the active orbital spaces considered. In accordance with Gendron et al.110 , C1 point group symmetry was assumed in calculations. As outlined in Section IV, our MPS-SI approach is integrated in the CASSI/RASSI framework48,80 of Molcas which allowed us to calculate the SO operator matrix elements from the contraction of WE-reduced spin transition density matrices with atomic mean-field integrals94,96 while EPR g-factors were calculated based on the approach of Bolvin9 . In the case of Np(VI)O2+ we carried out additional four-component EPR calculations 2 using a multi-reference CI (MRCI) approach11 based on the Dirac-Coulomb Hamiltonian and all-electron uncontracted cc-pVTZ111 and dyall.v3z112 basis sets for oxygen and neptunium, respectively. MP2 natural spinors obtained from correlating the Np 5s5p5d6s6p and O 2s2p electrons while keeping the remaining core electrons of Np and O frozen constitute the orbital basis for the MRCI expansion. In the MP2 calculation, a virtual spinor threshold was set to 20 hartree such that the virtual correlation space comprised all recommended core-valence and valence-correlation functions. The MRCI active space includes all occupied spinors with occupation numbers less than 1.98 as well as all virtual spinors with occupation numbers up to about 0.001. The active space was further split in three subspaces with a maximum allowed excitation level of singles and singles-doubles from the first space containing all spinors with occupation larger than 1.80 into a model space comprising the four partially 110 occupied nonbonding Np 5f spinors of Np(VI)O2+ . From the combined two lower spaces 2

a maximum of singles-doubles-triples excitations are allowed into the third space comprising in total 90 spinors. Table II summarizes our property calculations for NpO2+ and PuO2+ 2 2 . Starting with NpO2+ 2 , we note that our DMRG/NEVPT2-SO results are in close agreement with the corresponding CASSCF/CASPT2-SO data from Gendron et al.110 both for the g-factors, gk and g⊥ , respectively, and the vertical excitation energy ∆E from the 2 Φ5/2u ground to the 2 ∆3/2u ACS Paragon26 Plus Environment

Page 26 of 37

Page 27 of 37

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

Journal of Chemical Theory and Computation

first excited state. The DMRG/NEVPT2-SO gk = 4.235 value for NpO2+ 2 compares further well with our reference four-component MRCI gk = 4.283 value which indicates that the considered number of spin-free states is sufficient to achieve convergence of the SO coupling contributions. Moreover, we find that the ground- and excited state compositions of the SO coupled 2 Φ5/2u and 2 ∆3/2u states in terms of their φ, δ and π character are well reproduced in our DMRG/NEVPT2-SO calculation in comparison to the CASSCF/CASPT2-SO reference data, in particular the sizable admixture of 12% δ character in the SO coupled 2 Φ5/2u ground state which originates from a strong coupling of the spin-free 2 Φu and 2 ∆u states110 . Turning next to PuO2+ 2 , we find an excellent agreement of our DMRG-SCF-SO results for the electronic Ω = 4g ground and the lowest-lying excited states with the corresponding CASSCF-SO reference data by Gendron et al.110 . We are not only able to correctly reproduce the ground-state g-factors of gk = 6.076 and g⊥ = 0.000 but also find matching vertical excitation energies ∆E with deviations of at most 5 cm−1 for the first three excited states. Similar to NpO2+ 2 , the lowest spin-free states of the plutonyl ion originate from occupations of the non-bonding δ and φ orbitals with different combinations (parallel or antiparallel) of spin (MS ) and angular momentum (ML ) projections. In accordance with Hund’s rules, the spin-free ground state is a 3 Hg (δ 1 φ1 occupation) triplet state with parallel angular momentum projections (ML = ±2 ± 3 = ±5) and antiparallel spin projection MS = ∓1. Further combinations of MS and ML within the δ 1 φ1 occupation manifold lead to 3 Σ− g and 3 Πg states which we find at 3875 cm−1 and 6563 cm−1 above the 3 Hg ground state, respectively. Moreover, according to our DMRG-SCF calculations the first spin-singlet state 1

−1 Σ+ above the electronic ground state. g is located at 10384 cm

At the SO level, the 3 Hg term splits into its three components with Ω = 4g , 5g , and 6g which are then allowed to mix with states of the same Ω. Although the Ω = 4g ground state remains predominantly of 3 Hg character (96%), it exhibits an admixture of the 1 Γg state (2%) which is energetically separated from the 3 Hg state by about 15300 cm−1 in the spinfree framework. In addition, as can be seen from Table II, the first two SO-coupled excited −1 states of Ω = 0+ and 6985 cm−1 , respectively, originate g and Ω = 1g at ∆E = 4383 cm

from a strong admixture of different triplet and singlet states. Their state composition, as obtained from our DMRG-SCF-SO calculations, compares nicely with the CASSCF-SO data of Gendron et al.110 . ACS Paragon27 Plus Environment

Journal of Chemical Theory and Computation

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

Page 28 of 37

TABLE II. Relative energies ∆E [cm−1 ], g-factors (absolute values), and state composition of the 2+ two lowest-energy electronic states of NpO2+ 2 as well as the lowest electronic states of PuO2 as

obtained from DMRG-SCF-SO and DMRG/NEVPT2-SO calculations. The same active orbital space was employed as for the CASSCF-SO/CASPT2-SO data taken from the work of Gendron et al.110 . For details on the 4c-MRCI calculation, see text. NpO2+ 2 DMRG-SCF-SO NEVPT2-SO CASPT2-SOa

4c-MRCI

ground state

∆E

0

0

0

0



gk

4.228

4.235

4.233

4.283

gb⊥

0.002

0.002

0.002

0.000

composition

89% φ, 11% δ

88% φ, 12% δ 88% φ, 12% δ

-

1st excited state ∆E

3294

3086

3107

-

2∆

gk

2.025

2.035

2.037

-

gb⊥

0.002

0.017

0.005

-

composition

99% δ, 1% π

98% δ, 2% π 98% δ, 2% π

5/2u

3/2u

-

PuO2+ 2 DMRG-SCF-SO ground state

∆E

4g

gk

6.076

6.076

g⊥

0.000

0.000

composition 1st excited state ∆E 0+ g

composition

2nd excited state ∆E 1g

composition

3rd excited state ∆E 5g

composition

0

CASSCF-SOa

96% 3 Hg , 2% 1 Γg 4383

0

96% 3 Hg 4388

3 1 + 3 − 3 1 + 57% 3 Σ− g , 27% Πg , 14% Σg 57% Σg , 27% Πg , 14% Σg

6985

6985

3 1 3 − 3 1 36% 3 Σ− g , 48% Πg , 14% Πg 35% Σg , 50% Πg , 14% Πg

7150 99% 3 Hg

7152 99% 3 Hg

a

2+ Ref. 110; CAS(7,10)PT2-SO data for NpO2+ 2 and CAS(8,10)SCF-SO data for PuO2 .

b

A non-zero value for g⊥ indicates a minor symmetry-breaking in the optimized spin-free wave

functions due to the neglect of point group symmetry. ACS Paragon28 Plus Environment

Page 29 of 37

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

Journal of Chemical Theory and Computation

VI.

CONCLUSIONS AND OUTLOOK

In this work, we presented a state-interaction (SI) approach for nonorthogonal matrix product state (MPS) wave functions dubbed as MPS-SI where the MPSs can be obtained from several different (state-specific) DMRG-SCF wave function optimizations. Our MPSSI approach complements existing implementations for traditional CI-type wave functions while generalizing earlier SI approaches for DMRG wave functions to nonorthogonal wave functions as basis set for the MPS-SI approach. Allowing for nonorthogonality of the input MPS wave functions adds a greater flexibility to the preceeding orbital optimization step, which will be particularly beneficial for transition metal and lanthanide/actinide compounds, where state-specific optimizations allow us to explicitly consider different open-shell d- or f -orbital occupations. Moreover, our nonorthogonal MPS-SI approach no longer requires a state averaging of molecular states with different spin symmetry at the orbital optimization step if, for example, spin-orbit coupling elements between these states are to be calculated. In this work, we discussed the transformation of orbitals and MPS wave functions from a nonorthogonal to a biorthonormal basis by formulating a nonunitary orbital transformation algorithm for MPS wave functions. This opens up for an efficient calculation of overlap matrix elements between nonorthogonal MPS wave functions and matrix elements of the one- and two-particle reduced transition density matrix as well as of the spin-orbit operator. After diagonalization of the resulting spin-free and spin-orbit Hamiltonian matrices target properties are calculated in the basis of the corresponding orthogonal and non-interacting (spin-orbit coupled) eigenstates. We demonstrated the applicability of the nonorthogonal MPS-SI approach for two exam2+ 1 2 ple open-shell f n actinide molecules, namely NpO2+ 2 (f complex) and PuO2 (f complex)

for which we calculated ground- and excited-state properties including g-factors. Our spinorbit DMRG-SCF/NEVPT2 g-factors for the two lowest doublet states of NpO2+ 2 agree well with corresponding CASSCF/CASPT2 data reported in the literature. The ground-state g-factors are close to the values obtained from a four-component MRCI g-tensor calculation which includes spin-orbit coupling variationally from the outset. This indicates that the property values are (close to being) converged with respect to the number of spin-free states included in our MPS-SI approach. Subsequently, we carried out DMRG-SCF g-factor calculations for PuO2+ 2 by allowing the lowest spin-free singlet and triplet states to interact ACS Paragon29 Plus Environment

Journal of Chemical Theory and Computation

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

Page 30 of 37

through spin-orbit coupling. We found a similar agreement of our DMRG-SCF results to the CASSCF values as was the case for NpO2+ 2 which illustrates that our nonorthogonal MPS-SI approach works equally well for cases where the interacting states do not share the same molecular orbital basis. Combined with our recently developed Cholesky-decomposition DMRG-NEVPT2 implementation, our MPS-SI approach will be valuable not only for the study of (magnetic) properties of transition metal and/or heavy-element complexes in different spin states but also for the full exploration of the photophysics and, more generally, excited-state (surface hopping) dynamics of chromophores and light-harvesting materials. The latter requires, besides the calculation of non-adiabatic coupling elements to locate and characterize conical intersections, the determination of intersystem crossing probabilities between electronic states of different spin multiplicity induced by spin-orbit coupling.

ACKNOWLEDGMENTS This work was supported by the Schweizer Nationalfonds. SK is grateful to Prof. H. J. Aa. Jensen (SDU Odense) for providing him access to the four-component MRCI-EPR module. JA acknowledges support from the U.S. Department of Energy, Office of Basic Energy Sciences, Heavy Element Chemistry program, under grant DE-SC0001136 (formerly DE-FG02-09ER16066). JA thanks ETH Zurich for providing assistance for his stay as a visiting professor.

APPENDIX The action of an annihilation operator in a nonorthogonal basis An annihilation operator in a nonorthogonal spin-orbital basis acting on an ONV is considered here. We follow closely the derivation that can be found in Ref. 59. Assume an ONV for L spin-orbitals as it was given in Eq. (12). Recalling that annihilation of the vacuum state is impossible, aqσ |vaci =

X p

1/2 Sqp a ˜pσ |vaci = 0 , | {z } =0 ∀p

ACS Paragon30 Plus Environment

(A1)

Page 31 of 37

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

Journal of Chemical Theory and Computation

we write the action of an annihilation operator on the ONV as (anti-)commutator of the annihilation operator and the product of even (odd) numbers of creation operators   !  !  L L Y Y     k k pσ pσ aqσ  ∓ a†pσ a†pσ aqσ |ki =  aqσ  |vaci   p=1 p=1 | {z } =0   L Q † kpσ   aqσ , apσ |vaci for even L    p=1  = .     L   Q † kpσ   |vaci for odd L apσ  aqσ ,

(A2)

p=1

The (anti-)commutator appearing in Eq. (A2) can be written in terms of linear combinations of a basic anticommutator of an annihilation and an creation operator exploiting the corresponding operator identity i h ˆ ˆ ˆ A, B1 · · · BL n

    

o   ˆ ˆ ˆ A, B1 · · · BL 

=

L X p=1

o n ˆ Bˆp · · · BˆL . (−1)p−1 Bˆ1 · · · A,

(A3)

Inserting Eq. (A3) into the right-hand side of Eq. (A2) yields aqσ |ki =

L X

kpσ (−1)

"

p−1 P

kjσ

j=1

p=1

#



a†1σ

k1σ

 kLσ n kpσ o · · · a†Lσ |vaci . · · · aqσ , a†pσ | {z }

(A4)

=Sqp δσσ

Resolving the anticommutator by the identity in Eq. (10), as indicated in Eq. (A4), yields the final result aqσ |ki =

L X

kpσ (−1)

"

p−1 P j=1

kjσ

#

Sqp |k1 k2 . . . 0p . . . kL i .

(A5)

p=1

REFERENCES 1

Barbatti, M. WIREs Comput. Mol. Sci. 2011, 1, 620–633.

2

Mai, S.; Marquetand, P.; Gonz´alez, L. Int. J. Quant. Chem. 2015, 115, 1215–1231.

3

Suaud, N.; Bonnet, M.-L.; Boilleau, C.; Lab`eguerie, P.; Guih´ery, N. J. Am. Chem. Soc. 2009, 131, 715–722.

4

Chang, J.; Fedro, A. J.; van Veenendaal, M. Phys. Rev. B 2010, 82, 075124. ACS Paragon31 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

5

Marian, C. M. WIREs Comp. Mol. Sci. 2012, 2, 187–203.

6

Kaupp, M., B¨ uhl, M., Malkin, V. G., Eds. Calculation of NMR and EPR Parameters: Theory and Applications; Wiley-VCH, 2004.

7

Gerloch, M.; McMeeking, R. F. J. Chem. Soc., Dalton Trans. 1975, 2443–2451.

8

Bolvin, H.; Autschbach, J. In Handbook of Relativistic Quantum Chemistry; Liu, W., Ed.; Springer: Berlin, 2016.

9

Bolvin, H. ChemPhysChem 2006, 7, 1575–1589.

10

Ganyushin, D.; Neese, F. J. Chem. Phys. 2013, 138, 104113.

11

Vad, M. S.; Pedersen, M. N.; Nørager, A.; Jensen, H. J. A. J. Chem. Phys. 2013, 138, 214106.

12

Sharkas, K.; Pritchard, B.; Autschbach, J. J. Chem. Theory Comput. 2015, 11, 538–549.

13

van Lenthe, E.; Wormer, P. E. S.; van der Avoird, A. J. Chem. Phys. 1997, 107, 2488– 2498.

14

Repisky, M.; Komorovsky, S.; Malkin, E.; Malkina, O. L.; Malkin, V. G. Chem. Phys. Lett. 2010, 488, 94–97.

15

Verma, P.; Autschbach, J. J. Chem. Theory Comp. 2013, 9, 1932–1948.

16

Lushington, G. H.; Grein, F. Int. J. Quant. Chem. 1996, 60, 1679–1684.

17

Brownridge, S.; Grein, F.; Tatchen, J.; Kleinschmidt, M.; Marian, C. M. J. Chem. Phys. 2003, 118, 9552–9562.

18

Neese, F. J. Chem. Phys. 2005, 122, 034107.

19

Vancoillie, S.; Malmqvist, P.-˚ A.; Pierloot, K. ChemPhysChem 2007, 8, 1803–1815.

20

Neese, F. Mol. Phys. 2007, 105, 2507–2514.

21

Tatchen, J.; Kleinschmidt, M.; Marian, C. M. J. Chem. Phys. 2009, 130, 154106.

22

Roemelt, M. J. Chem. Phys. 2015, 143, 044112.

23

Sayfutyarova, E. R.; Chan, G. K.-L. J. Chem. Phys. 2016, 144, 234301.

24

Szalay, P. G.; M¨ uller, T.; Gidofalvi, G.; Lischka, H.; Shepard, R. Chem. Rev. 2012, 112, 108–181.

25

Roca-Sanjuan, D.; Aquilante, F.; Lindh, R. WIREs Comp. Mol. Sci. 2012, 2, 585–603.

26

Stein, C. J.; von Burg, V.; Reiher, M. J. Chem. Theory Comput. 2016, 12, 3764–3773.

27

Olsen, J. Int. J. Quantum Chem. 2011, 111, 3267–3272.

28

Stein, C. J.; Reiher, M. J. Chem. Theory Comput. 2016, 12, 1760–1771.

29

Aquilante, F. et al. J. Comput. Chem. 2015, 37, 506–541. ACS Paragon32 Plus Environment

Page 32 of 37

Page 33 of 37

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

Journal of Chemical Theory and Computation

30

White, S. R. Phys. Rev. Lett. 1992, 69, 2863–2866.

31

White, S. R. Phys. Rev. B 1993, 48, 10345–10356.

32

Schollw¨ock, U. Rev. Mod. Phys. 2005, 77, 259–315.

33

Schollw¨ock, U. Ann. Phys. 2011, 326, 96–192. 34 ¨ Noack, R. M.; Solyom, J.; Tincani, L. Applications of quantum information in Legeza, O.; the density-matrix renormalization group. Computational Many-Particle Physics. 2008; pp 653–664.

35

Chan, G. K.-L.; Dorando, J. J.; Ghosh, D.; Hachmann, J.; Neuscamman, E.; Wang, H.; Yanai, T. In Frontiers in Quantum Systems in Chemistry and Physics; Wilson, S., Grout, P. J., Maruani, J., Delgado-Barrio, G., Piecuch, P., Eds.; Progress in Theoretical Chemistry and Physics; Springer Netherlands, 2008; pp 49–65.

36

Marti, K. H.; Reiher, M. Z. Phys. Chem. 2010, 224, 583–599.

37

Marti, K. H.; Reiher, M. Phys. Chem. Chem. Phys. 2011, 13, 6750–6759.

38

Chan, G. K.-L.; Sharma, S. Annu. Rev. Phys. Chem. 2011, 62, 465–481.

39

Wouters, S.; van Neck, D. Eur. Phys. J. D 2014, 68, 272.

40

Kurashige, Y. Mol. Phys. 2014, 112, 1485–1494.

41

Yanai, T.; Kurashige, Y.; Mizukami, W.; Chalupsky, J.; Lan, T. N.; Saitow, M. Int. J. Quantum Chem. 2015, 115, 283–299.

42

¨ Szalay, S.; Pfeffer, M.; Murg, V.; Barcza, G.; Verstraete, F.; Schneider, R.; Legeza, O. Int. J. Quantum Chem. 2015, 115, 1342–1391.

43

Knecht, S.; Hedeg˚ ard, E. D.; Keller, S.; Kovyrshin, A.; Ma, Y.; Muolo, A.; Stein, C. J.; Reiher, M. Chimia 2016, 70, 244–251.

44

Chan, G. K.-L.; Keselman, A.; Nakatani, N.; Li, Z.; White, S. R. J. Chem. Phys. 2016, 145, 014102.

45

Zgid, D.; Nooijen, M. J. Chem. Phys. 2008, 128, 144116.

46

Ghosh, D.; Hachmann, J.; Yanai, T.; Chan, G. K.-L. J. Chem. Phys. 2008, 128, 144117.

47

Knecht, S.; Legeza, O.; Reiher, M. J. Chem. Phys. 2014, 140, 041101.

48

Malmqvist, P.-A.; Roos, B. O. Chem. Phys. Lett. 1989, 155, 189–194.

49

Malmqvist, P.-A. Int. J. Quantum Chem. 1986, 30, 479–494.

50

Olsen, J.; Godefroid, M. R.; J¨onsson, P.; Malmqvist, P.-A.; Fischer, C. F. Phys. Rev. E 1995, 52, 4499–4508.

51

Olsen, J. J. Chem. Phys. 2015, 143, 114102. ACS Paragon33 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

52

Chen, Z.; Chen, X.; Wu, W. J. Chem. Phys. 2013, 138, 164120.

53

Olsen, J.; Roos, B. O.; Jørgensen, P.; Jensen, H. J. A. J. Chem. Phys. 1988, 89, 2185– 2192.

54

Roos, B. O. The Multiconfigurational Self-Consistent Field Theory. Lecture Notes in Quantum Chemistry: European Summer School in Quantum Chemistry. Berlin Heidelberg, 1992; pp 177–254.

55

Roos, B. O. In Radiation Induced Molecular Phenomena in Nucleic Acids; Shukla, M. K., Leszczynski, J., Eds.; Springer, 2008; pp 125–156.

56

Keller, S.; Dolfi, M.; Troyer, M.; Reiher, M. J. Chem. Phys. 2015, 143, 244118.

57

Keller, S.; Reiher, M. J. Chem. Phys. 2016, 144, 134101.

58

Moshinsky, M.; Seligman, T. H. Ann. Phys. 1971, 66, 311–334.

59

Helgaker, T., Jørgensen, P., Olsen, J., Eds. Molecular Electronic-Structure Theory, 1st ed.; Wiley & Sons, Chichester (England), 2000.

60

L¨owdin, P.-O. J. Chem. Phys. 1950, 18, 365–375.

61

Heitler, W.; London, F. Z. Phys. 1927, 44, 455–472.

62

L¨owdin, P. J. Mol. Struct. (Theochem) 1991, 229, 1–14.

63

Wu, W.; Su, P.; Shaik, S.; Hiberty, P. C. Chem. Rev. 2011, 111, 7557–7593.

64

Rashid, Z.; van Lenthe, J. H. J. Chem. Phys. 2013, 138, 054105.

65

Chen, Z.; Chen, X.; Wu, W. J. Chem. Phys. 2013, 138, 164119.

66

Thom, A. J. W.; Head-Gordon, M. J. Chem. Phys. 2009, 131, 124113.

67

Sundstrom, E. J.; Head-Gordon, M. J. Chem. Phys. 2014, 140, 114103.

68

L¨owdin, P.-O. Phys. Rev. 1955, 97, 1474–1489.

69

Prosser, F.; Hagstrom, S. Int. J. Quantum Chem. 1968, 2, 89–99.

70

Weltin, E. E. Int. J. Quantum Chem. 1976, 10, 163–174.

71

Dalgaard, E. Theoret. Chim. Acta 1983, 64, 181–186.

72

Verbeek, J.; van Lenthe, J. H. Int. J. Quantum Chem. 1991, 40, 201–210.

73

Verbeek, J.; van Lenthe, J. H. J. Mol. Struct: Theochem 1991, 229, 115–137.

74

Koch, H.; Dalgaard, E. Chem. Phys. Lett. 1993, 212, 193–200.

75

Amovilli, C. In Quantum Systems in Quantum Chemistry and Physics; McWeeny, R., Maruani, J., Smeyers, Y. G., Wilson, S., Eds.; Kluwer Academics Publishers, 1997; pp 343–347.

76

Dijkstra, F.; van Lenthe, J. H. Int. J. Quantum Chem. 1998, 67, 77–83. ACS Paragon34 Plus Environment

Page 34 of 37

Page 35 of 37

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

Journal of Chemical Theory and Computation

77

Slater, J. Phys. Rev. 1929, 34, 1293–1322.

78

Condon, E. Phys. Rev. 1930, 36, 1121–1133.

79

McDouall, J. J. W. Theor. Chim. Acta 1992, 83, 339–350.

80

Malmqvist, P.-A.; Roos, B. O.; Schimmelpfennig, B. Chem. Phys. Lett. 2002, 357, 230– 240.

81

Jordan, P.; Wigner, E. Z. Phys. 1928, 47, 631–651.

82

Legeza, O.; S´olyom, J. Phys. Rev. B 2003, 68, 195116.

83

Rissler, J.; Noack, R. M.; White, S. R. Chem. Phys. 2006, 323, 519–531.

84

Legeza, O.; S´olyom, J. Phys. Rev. Lett. 2006, 96, 116401.

85

Rose, M. E. Elementary Theory of Angular Momentum; Dover Publications, 1995.

86

Aquilante, F. et al. J. Comp. Chem. 2015, 37, 506–541.

87

Douglas, M.; Kroll, N. M. Ann. Phys. 1974, 82, 89–155.

88

Hess, B. A. Phys. Rev. A 1986, 33, 3742–3748.

89

Reiher, M. Theor. Chem. Acc. 2006, 116, 241–252.

90

H. J. Aa. Jensen, oral presentation given at the REHE conference in M¨ uhlheim (Germany), 2005.

91

Kutzelnigg, W.; Liu, W. J. Chem. Phys. 2005, 123, 241102.

92

Saue, T. ChemPhysChem 2011, 12, 3077–3094.

93

Peng, D.; Reiher, M. Theor. Chem. Acc. 2012, 131, 3–20.

94

Heß, B. A.; Marian, C. M.; Wahlgren, U.; Gropen, O. Chem. Phys. Lett. 1996, 251, 365.

95

Tatchen, J.; Marian, C. M. Chem. Phys. Lett. 1999, 313, 351–357.

96

B. Schimmelpfennig, AMFI is an atomic mean-field spin-orbit integral program, University of Stockholm, 1996.

97

Kurashige, Y.; Yanai, T. J. Chem. Phys. 2011, 135, 094104.

98

Kurashige, Y.; Chalupsky, J.; Lan, T. N.; Yanai, T. J. Chem. Phys. 2014, 141, 174111.

99

Wouters, S.; Van Speybroeck, V.; Van Neck, D. J. Chem. Phys. 2016, 145, 054120.

100

Sharma, S.; Chan, G. K.-L. J. Chem. Phys. 2014, 141, 111101.

101

Sokolov, A. Y.; Chan, G. K.-L. J. Chem. Phys. 2016, 144, 064102.

102

Roemelt, M.; Guo, S.; Chan, G. K.-L. J. Chem. Phys. 2016, 144, 204113.

103

Guo, S.; Watson, M. A.; Hu, W.; Sun, Q.; Chan, G. K.-L. J. Chem. Theory Comput. 2016, 12, 1583–1591.

ACS Paragon35 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

104

Freitag, L.; Knecht, S.; Angeli, C.; Reiher, M. J. Chem. Theory Comput. 2016, submitted, arXiv:1608.02006.

105

Rota, J.-B.; Knecht, S.; Fleig, T.; Ganyushin, D.; Saue, T.; Neese, F.; Bolvin, H. J. Chem. Phys. 2011, 135, 114106.

106

Roos, B. O.; Lindh, R.; Malmqvist, P.-˚ A.; Veryazov, V.; Widmark, P.-O. J. Phys. Chem. A 2004, 108, 2851.

107

Roos, B. O.; Lindh, R.; Malmqvist, P.-˚ A.; Veryazov, V.; Widmark, P.-O. Chem. Phys. Lett. 2005, 409, 295–299.

108

Wolf, A.; Reiher, M.; Hess, B. A. J. Chem. Phys. 2002, 117, 9215–9226.

109

Huber, K. P.; Herzberg, G. Molecular Spectra and Molecular Structure: Constants of Diatomic Molecules; Van Nostrand Reinhold, New York, 1979.

110

Gendron, F.; Pritchard, B.; Bolvin, H.; Autschbach, J. Inorg. Chem. 2014, 53, 8577–8592.

111

Dunning Jr., T. H. J. Chem. Phys. 1989, 90, 1007.

112

Dyall, K. G. Theor. Chem. Acc. 2007, 117, 491.

ACS Paragon36 Plus Environment

Page 36 of 37

Page 37 of 37

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 ACS Paragon Plus Environment