Combining Internally Contracted States and Matrix Product States To

Subscriber access provided by Vanderbilt Libraries

Article

Combining internally contracted states and matrix product states to perform multireference perturbation theory Sandeep Sharma, Gerald Knizia, Sheng Guo, and Ali Alavi J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00898 • Publication Date (Web): 06 Jan 2017 Downloaded from http://pubs.acs.org on January 9, 2017

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

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

Page 1 of 14

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

Combining internally contracted states and matrix product states to perform multireference perturbation theory Sandeep Sharma∗ Max Planck Institute for Solid State Research, Heisenbergstraße 1, 70569 Stuttgart, Germany Department of Chemistry and Biochemistry, University of Colorado Boulder, Boulder, CO 80302, USA

Gerald Knizia Department of Chemistry, Pennsylvania State University, 401A Chemistry Bldg, University Park, PA 16802, USA

Sheng Guo Department of Chemistry, Princeton University, Princeton, NJ 08544, USA

Ali Alavi† Max Planck Institute for Solid State Research, Heisenbergstraße 1, 70569 Stuttgart, Germany Deptartment of Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom We present two efficient and intruder-free methods for treating dynamic correlation on top of general multi-configuration reference wave functions—including such as obtained by the density matrix renormalization group (DMRG) with large active spaces. The new methods are the second order variant of the recently proposed multi-reference linearized coupled cluster method (MRLCC) [S. Sharma, A. Alavi, J. Chem. Phys. 143, 102815 (2015)], and of N-electron valence perturbation theory (NEVPT2), with expected accuracies similar to MRCI+Q and (at least) CASPT2, respectively. Great efficiency gains are realized by representing the first-order wave function with a combination of internal contraction (IC) and matrix product state perturbation theory (MPSPT). With this combination, only third order reduced density matrices (RDMs) are required. Thus, we obviate the need for calculating (or estimating) RDMs of fourth or higher order; these had so far posed a severe bottleneck for dynamic correlation treatments involving the large active spaces accessible to DMRG. Using several benchmark systems, including first and second row containing small molecules, Cr2 , pentacene and oxo-Mn(Salen), we show that active spaces containing at least 30 orbitals can be treated using this method. On a single node, MRLCC2 and NEVPT2 calculations can be performed with over 550 and 1100 virtual orbitals, respectively. We also critically examine the errors incurred due to the three sources of errors introduced in the present implementation calculating second order instead of third order energy corrections, use of internal contraction and approximations made in the reference wavefunction due to DMRG. I.

INTRODUCTION

Many processes in chemistry and biology involve catalysis by transition metal centers. These include systems ranging from ubiquitous catalysts of organic synthesis to the active sites of most enzymes. Unfortunately, the electronic structure of such systems is often strongly correlated, which makes the theoretical study of their nature and reactivity difficult. Strong correlation is characterized by the necessity of considering more than one Slater determinant to describe the electronic wavefunction even qualitatively. This property renders standard methods like Kohn-Sham Density Functional Theory (DFT) or single-reference wave function methods like MP2 or CCSD(T) unreliable, if not straight out inapplicable. To describe such systems, multi-reference (MR) methods are needed. However, currently available MR methods still have various defects, which frequently force researchers into uneasy compromises regarding active space size, quantitative accuracy, reliability, intruder states, or the need to strongly truncate the studied sys-

tem to make computational expense manageable. We here present new variants of multi-reference linearized coupled cluster (MRLCC2)1–3 and N -electron valence perturbation theory (NEVPT2),4–7 which greatly alleviate these problems. By judiciously combining internal contraction (IC)8–12 and matrix product state perturbation theory (MPSPT),1–3 we arrive at MR methods capable of efficiently and accurately describing large and complex strongly correlated systems quantitatively (i.e., including dynamic correlation, vide infra) if combined with a suitable active space treatment, such as provided by DMRG reference functions. To quantitatively describe a strongly correlated system, both dynamic correlation and static correlation must be accounted for. However, both problems have rather different demands on the wave function ansatz. To account for static correlation, a wave function ansatz needs high flexibility in many-body configuration space, but only low flexibility in one particle space—a small number of active orbitals is normally sufficient. This problem is commonly treated with brute force full config-

ACS Paragon Plus Environment


Page 2 of 14 2

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

uration interaction (FCI),13 which is practical for up to approximately 16 active orbitals. However, recently significant progress has been made with alternative methods to FCI; in particular, due to developments in the density matrix renormalization group (DMRG) algorithm, FCI Quantum Monte Carlo (FCIQMC),14–16 and selected configuration interaction approaches17–19 , static correlation can now be treated accurately in much larger active spaces. We here focus on DMRG,20,21 a variational method which minimizes the energy of a wavefunction parametrized as a matrix product state (MPS).22,23 DMRG can handle active spaces of around 30–40 orbitals, and in some cases even up to 10024–34 . However, by themselves the mentioned methods are not efficient for obtaining quantitative accuracy and for this dynamic correlation must also be calculated. To account for dynamic correlation, a wave function ansatz capable of efficiently handling a large set of (virtual) orbitals is required, while the effective many-body structure of the wave function can be comparatively simple and is amendable to convenient parameterization. An efficient approach for this is to treat excitations out of a qualitatively correct, strongly correlated reference function. Methods of this class are hierarchies based on truncated multireference configuration interaction (MRCI)10,35,36 , various flavors of perturbation theory6,37–39 and coupled cluster theory.40,41 However, the established methods of these classes are based on CI treatments of the static correlation, which either makes them cumbersome to apply (RAS, GAS)42–44 , or limits them to small active spaces (CAS). Some alternative approaches fall outside this static-dynamic partitioning framework (e.g., methods based on geminals45,46 , Jastrow factors47 or reduced density matrices48 ); while promising, they have so far not matured into routinely applied methods. For these reasons, a major challenge in recent years has been finding ways of combining the powerful DMRG approach to static correlation with one of the techniques more suitable for dynamic correlation (Sec. II). It was only recently realized that second-order perturbation theories (PT2s) can be formulated with matrix product states (MPS) by recasting them into an effective optimization problem.1–3 This approach, called matrix product state perturbation theory (MPSPT), converges to the fully uncontracted version of a PT2, and can be performed with relatively minor modifications to an existing DMRG program.1–3 MPSPT can treat complex zeˆ 0 involving two-body interacroth order Hamiltonians H tions; this allows formulating PT2 methods which are ˆ 0 -based PT2s such much more accurate than one-body-H as MP2 or CASPT2. We recently introduced multireference linearized coupled cluster (MRLCC) as one such method.1,3 While formally a PT, the third order MRLCC3 was shown to be extremely accurate (comparable to MRCI+Q), stable with no intruder state problems, exactly size-consistent, and to a large extent able to compensate for sub-optimal reference states.1,3

However, the previous implementation of MRLCC is computationally still rather expensive; this is a result of parameterizing the perturbed state using an MPS, which cannot efficiently treat large virtual orbital spaces. Here we solve this problem by parameterizing the perturbed wave function using a mixture of Celani-Werner-style partial internal contraction (IC)11,12 and MPSPT; and limiting ourselves to calculating only the second order correction to the energies. The resulting method can effectively handle large virtual spaces, while at the same time being applicable to systems with 30 active orbitals and more. As a side product, we also obtain an efficient variant of NEVPT2, which can be applied to systems with 1000 virtual orbitals and 30 active orbitals at the same time without the need for invoking the strongcontraction approximation. The article is organized as follows: Sec. II discusses MR methods of describing dynamic correlation, including their benefits and drawbacks, and previous work aiming to make them compatible with DMRG. Sec. III introduces the NEVPT2 and MRLCC2 theories, and covers our approach to their efficient implementation with IC and MPSPT. Sec. IV covers accuracy tests and benchmark applications; in particular, applications to first row dimers, Cr2 , pentacene and oxo-Mn(Salen), are shown, and timing tests for MRLCC2 and NEVPT2 with different active spaces and basis set sizes are presented. Sec. V covers concluding remarks and outlook for future research.

II. PARAMETERIZING DYNAMIC CORRELATION: PROBLEMS AND SOLUTIONS

The naive treatment of dynamical correlation in MR problems is extremely challenging,49,50 but a major breakthrough was made with the introduction of internal contraction (IC) by Meyer.8 IC, now widely used in MR theories7,9,41,51–53 aggressively truncates the many body basis and reduces the cost of the calculation from exponential to merely polynomial in the number of active space orbitals. The difficulty with IC is that the working equations become extremely cumbersome and can often only be derived using domain specific computer algebra systems. The working equations can be simplified by treating the active orbitals and the core orbitals on an equal footing as was done in the MRCI program developed by Werner-Knowles (WK scheme)9,10 . But large efficiency gains can be achieved by recognizing that density matrices should only contain active space indices. This was done in the case of CASPT2 by Celani and Werner (CW scheme)11,12 , and much later for MRCI by Shamasundar et al36 . The straightforward use of IC in perturbation theory or truncated configuration interaction with DMRG is complicated by the fact that fourth (and sometimes even higher) order reduced density matrices (RDM) in the active space are required. Calculation and storage of


Page 3 of 14

Journal of Chemical Theory and Computation 3

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

such high order RDMs is limited to an active space of around 25 orbitals beyond which it become prohibitively expensive54–57 . To circumvent this difficulty, some researches have resorted to approximating the higher body RDM by reconstructing its disconnected part by antisymmetric multiplication of lower body RDM and ignoring the density cumulant which is the connected part53,58–62 . The RDM reconstruction can be performed by setting the three body and four body cumulants to zero as was done for n-electron valence perturbation theory (NEVPT2)60 ; which resulted in severe numerical problems rendering the theory virtually unusable. Canonical Transformation (CT) theory53,58,59 also sets three and higher body cumulants to zero and is known to suffer from intruder state problems (the CT intruder states are different than the ones in perturbation theory). A much milder approximation is to use the exact three body cumulants but setting the four body cumulants to zero61,62 . This increases the numerical stability considerably but artifacts still remain that have to be eliminated by using level shifts. The underlying difficulty in neglecting cumulants beyond a certain rank is that, while the cumulants decay exponentially rapidly with the rank in weak correlation, they either decay very slowly or not at all in highly multireference situations63 . In this work we follow a different route, by using partial internal contraction. All terms that require 4-RDMs are treated fully uncontracted, and the rest of the terms, requiring three body or lower RDMs are treated using IC. This partitioning is in fact also used in the MRCI of Shamasundar et al.36 . The reason this is a practical route is as follows: all terms which are difficult to treat using IC and require 4-RDMs (IC states with three active space indices and one virtual or core index) are relatively easily treated using uncontracted theories, because they involve only one core or virtual orbital, and are effectively singly excited states. However, if these uncontracted terms are represented using wave functions, as done in previous work,11,12,36 their number still increases exponentially in the number of active orbitals. For this reason, we here treat such terms with MPSPT,1–3 which does not have this problem. This combination of IC and MPSPT is the key development of the current work, which makes efficient MRLCC2 and NEVPT2 methods accessible to the large active spaces treatable with DMRG. We will now discuss its theory and implementation.

III.

THEORY AND IMPLEMENTATION

orbitals (i.e., the union of the three subsets) are indexed by m, n, o, p. We work exclusively with spin free quantities, by integrating out the spin degrees of freedom. Thus, the Hamiltonian is mn ˆ op ˆ = tn E ˆm H m n + Wop Emn .

(1)

Here and in the following, repeated indices imply conˆ and traction, operators are written with a hat (e.g. E) numerical tensors are written without a hat (e.g. W ). tnm mn are the one electron integrals, Wop = hmn|opi are the two electron integrals, and X X ˆm = ˆ op = E a ˆ†mσ a E ˆnσ a ˆ†oσ a ˆ†pτ a ˆnτ a ˆmτ (2) n mn σ

στ

are the spin-free single- and double excitation operators, in which σ, τ ∈ {α, β} denote the spin degrees of freedom.

ˆ0 Definitions of H

In the Rayleigh-Schroedinger perturbation theory, the zeroth order wavefunction |Ψ0 i must be an eigenfunction ˆ 0 . However, for a given of the zeroth order Hamiltonian H ˆ 0, |Ψ0 i, this requirement does not uniquely determine H ˆ 0 will lead to a difand in general each viable choice of H ferent perturbation theory. For example, the CASSCF (or CAS-CI) wavefunction is, at the same time, an eigenfunction of the Fock operator, Dyall’s Hamiltonian,64 and ˆ 0 leads to Fink’s Hamiltonian;38,65 choosing these as H CASPT2 theory, NEVPT2 theory, and MRLCC2 theory, respectively, all of which have vastly different properties. In this work we focus on the Dyall- and Fink- Hamiltonians, which are described below. These have previously been shown to lead to accurate multi-reference perturbation methods which do not suffer from the intruder state problems plaguing CASPT2.66,67 Concretely, Dyall’s Hamiltonian64 is given by s ˆr rs ˆ tu b â ˆD = f j E î H i j + fr Es + Wtu Ers + fa Eb ,

(3)

where f is a Fock matrix. f contains modified one body integrals, approximating the interaction with core- and active electrons as a mean-field. It is defined as jk jk jr jr s Γr fij = tji + 2Wik − Wki + Wis − Wsi s b b ak ak ar ar fa = ta + 2Wbk − Wkb + Wbs − Wsb Γr rk rk frs = tsr + 2Wsk − Wsj ,

Notation

The molecular orbitals are divided into three subsets: the core orbitals indexed by i, j, k, l, which are always doubly occupied in the reference wavefunction |Ψ0 i (CASSCF or CAS-CI); the active orbitals r, s, t, u, which can have any occupancy in |Ψ0 i; and the virtual orbitals a, b, c, d, which are unoccupied in |Ψ0 i. General

ˆsr |Ψ0 i denotes the 1-RDM. When using where Γrs = hΨ0 |E CASSCF canonical orbitals, f becomes diagonal, and the diagonal elements can be viewed as orbital energies. A diagonal Fock matrix (implicitly assumed in the work of Angeli et al.4–7 ) simplifies NEVPT2 theory. However, in the present work we avoid this assumption, and the expressions derived here are more general.



Page 4 of 14 4

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

Fink’s Hamiltonian38,65 is mn ˆ op ˆ F = tn E ˆm H m n + Wop Emn

∆nex =0

,

(4)

where ∆nex = 0 signifies that all terms which would change the total number of electrons in the core, active or virtual spaces are removed. An interesting property of Fink’s Hamiltonian is that the zeroth order energy is equal to the energy of the reference wavefunction, and the first order energy is zero. This is not true in NEVPT2 theory; similarly to MøllerPlesset perturbation theory (MP2), in NEVPT2 the reference energy is given by the sum of the zeroth and first order energies, neither of which vanish. The first order wave function |Ψ1 i

In perturbation theory, the first order correction (Ψ1 ) to the reference wavefunction (Ψ0 ) is determined by solving the linear equation ˆ 0 − E0 |Ψ1 i = −Q ˆ Vˆ |Ψ0 i, H (5) ˆ −H ˆ 0 is the perturbation, and Q ˆ = in which Vˆ := H 1 − |Ψ0 ihΨ0 | is a projector. To define the representation of |Ψ1 i, we first split the full N -electron Fock-space into configuration spaces characterized by a unique occupation pattern (∆Nc , ∆Na , ∆Nv ). The three integers hereby denote the change in the total number of electrons in core (Nc ), active (Na ), and virtual (Nv ) orbitals relative to |Ψ0 i. With a two-body perturbation Vˆ , eight such spaces can be reached by Vˆ |Ψ0 i, and are therefore required in the expansion of |Ψ1 i (Eq. (5)). These spaces are listed in Table I, and their representation will be discussed next. Neither Dyall’s nor Fink’s Hamiltonian ˆ 0 ) change the number of electrons in the three orbital (H spaces; consequently, both are block-diagonal with respect to the eight configuration spaces. Therefore Eq. (5) can be solved in each space independently. To solve Eq. (5), we first choose a suitable basis of many-body states {Φµ }, covering the eight configuration classes. We then write the first order wave function as a linear combination |Ψ1 i = dµ |Φµ i. The choice of the basis states {Φµ } is not unique; various approaches are possible. For example, one could choose to include all individual determinants (or CSFs) which have non-zero overlap with Ψ1 ; this is known as an uncontracted formulation of the theory. Unfortunately, with this choice the cost of solving Eq. (5) generally scales exponentially with the number of active orbitals, quickly becoming impractical. Alternatively, the excitation classes can be represented using internal contraction (IC). Here the many-body basis {|Φµ i} is ˆµ formed by acting with orbital excitation operators O

TABLE I: The eight classes of perturber states reachable by applying a two-body perturbation Vˆ to the zeroth order state |Ψ0 i, and which therefore contribute to |Ψ1 i. The classes are formed by changing the total number of electrons in core, active, and virtual orbitals, relative to the occupation pattern of |Ψ0 i (indicated by positive/negative numbers). Column 5 and 6 denote the basis of the many-body states and the order of the active-space RDM needed, if the space were parameterized using internal contraction (IC). MPSPT denotes matrix product state perturbation theory. In the methods discussed here, the contributions to |Ψ1 i within spaces I–VIII can be solved for independently of each other.

I II III IV V VI

Occupation-change Core Active Virtual −2 0 +2 −1 −1 +2 −2 +1 +1 0 −2 +2 −2 +2 0 −1 0 +1

VII −1 VIII 0

+1 −1

0 +1

Basis in Needed Solution IC case RDM strategy a ˆb ˆ Ei Ej |Ψ0 i 1 IC îa E ˆrb |Ψ0 i E 2 IC îa E ˆjr |Ψ0 i E 2 IC ˆra E ˆsb |Ψ0 i E 3 IC r ˆs ˆ Ei Ej |Ψ0 i 3 IC îa E ˆsr |Ψ0 i, {E 3 IC ˆsa E îr |Ψ0 i, E îa |Ψ0 i} E ˆtr E ˆjs |Ψ0 i E 4∗ MPSPT r â ∗ ˆ Es Et |Ψ0 i 4 MPSPT

∗ In this work no 4-RDM is required because we do not use internal contraction for classes VII and VIII.

ˆb ˆµ = E âE (e.g. O i j ) onto the entire zeroth order state |Ψ0 i, rather than its individual configurations: ˆ µ |Ψ0 i. |Φµ i = O

(6)

A list of such states suitable for representing |Ψ1 i is shown in column 6 of Table I. With IC, the parameterization of |Ψ1 i becomes much more efficient than in uncontracted theories; the cost of solving Eq. (5) reduces from exponential to polynomial in the number of active orbitals. The drawback of this approach is that the IC basis states (6) are complicated, non-orthogonal manybody states. This leads to complex expressions for operator matrix elements between IC basis states; these generally contain series of tensor contractions involving high order RDMs, as shown in the next section.

Internal contraction

To demonstrate the positive and negative aspects of IC, we here derive the equations for perturber space IV (see Table I). We write both Ψ1 ’s and Vˆ |Ψ0 i’s space IV contribution as linear combination of IC states, rs ˆ a ˆ b |ΨIV 1 i = dab Er Es |Ψ0 i tu ˆ c ˆ d (Vˆ |Ψ0 i)IV = Wcd Et Eu |Ψ0 i.


(7) (8)

Page 5 of 14


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

tu Here dsr ab are the unknown parameters of |Ψ1 iIV , and Wcd are the two-body integrals. Substituting into Eq. 5, we get tu ˆ c ˆ d ˆ 0 − E0 drs E â ˆb H (9) ab r Es |Ψ0 i = −QWcd Et Eu |Ψ0 i.

To obtain closed numerical equations for dsr ab , we leftr ′ ˆ s′ ˆ project onto the bra states hΨ0 |Ea′ Eb′ , arriving at ˆ 0 − E0 ) E ˆra E ˆsb |Ψ0 idrs âr′′ E ˆbs′′ (H hΨ0 |E ab = r ′ ˆ s′ ˆ c ˆ d ˆ − hΨ0 |E ′ E ′ E E |Ψ0 iW tu . a

b

t

u

cd

With A and S as defined in Eqs. (12) and (11), this can be equivalently written in tensor form as

error-prone when done manually. For this reason, various computer programs have been developed to automate the evaluation process. The resulting expressions only contain one- and two-body integral tensors (t and W ), δ-symbols, and RDMs containing exclusively active indices. In this work we used the freely available program SecondQuanizationAlgebra, which was originally developed for implementing Canonical Transformation theory,68 and has since been modified and extended.61 For example, the program evaluates the left hand side of ˆ 0 being Fink’s Eq. 10, where A is given by (14) with H Hamiltonian, into a sum of eleven different tensor contraction terms. The first three of these are ′ ′ ′ ′

′ ′ ′ ′ a s b rasb Arrasb d

=−

r ′ a′ s′ b′ tcud Stcud w .

(10)

′ ′ ′ ′

a s b rasb cr a s b = Arrasb d ′ ′

The tensor S represents an overlap matrix of IC states r ′ a ′ s′ b ′ âr′′ E ˆbs′′ E ˆtc E ûd |Ψ0 i. Stcud = hΨ0 |E

′ ′ ′ ′ asb ˆ 0 − E0 ) E âE ˆb ˆ r′′ E ˆ s′′ (H Arrasb = hE r si a b âr′′ E ˆbs′′ (H ˆ 0 − E0 ), E ˆra E ˆsb i = hE âE ˆb ˆ ˆ s′′ E ˆ r′′ E + hE b r s (H0 − E0 )i a ′ ′ âr′ E ˆbs′ (H ˆ 0 − E0 ), E ˆra E ˆsb i = hE

′

(12) (13) (14)

To cancel the second term of (13), we used that |Ψ0 i ˆ 0 with eigenvalue E0 . For is an eigenfunction of H evaluating the commutator in Eq (14), only up to three-body RDMs are needed. Using similar procedures, IC equations for the other classes in Table I can be derived; the maximum order RDM required in these equations is shown in column 6. Since computing and storing a 4-RDM would be costly, we treat only the indicated six of the eight perturbation classes with IC.

′ ′

′

′

′ ′ ′

′

rasb r s t b + 4Wtta Γsrt δa + · · · , ′b d

(15)

The full equations for all six perturbation classes treated with IC, for both the Fink’s and Dyall’s zeroth order Hamiltonians, are given in the supplementary material. Once the equations have been derived, they have to be factorized into a sequence of binary tensor contractions (contractions in which one output tensor is formed by contracting the common indices of two input tensors). At this moment our program does not attempt to find the optimal factorization for a set of tensor contractions; while a powerful factorization algorithm has been reported in the literature,69 it is so far not applicable to our case, and not easy to implement. Instead we factorize each multi-tensor-contraction independently using a simple heuristic.? Several computer libraries capable of efficiently evaluating the final sequence of binary contractions have recently become available.70–73 We here use a prototyping library developed by one us, which will be described elsewhere. Solution of the linear equations In general, the IC basis functions (6) of a given excitation class are not orthogonal to each other. For example, for perturber class IV, the overlap matrix is r ′ a ′ s′ b ′ ˆ r′′ E ˆ s′′ E ˆcE ˆd Stcud = hE a b t ui ′ ′

Wick’s theorem and tensor contractions To arrive at practical working equations for the Ψ1 contributions in the various perturbation classes, the second-quantized operator expressions, such as âr′′ E ˆbs′′ (H0 − E0 ), E ˆra E ˆsb i hE

′

+ 4tab Γrsrs drasb δab

(11)

Using Wick’s theorem, this expression can be simplified into a sum of products of one 2-RDM with only active indices rstu, and two Kronecker-δs connecting the virtual indices abcd (see Eq. (16)). Similarly, A (Eq. (12)) can be expressed as a sum of tensor contractions via Wick’s theorem; these contain interaction tensors (W and t), δ-symbols with virtual and/or core indices, and up to four-body RDMs (with up to eight active indices—four ˆ 0 and four from |ΨIV i). However, the 4-RDM can from H 1 be eliminated, without approximation, as follows:

′

r s rasb a b = −4ttr Γts d δa δb

′

′

′ ′

′

′

r s a b sr b a = Γtu δc δd + Γtu δc δd .

(16)

An orthogonal basis can be generated by limiting the ˆtc E ûd |Ψ0 i to c ≥ d, and diagonalizing virtual indices in E ′ ′ the two body density matrix Γrtus (just the first term in Eq (16)). With this orthogonal basis, the metric S can be eliminated from Eq. (10), simplifying it to

of Eq. (14), must be transformed into a series of tensor contractions using Wicks’ theorem. While this processes is relatively straight-forward, it is extremely tedious and


′ ′ ′ ′ ′ ′ a b ˜αab A˜α =w ˜α a b , αab d

(17)


Page 6 of 14 6

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

where we have defined the orthogonal-basis quantities ′ ′ ′ ′ ′ ′ ′ α ′ a ′ b′ asb A˜αab =U α r s Arrasb U αrs d˜αab =U αrs drasb

w ˜ αab =U αrs wrasb .

(18)

Here U αrs is the αth eigenvector of the 2-RDM. We use the conjugate gradient (CG) method for solving Eq. (17). α′ a′ b′ ˜αab In each CG iteration, A˜αab d must be evaluated, which is non-trivial. Using Wick’s′ theorem, we only ara′ s′ b′ rasb rive at tensor contractions for Arrasb d in the nonorthogonal basis (note the missing ∼ s). Thus, we first convert d˜ to d by left multiplying by U † . Once d is calculated, c = Ad is evaluated in the non-orthogonal IC basis as described in the previous section. Finally, c is converted to c˜ by left multiplying by U . Thus, the entire CG scheme is performed in the orthogonal basis, and only the contraction c = Ad is evaluated in the non-orthogonal IC basis. A simpler solution scheme is possible if using Dyall’s Hamiltonian and canonical orbitals. In this case, the ′ ′ ′ ′ ′ ′ ′ ′ a s b rasb product cr a s b = Arrasb d can be simplified to ′ ′ r ′ s′ cr as b = (Krs + faa + fbb )drasb . Note that the virtual orbital indices are identical in tensors c and d. This simplification allows a closed-form solution of the linear equation, making an iterative CG scheme unnecessary. For this, only the K matrix must be diagonalized, which is of the same size as the overlap matrix S. Matrix product state perturbation theory

As noted in Table I, if we were to employ IC for perturber classes VII and VIII, this would require computing and storing the 4-RDM of the reference function |Ψ0 i. This would be costly, and would severely limit the applicability of our methods. To avoid this problem, we represent these perturber classes with matrix product state perturbation theory (MPSPT) instead of IC. The MPSPT approach requires no RDMs, and is capable of representing arbitrary uncontracted wave functions of the respective perturber classes. The core idea of this approach is to write the linear equation (5) as a minimization problem using the Hylleraas functional: H[Ψ1 ] = hΨ1 |H0 − E0 |Ψ1 i + 2hΨ1 |QV |Ψ0 i.

(19)

This reformulation makes Eq. (5) susceptible to the optimization techniques of matrix product state (MPS) ansatz wave functions. In this work both the reference wavefunction |Ψ0 i and the first order correction |Ψ1 i are parametrized using a MPS ansatz. In an MPS ansatz, the wavefunction is expressed as a product of a set of three-dimensional tensors, which is shown graphically in Fig. 1. The virtual bonds (horizontal lines in Fig. 1) encode entanglement between the different orbitals, and by increasing this bond dimension (M ), any wavefunction can be exactly represented

FIG. 1: A matrix product state (MPS) can be represented nl graphically using a series of three-dimensional tensors Ail−1 il . In these, the physical index nl (pointing upwards) denotes the occupation of the orbital l, and the other two indices il , il−1 (pointing horizontally, and known as virtual indices), are sequentially contracted.

TABLE II: Scaling of the computational cost for calculating ˆ 1 ) and two the overlap and expectation values of one body (O ˆ body (O2 ) operators between two MPS |Φ1 i and |Φ2 i. Here k is the number of orbitals, and M1 and M2 are the virtual bond dimensions of |Φ1 i and |Φ2 i respectively. Operation hΦ1 |Φ2 i ˆ 1 |Φ2 i hΦ1 |O ˆ 2 |Φ2 i hΦ1 |O

Computational Cost kM12 M2 + kM1 M22 k2 M12 M2 + k2 M1 M22 k3 M12 M2 + k3 M1 M22

using an MPS. An important property of the MPS ansatz is that overlaps and operator expectation values between two different MPS wave functions can be calculated at a polynomial cost in the number of orbitals (k) and the virtual bond dimension (M ). A summary of the costs of each such calculation is shown in Table II. The ability to perform the various linear algebra operations shown in Table II in polynomial time makes it possible to calculate the Hylleraas functional, and to subsequently optimize it by varying the tensors parameterizing the MPS |Ψ1 i. The optimization is performed using a sweep algorithm, just as DMRG optimizations of ground state wavefunctions. In a sweep iteration, each single tenl sor Anil−1 il of |Ψ1 i is optimized individually, while keeping other A tensors fixed. This is repeated until convergence of all As. So this sweep algorithm solves the optimization problem of a large, multi-linear function of k tensors by iteratively solving a single linear equation for a single tensor1–3 The simple MPSPT algorithm requires a sweep over the entire set of orbitals and optimizes each tensor separately. This can become prohibitively expensive when working with several hundred virtual orbitals. But by realizing that the first order state has a maximum of a single electron in the virtual space, one can once and for all fix the tensors corresponding to the virtual orbitals, and subsequently no sweeps over these orbitals are necessary. These virtual orbital tensors are shown in red in the figure 2; they are fixed at the beginning of the calculation and are not changed for the rest of the calculation. Note that the maximum virtual bond dimension required for these tensors is nv + 1, where nv is the number of virtual orbitals, which is less than a typical M of a few thousand that is routinely used in calculations. A similar simplification is possible for the perturber


Page 7 of 14


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


Page 8 of 14 8

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

TABLE III: Atomization energies (AE), ionization potential (IP) and electron affinities (EA) of various benchmark molecules are calculated using uncontracted MRLCC3 (UMRLCC3), uncontracted MRLCC2 (UMRLCC2) and partially internally contracted MRLCC2 implemented in the present work (p-ICMRLCC2). Where possible we have also tabulated accurate estimates of the FCI energies; AE were calculated using FCIQMC74 , whereas IP and EA were calculated using DMRG. The FCIQMC energies and geometries of molecules used for AE calculations are obtained from work by Cleland et al.74 and the geometries for molecules used for calculating IP and EA were obtained from work by Knizia et al.75 . All numbers reported here are converged to better than 1mEh accuracy. FCI Molecule DZ TZ QZ Atomization Energy C2 -0.207 -0.223 -0.228 N2 -0.320 -0.346 -0.356 F2 -0.044 -0.057 -0.059

-0.204 -0.220 -0.226 -0.317 -0.343 -0.354 -0.040 -0.052 -0.056

-0.204 -0.216 -0.221 -0.314 -0.337 -0.347 -0.051 -0.064 -0.067

-0.204 -0.216 -0.221 -0.314 -0.336 -0.346 -0.051 -0.063 -0.066

Ionization H2 O NH3 Cl2 OH

-0.674 -0.612 -0.411 -0.470

0.674 -0.612 -0.411 -0.452

Potential -0.674 -0.618 -0.449

Electron Affinity CH3 -0.048 CN 0.100 NO -0.053 SH -

DZ

UMRLCC3 TZ QZ

-

-

-0.672 -0.617 -0.412 -0.448

-0.684 -0.624 -0.415 -0.467

-

-

-0.051 -0.031 -0.021 0.099 0.124 0.132 -0.054 -0.029 -0.017 0.040 0.064 0.075

TABLE IV: Effect of internal contraction and accuracy of |Ψ0 i on absolute MRLCC2 energies. The rows show reference energy (E0 + 2099.0), and MRLCC2 correction due to classes I-VI of Tab. I(E2,IC ), and classes VII-VIII (E2,MPTSPT ), in atomic units. The UC (uncontracted) column contains reference energies obtained with uncontracted basis states for all classes. In the other four columns, an IC basis is used for perturber classes I-VI, and classes VII-VIII are represented with MPSPT. In these four columns, |Ψ0 i is an MPS with maximum virtual bond dimension M of 500, 80, 50 and 30 respectively. This illustrates the effect of approximating |Ψ0 i on the correlation energy. UC

E0 E2 (I–VI)a E2 (VII–VIII)b

Bond dimension M of |Ψ0 i 500 80 50 30 -0.2499 -0.2499 -0.2487 -0.2461 -0.2192 -1.0597 -1.0578 -1.0582 -1.0587 -1.0629 -0.0672 -0.0672 -0.0675 -0.0679 -0.0709

a

: E2 contrib. from I–VI, represented with IC

b

: E2 contrib. from VII–VIII, represented with MPSPT

The MPSPT part of the calculation dominates the overall cost of the NEVPT2 calculation, but is of similar cost (within a factor of 3) as internal contraction in the MRLCC2 calculations. Even though the cost of performing MPSPT is about the same in NEVPT2 calculations as in MRLCC2 calculations, the NEVPT2 cal-

-0.689 -0.628 -0.420 -0.474

DZ

UMRLCC2 TZ QZ

-0.684 -0.617 -0.412 -0.452

-0.689 -0.620 -0.416 -0.477

-0.049 -0.030 -0.021 0.100 0.124 0.133 -0.048 -0.022 -0.009 0.039 0.061 0.071

p-ICMRLCC2 DZ TZ QZ

-0.684 -0.617 -0.413 -0.470

-0.689 -0.620 -0.416 -0.470

-0.049 -0.031 -0.022 0.100 0.124 0.133 -0.048 -0.022 -0.010 0.039 0.061 0.071

culations can be stretched much further and performed more efficiently by using the fact that the perturbative contribution from each external orbital can be calculated independently of one another. We can perform simultaneous independent calculations where each one only correlates a subset of virtual orbitals. In the calculations on Pentacene with cc-pVTZ basis set and oxo-Mn(Salen) with cc-pVQZ basis set, there were respectively 722 and 1101 virtual orbitals, four independent MPSPT calculations were performed to treat the perturber class VII. It should be noted that, because 4-RDMs are not necessary large active spaces containing 30 orbitals (Cr2 with (12e,30o) active space) can be treated using MPSPT.

oxo-Mn(Salen)

Managanese salen derivatives like the Jacobson’s catalyst82,83 are used for performing enetioselective epoxiation of olefins. The mechanism of the epoxidation is not clear, with various experiments providing evidence supporting at least three different possibilities84,85 . To provide further clarification, several theoretical studies86–96 have tried to understand the electronic structure of these clusters. The chlorine containing neutral oxo-Mn(salen) model cluster studied in this work is shown in Figure 3. Based on DFT and coupled cluster calculations94 it was shown that the cluster can be in a singlet, triplet and


Page 9 of 14





Page 10 of 14

Page 11 of 14


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

although universal (much like Torheyden’s post F12 correction scheme106 ), has a severe memory bottleneck and can only be used for relatively small basis sets107 . A more efficient route for using F12-geminals in our method will be the approach used by Ten-no108 and later by Shiozaki et al.109,110 for multireference problems. PNO were used in? the multireference problems for the first time by Fink et al.111 , and after a long lull PNO based methods have seen a huge revival due to work of Neese et al.112,113 . Our next project will aim to make use of these recent developments by extending this method by combining it with F12-methods and PNOs.

edges the supported from University of Colorado. SG was supported via grants from the US National Science Foundation, NSF-CHE1650436 and NSF-CHE1657286.

Acknowledgments

Supporting information can be obtained free of charge from the ACS Publications website. It contains the tensor contractions needed to implement the working equations of the ICMRLCC and NEVPT2 theories.

The calculations made use of the facilities of the Max Planck Society’s Rechenzentrum Garching. S.S. acknowl-

∗ † 1

2

3 4

5

6

7

8

9

10

11 12 13

14

15

16

17 18

19

20 21 22

Electronic address: [email protected] Electronic address: [email protected] Sharma, S.; Jeanmairet, G.; Alavi, A. J. Chem. Phys. 2016, 144, 034103. Sharma, S.; Chan, G. K.-L. J. Chem. Phys. 2014, 141, 111101. Sharma, S.; Alavi, A. J. Chem. Phys. 2015, 143, 102815. Angeli, C.; Borini, S.; Cestari, M.; Cimiraglia, R. J. Chem. Phys. 2004, 121, 4043–9. Angeli, C.; Bories, B.; Cavallini, A.; Cimiraglia, R. J. Chem. Phys. 2006, 124, 054108. Angeli, C.; Cimiraglia, R.; Evangelisti, S.; Leininger, T.; Malrieu, J.-P. J. Chem. Phys. 2001, 114, 10252. Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. J. Chem. Phys. 2002, 117, 9138. Meyer, W. Modern Theoretical Chemistry; Plenum Press New York, 1977. Werner, H.-J.; Knowles, P. J. Theor. Chim. Acta. 1990, 78, 175–187. Knowles, P. J.; Werner, H.-J. Theor. Chim. Acta. 1992, 84, 95. Celani, P.; Werner, H.-J. J. Chem. Phys. 2000, 112, 5546. Celani, P.; Werner, H.-J. J. Chem. Phys. 2003, 119, 5044. Knowles, P.; Handy, N. Chem. Phys. Lett. 1984, 111, 315–321. Cleland, D.; Booth, G. H.; Alavi, A. J. Chem. Phys. 2010, 132, 041103. Booth, G. H.; Thom, A. J. W.; Alavi, A. J. Chem. Phys. 2009, 131, 054106. Booth, G. H.; Gruneis, A.; Kresse, G.; Alavi, A. Nature 2013, 493, 365–370. Huron, B. J. Chem. Phys. 1973, 58, 5745. Schriber, J. B.; Evangelista, F. A. J. Chem. Phys. 2016, 144, 161106. Holmes, A. A.; Tubman, N. M.; Umrigar, C. J. J. Chem. Theory Comput. 2016, White, S. R. Phys. Rev. Lett. 1992, 69, 2863. White, S. R. Phys. Rev. B 1993, 48, 10345. Schollw¨ ock, U. Ann. Phys. 2011, 326, 96.

VI.

23

24

25

26

27

28

29

30 31

32

33

34

35

36

37

38 39 40

41

42

SUPPORTING INFORMATION

Chan, G. K.-L.; Sharma, S. Ann. Rev. Phys. Chem. 2011, 62, 465. ¨ Noack, R. M.; S´ Legeza, O.; olyom, J.; Tincani, L. In Computational Many-Particle Physics; Fehske, H., Schneider, R., Weiße, A., Eds.; Lecture Notes in Physics; Springer Berlin Heidelberg, 2008; Vol. 739; pp 653–664. Marti, K. H.; Reiher, M. Phys. Chem. Chem. Phys. 2011, 13, 6750–6759. Kurashige, Y.; Yanai, T. J. Chem. Phys. 2009, 130, 234114. White, S. R.; Martin, R. L. J. Chem. Phys. 1999, 110, 4127. Chan, G. K. L.; Head-Gordon, M. J. Chem. Phys. 2002, 116, 4462. Wouters, S.; Poelmans, W.; Ayers, P. W.; Van Neck, D. Comp. Phys. Comm. 2014, 185, 1501–1514. Zgid, D.; Nooijen, M. J. Chem. Phys. 2008, 128, 014107. Nakatani, N.; Chan, G. K.-L. J. Chem. Phys. 2013, 138, 134113. Sharma, S.; Chan, G. K.-L. J. Chem. Phys. 2012, 136, 124121. Keller, S.; Dolfi, M.; Troyer, M.; Reiher, M. J. Chem. Phys. 2015, 143, 244118. Chan, G. K.-L.; Keselman, A.; Nakatani, N.; Li, Z.; White, S. R. J. Chem. Phys. 2016, 35. Werner, H.-J.; Reinsch, E. A. J. Chem. Phys. 1982, 76, 3144. Shamasundar, K. R.; Knizia, G.; Werner, H.-J. J. Chem. Phys. 2011, 135, 054101. Andersson, K.; Malmqvist, P. A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. J. Phys. Chem. 1990, 94, 5483– 5488. Fink, R. F. Chem. Phys. 2009, 356, 39–46. Hirao, K. Chem. Phys. Lett. 1992, 190, 374–380. Lyakh, D. I.; Musial, M.; Lotrich, V. F.; Bartlett, R. J. Chem. Rev. 2012, 112, 182–243. Evangelista, F. A.; Hanauer, M.; K¨ ohn, A.; Gauss, J. J. Chem. Phys. 2012, 136, 204108. Olsen, J.; Roos, B. O.; Jorgensen, P.; Jensen, H. J. A. J.



Page 12 of 14 12

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

43

44

45

46

47 48 49

50

51

52

53

54

55

56

57

58

59

60

61

62

63 64 65 66

67

68

69

70

71

72

Chem. Phys. 1988, 89, 2185. Malmqvist, P. A.; Rendell, A.; Roos, B. O. J. Phys. Chem. 1990, 94, 5477–5482. Ma, D.; Li Manni, G.; Gagliardi, L. J. Chem. Phys. 2011, 135, 044128. ´ Jeszenszki, P.; Zoboki, T. J. Surj´ an, P. R.; Szabados, A.; Math. Chem. 2012, 50, 534–551. Johnson, P. A.; Limacher, P. A.; Baerdemacker, S. D.; Neck, D. V.; Bultinck, P. Comp. Theor. Chem. 2013, 1003, 101–113. Neuscamman, E. Phys. Rev. Lett. 2012, 109, 203001. Mazziotti, D. A. Chem. Rev. 2012, 112, 244–262. Buenker, R. J.; Peyerimhoff, S. D. Theor. Chim. Acta. 1974, 35, 33–58. Buenker, R. J.; Peyerimhoff, S. D. Theor. Chim. Acta. 1975, 39, 217–228. Mahapatra, U. S.; Datta, B.; Bandyopadhyay, B.; Mukherjee, D. Adv. Quantum Chem. 1998, 30, 163–193. Andersson, K.; Malmqvist, P.-A.; Roos, B. O. J. Chem. Phys. 1992, 96, 1218. Yanai, T.; Kurashige, Y.; Neuscamman, E.; Chan, G. K.L. J. Chem. Phys. 2010, 132, 24105. Kurashige, Y.; Yanai, T. J. Chem. Phys. 2011, 135, 94104. Guo, S.; Watson, M. A.; Hu, W.; Sun, Q.; Chan, G. K.-L. J. Chem. Theory Comput. 2016, 12, 1583–91. Wouters, S.; Limacher, P. A.; Neck, D. V.; Ayers, P. W. J. Chem. Phys. 2012, 136, 134110. Roemelt, M.; Guo, S.; Chan, G. K.-L. J. Chem. Phys. 2016, 144, 204113. Yanai, T.; Chan, G. K.-L. J. Chem. Phys. 2006, 124, 194106. Neuscamman, E.; Yanai, T.; Chan, G. K.-L. Int. Rev. Phys. Chem. 2010, 29, 231. Zgid, D.; Ghosh, D.; Neuscamman, E.; Chan, G. K.-L. J. Chem. Phys. 2009, 130, 194107. Saitow, M.; Kurashige, Y.; Yanai, T. J. Chem. Phys. 2013, 139, 44118. Saitow, M.; Kurashige, Y.; Yanai, T. J. Chem. Theory Comput. 2015, 11, 5120. Hanauer, M.; K¨ ohn, A. Chem. Phys. 2012, 401, 50–61. Dyall, K. G. J. Chem. Phys. 1995, 102, 4909–4918. Fink, R. F. Chem. Phys. Lett. 2006, 428, 461–466. Roos, B. O.; Andersson, K. Chem. Phys. Lett. 1995, 245, 215–223. Roos, B. O.; Andersson, K.; F¨ ulscher, M. P.; SerranoAndrés, L.; Pierloot, K.; Merch´ an, M.; Molina, V. J. Mol. Struct:THEOCHEM 1996, 388, 257–276. Neuscamman, E.; Yanai, T.; Chan, G. K.-L. J. Chem. Phys. 2009, 130, 124102. Engels-Putzka, A.; Hanrath, M. J. Chem. Phys. 2011, 134, 124106. Epifanovsky, E.; Wormit, M.; Ku´s, T.; Landau, A.; Zuev, D.; Khistyaev, K.; Manohar, P.; Kaliman, I.; Dreuw, A.; Krylov, A. I. J. Comp. Chem. 2013, 34, 2293– 2309. Solomonik, E.; Matthews, D.; Hammond, J.; Demmel, J. Cyclops Tensor Framework: Reducing Communication and Eliminating Load Imbalance in Massively Parallel Contractions. Parallel Distributed Processing (IPDPS), 2013 IEEE 27th International Symposium on. 2013; pp 813–824. Calvin, J. A.; Lewis, C. A.; Valeev, E. F. Scalable taskbased algorithm for multiplication of block-rank-sparse

73 74

75

76

77 78

79

80

81 82

83

84

85

86

87 88

89

90

91

92

93

94

95

96

97 98

99 100

101 102

103

104 105

matrices. Proceedings of the 5th Workshop on Irregular Applications Architectures and Algorithms - IA3 ’15. New York, New York, USA, 2015; pp 1–8. Hirata*, S. J. Phys. Chem. A 2003, 107, 9887–9897. Cleland, D.; Booth, G. H.; Overy, C.; Alavi, A. J. Chem. Theory Comput. 2012, 8, 4138–4152. Knizia, G.; Adler, T. B.; Werner, H.-J. J. Chem. Phys. 2009, 130, 54104. Balabanov, N. B.; Peterson, K. A. J. Chem. Phys. 2005, 123, 64107. Reiher, M.; Wolf, A. J. Chem. Phys. 2004, 121, 2037–47. Werner, H.-J.; Knowles, P. J. J. Chem. Phys. 1985, 82, 5053. Knowles, P. J.; Werner, H.-J. Chem. Phys. Lett. 1985, 115, 259–267. Woon, D. E.; Dunning, T. H. J. Chem. Phys. 1993, 98, 1358. Dunning, T. H. J. Chem. Phys. 1989, 90 . Jacobsen, E. N.; Zhang, W.; Muci, A. R.; Ecker, J. R.; Deng, L. J. Am. Chem. Soc. 1991, 113, 7063–7064. Zhang, W.; Loebach, J. L.; Wilson, S. R.; Jacobsen, E. N. J. Am. Chem. Soc. 1990, 112, 2801–2803. Linker, T. Angew. Chem. Int. Ed. Engl. 1997, 36, 2060– 2062. Fu, H.; Look, G. C.; Zhang, W.; Jacobsen, E. N.; Wong, C. H. J. Org. Chem. 1991, 56, 6497–6500. Teixeira, F.; Mosquera, R. A.; Melo, A.; Freire, C.; Cordeiro, M. N. D. S. Phys. Chem. Chem. Phys. 2014, 16, 25364–25376. Stein, C. J.; Reiher, M. J. Chem. Theory Comput. 2016, Bogaerts, T.; Wouters, S.; VanDerVoort, P.; VanSpeybroeck, V. Chem. Cat. Chem. 2015, 7, 2711–2719. Ma, D.; Li Manni, G.; Gagliardi, L. J. Chem. Phys. 2011, 135, 044128. Sears, J. S.; Sherrill, C. D. J. Chem. Phys. 2006, 124, 144314. Ivanic, J.; Collins, J. R.; Burt, S. K. J. Phys. Chem. A 2004, 108, 2314–2323. Wouters, S.; Bogaerts, T.; Van Der Voort, P.; Van Speybroeck, V.; Van Neck, D. J. Chem. Phys. 2014, 140, 241103. Linde, C.; ˚ Akermark, B.; Norrby, P.-O.; Svensson, M. J. Am. Chem. Soc. 1999, 121, 5083–5084. Abashkin, Y. G.; Collins, J. R.; Burt, S. K. Inorg. Chem. 2001, 40, 4040–4048. Abashkin, Y. G.; Burt, S. K. J. Phys. Chem. B 2004, 108, 2708–2711. Abashkin, Y. G.; Collins, J. R.; Burt, S. K. Inorg. Chem. 2001, 40, 4040–4048. Zgid, D.; Nooijen, M. J. Chem. Phys. 2008, 128, 144116. Ghosh, D.; Hachmann, J.; Yanai, T.; Chan, G. K. L. J. Chem. Phys. 2008, 128, 144117. Sharma, S. J. Chem. Phys. 2015, 142, 024107. Olivares-Amaya, R.; Hu, W.; Nakatani, N.; Sharma, S.; Yang, J.; Chan, G. K.-L. J. Chem. Phys. 2015, 142, 034102. Sun, Q. J. Comp. Chem. 2015, 36, 1664–1671. Kurashige, Y.; Chan, G. K.-L.; Yanai, T. Nat Chem 2013, 5, 660–666. Sharma, S.; Sivalingam, K.; Neese, F.; Chan, G. K.-L. Nat Chem 2014, 6, 927–933. Yanai, T.; Shiozaki, T. J. Chem. Phys. 2012, 136, 84107. Sharma, S.; Yanai, T.; Booth, G. H.; Umrigar, C. J.; Chan, G. K.-L. J. Chem. Phys. 2014, 140, 104112.


Page 13 of 14


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

106

107

108 109

110

111

112

Torheyden, M.; Valeev, E. F. J. Chem. Phys. 2009, 131, 171103. Kersten, J. A. F.; Booth, G. H.; Alavi, A. J. Chem. Phys. 2016, 145, 054117. Ten-no, S. Chem. Phys. Lett. 2007, 447, 175–179. Shiozaki, T.; Werner, H.-J. J. Chem. Phys. 2010, 133, 141103. Shiozaki, T.; Knizia, G.; Werner, H.-J. J. Chem. Phys. 2011, 134, 34113. Fink, R.; Staemmler, V. Theor. Chim. Acta. 1993, 87, 129–145. Riplinger, C.; Neese, F. J. Chem. Phys. 2013, 138, 034106.

113

114

Guo, Y.; Sivalingam, K.; Valeev, E. F.; Neese, F. J. Chem. Phys. 2016, 144, 094111. For a contraction with n input tensors (n > 2), our heuristic algorithm first choses the two input tensors which, if contracted over their common indices, produce the smallest intermediate tensor. This binary contraction is then evaluated, and the two input tensors are replaced by the single output tensor, resulting in a n − 1 tensor contraction. This is repeated until a binary contraction is reached, which is evaluated with transpositions and matrix multiplications.



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

Page 14 of 14

Internal contraction + Matrix product states

1.96

d-d excitation

0.16 Singlet state

1.00

1.04 ACS Paragon Plus Environment

Triplet state

Combining Internally Contracted States and Matrix Product States To

Recommend Documents