Embedded Multireference Coupled Cluster Theory - Journal of

Subscriber access provided by READING UNIV

Article

Embedded Multireference Coupled Cluster Theory David James Coughtrie, Robin Giereth, Daniel Kats, Hans-Joachim Werner, and Andreas Köhn J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01144 • Publication Date (Web): 18 Jan 2018 Downloaded from http://pubs.acs.org on January 18, 2018

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 61 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

Embedded Multireference Coupled Cluster Theory David J. Coughtrie,∗ Robin Giereth, Daniel Kats, Hans–Joachim Werner, and Andreas Köhn∗ Institute for Theoretical Chemistry, University of Stuttgart, 70569 Stuttgart, Germany E-mail: [email protected]; [email protected] Phone: +49 (0)711 / 685-64697. Fax: +49 (0)711 / 685-64442

Abstract Internally contracted multireference coupled cluster (icMRCC) theory is embedded within multireference perturbation theory (MRPT) to calculate energy differences in large strongly correlated systems. The embedding scheme is based on partitioning the orbital spaces of a complete active space self consistent field (CASSCF) wave function, with a truncated virtual space constructed by transforming selected projected atomic orbitals (PAOs). MRPT is applied to the environment using a subtractive embedding approach, that also allows for multilayer embedding. Benchmark calculations are presented for biradical bond dissociation, spin splitting in a heterocyclic carbene and hydrated Fe(II), and for the super-exchange coupling constant in solid nickel oxide. The method is further applied to two large transition metal complexes with a triple-ζ basis set - an iron complex with 175 atoms and 2939 basis functions, and a nickel complex with 231 atoms and 4175 basis functions.

1

ACS Paragon Plus Environment

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

1

Introduction

Large multireference systems are widespread and of enduring interest, due to their intriguing chemical behaviour. Examples include transition metal complexes, 1–4 single molecule magnets, 5–7 bond breaking, 8–10 poly-radical species, 8,11 conical intersections in photochemistry, 12–14 and bioluminescence. 15,16 Performing calculations on these systems can provide insight into their properties and mechanisms, complementing experimental studies and helping to rationalize their behaviour. Recent years have seen the advent of algorithms for treating large-scale single reference systems with accurate wave-function methods, in particular coupled cluster (CC) theory. 17–23 These implementations achieve linear scaling with respect to the number of electrons in the system, by using localised occupied orbitals and adapted subsets of the virtual orbital space, such as pair natural orbitals (PNOs). 19–30 Multireference systems cannot be treated with these methods, however, because they exhibit several low-lying energy levels that cannot be qualitatively described with a single configuration. Development of linear scaling multireference techniques is progressing, for example linear scaling multireference perturbation theory (MRPT) has been applied to large transition metal complexes. 31–33 Furthermore a recent local PNO implementation of Mukherjee’s multireference coupled cluster (LPNO-Mk-MRCC) method, 34,35 portends a future where even MRCC 36 theory can be used for large-scale systems. At present however, LPNO-MkMRCC inherits the limitations of the underlying Mk-MRCC method, 37,38 principally that only small active spaces can be used, and that the energy is not invariant to orbital rotations within the active space. An established alternative for treating large systems is embedding. 39 Frequently a system contains a target region of in principle limited size, but which is surrounded by a much larger ‘inactive’ environment. This environment cannot be simply neglected, however, as it interacts with the target region and constrains its structure. Instead, the environment may be approximately described with less accurate but more affordable means, allowing 2


Page 2 of 61

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


computationally expensive theories to be focused on the target region of interest. Many embedding schemes are available, such as the widespread quantum mechanics/molecular mechanics (QM/MM) approach. 40–42 The environment region is calculated classically using molecular mechanics force fields, and the central target region treated quantum mechanically, often using density functional theory (DFT) or MRPT, 12,43,44 though higher level correlated methods can also be used. A closely related method is the subtractive ONIOM scheme. 45–47 Alternatively, exact embedding can be achieved through partitioning of the total electron density into environment and target regions. The Kohn–Sham equations for one region are solved in the presence of an embedding potential arising due to the charge density of the other, leading to DFT-in-DFT, wave-function theory (WFT)-in-DFT and WFT-in-WFT schemes. 48–52 The projector embedding technique 53–55 is an elegant example of this approach, that also allows for truncation of the virtual space, 56,57 and was recently applied to a large enzyme within a wider QM/MM environment. 58 Density matrix embedding theory (DMET) is a related WFT-in-WFT embedding scheme, 59–64 in which a set of bath orbitals describing the entanglement of the target region with the environment is constructed via a Schmidt decomposition. Partitioning of orbital subspaces has also been used successfully for WFT-in-WFT embedding. 50,65–68 Typically this involves performing an initial self-consistent field (SCF) calculation, localising the orbitals, and then applying a higher level correlated method to those orbitals in the target region. This is the approach followed by the regional local coupled cluster method of Mata et al., 69 the multi-level extension of the cluster-in-molecule method, 70,71 the multi-level coupled cluster theory of Myhre and Koch, 72–76 and is also a component of the method of increments. 77–83 In this work we draw on several aspects of the methods described above in order to embed internally contracted MRCC theory (icMRCC) 84–89 within a complete active space SCF (CASSCF) 90 wave function, and subsequently within MRPT. icMRCC has the advanta-

3



geous properties of size extensivity, energy invariance with respect to active orbital rotations, and can treat moderately sized active spaces at feasible cost. We peruse a WFT-inWFT approach based on the orbital space partitioning of a CASSCF calculation, as this allows for embedding within a multireference wave function. Projected atomic orbitals 91 (PAOs) are used to reduce the size of the virtual space without introducing significant errors in calculated energy differences. This is combined with a simple subtractive embedding procedure, in the fashion of the ONIOM model, that allows MRPT to be applied to the environment and for multi-layer embedding. We refer to the overall approach as an embedding scheme, rather than a multi-level method (such as the regional local coupled cluster approach). Treating the environment at only the CASSCF level can be considered a multi-level calculation, as it involves just a single calculation that treats individual interactions at different levels depending on the location of their participating orbitals. In contrast, applying MRPT to the environment requires combining several distinct calculations using different levels of theory, typical of embedding. Benchmark calculations on energy differences are presented to assess the embedding error, for the diradical homolytic bond cleavage of two small alkanes, spin splitting in a small carbene molecule, calculation of the super-exchange coupling constant in solid nickel oxide, and the spin splitting in a hydrated Fe(II) complex. Finally we present embedded icMRCC calculations for two large transition metal complexes, consisting of several hundred atoms and using several thousand basis functions, demonstrating the applicability of this approach to large-scale multireference systems of real world interest. This paper is arranged as follows: section 2 gives an overview of the relevant theory as well as details of the embedding scheme, section 3 presents benchmark calculations, and section 4 presents results for the two large transition metal complexes.

4


Page 4 of 61

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


2 2.1

Theory icMRCC

MRCC theory applies the exponential cluster operator to a CASSCF reference wave function. In the CASSCF wave function the orbitals are classified into three types: inactive orbitals that are always doubly occupied, active orbitals whose occupation can vary between two and zero, and virtual orbitals that are unoccupied. The wave function is a linear combination of all possible determinants |Φi i that can be made from distributing the active electrons (those left over once the inactive space is filled) amongst the active orbitals,

|Ψ0 i =

n X

|Φi ici ,

i=1

where ci are the coefficients of each determinant. The goal is to include dynamic correlation by adding electronically excited determinants to this wave function. There are many versions of MRCC theory that have been developed over the last four decades. 36 Here we use the internally contracted variant (icMRCC), 84–89 because it retains the desirable properties of size consistency and extensivity, is invariant to orbital rotations within the active space, and is not prohibitively expensive. The wave function ansatz is ˆ

ˆ

|Ψi = eT |Ψ0 i = eT

n X

|Φi ici .

i=1

The cluster operator is restricted to only allow single and double excitations (icMRCCSD), with perturbative approximations used for the triple excitations (icMRCCSD(T)). 92 The method has been applied to the calculation of a potential surface, 93 transition metal compounds, 94 and excitation energies. 95 The cluster amplitudes within Tˆ and the coefficients ci are optimised until they are mutually consistent; however, the molecular orbitals remain fixed. The working equations needed to perform this optimisation are significantly more complicated than for single 5



reference CC theory, because certain classes of excitation operator no longer commute. A consequence of this extra complexity is that the theory is too cumbersome to be implemented by hand, and we must therefore rely on automated techniques to simplify and evaluate the working equations. We use the general contraction code (GeCCo), a program developed within our group, that can manipulate expressions algebraically as well as evaluate them numerically. Our present implementation of icMRCC is intended to assess and demonstrate the formal properties of the theory, but has not yet been optimised for fast calculations. The scaling behaviour of icMRCC is particularly relevant to the present work. 86,88 To aid discussion we denote the sizes of the inactive, active and virtual orbital spaces as Minact , Mact and Mvirt respectively. For a small fixed size active space where Mvirt Mact (for 4 2 ), and Mvirt instance when Mact ≤ 5, as is the case here), icMRCCSD scales as O(Minact 3 4 icMRCCSD(T) scales as O(Minact Mvirt ), analogous to their single reference counterparts.

Our goal is to therefore reduce the size of the virtual and inactive spaces used in the correlated calculation.

2.2

Multireference Perturbation Methods

An immediate advantage of embedding is that much more affordable methods can be used for the environment of the system, such as MRPT. In MRPT a number of choices for the ˆ 0 exist, leading to a variety of methods. Complete active space zeroth order Hamiltonian H ˆ 0, perturbation theory to 2nd order 96,97 (CASPT2) uses a generalized Fock operator for H that consists of effective one-electron operators only. 98 This sometimes suffers from the intruder state problem, however, in which virtual orbitals have energies similar to those of the active space, leading to small denominators and unphysically large second order ˆ 0 can ameliorate this problem, 98–101 however some arbienergies. Using a level shift in H trariness is introduced into the final energy. The n-electron valence states perturbation ˆ 0 , which takes theory to 2nd order 102–104 (NEVPT2), uses the Dyall Hamiltonian for H 6


Page 6 of 61

Page 7 of 61 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


into account two-electron interactions within the active space. A significant advantage of NEVPT2 is that it is free of intruder states. 102 Importantly for our purposes, linear scaling local variants of these methods have been recently developed based on PNOs, allowing them to be applied to systems of hundreds of atoms. Specifically the domain-based local DLPNO-NEVPT2 31 method available in the ORCA program, 105 and the PNO-CASPT2 32,106 method available in the Molpro program. 107,108 We use the latter for full calculations on the large complexes in section 4.

2.3

Embedding within a CASSCF wave function

Localised molecular orbitals can be sensibly partitioned into a target set necessary for describing a reaction, and an environment set that does not contribute significantly to energy differences. The former are typically spatially localised around the active orbitals, though may not be the highest in energy. Localisation of the occupied orbital spaces (without mixing core, inactive, and active orbitals) is performed using the intrinsic bonding orbital (IBO) method of Knizia. 109 The IBO localisation scheme is a more stable variant of the Pipek-Mezey scheme, producing compact chemically intuitive orbitals that only weakly depend on the orbital basis set. The orbitals to be embedded are chosen directly using chemical intuition. Automated ‘black box’ selection procedures are available in which a set of embedded atoms is chosen and a partial charge criteria used to determine the associated orbitals. 50,58,69 We do not use them here, however, because the high computational cost of icMRCC calculations using our present non-optimal implementation requires a more limited selection. The entire active space is always retained within the embedding region, along with those inactive orbitals that have been chosen and a subset of virtual functions (whose construction is discussed below). No further relaxation of these embedded orbitals is performed. The remaining (inactive) environment orbitals, including the low-energy atomic core, are frozen and left uncorrelated. 110,111 The Hamiltonian for the embedded electrons thus be7



Page 8 of 61

comes (in spin-free second quantized form, using atomic units)

ˆ = E0 + H

emb X pq

emb

X ˜ p Eˆ q + 1 g pq Eˆ rs h q p 2 pqrs rs pq

(1)

pq = (pr|qs), where the indices p, q, r and s label embedded molecular orbitals of any kind, grs rs and Eˆpq and Eˆpq are spin-summed single and double excitation operators respectively. E0 is

the energy of the frozen electrons,

E0 = hnuc + 2

froz X

hii

+

froz X

ji ij 2gji − gji

ij

i

where i and j label the frozen orbitals, hii is the scalar part of the one-body Hamiltonian, and hnuc is the interaction energy between the nuclei. The second term in equation 1 describes the mean-field interaction of the embedded electrons with the frozen environment, where ˜ p = hp + h q q

froz X

ip pi 2giq − giq .

i

This procedure is exactly the same as that of the ‘frozen-core approximation’, but is applied to the environment valence orbitals as well. The virtual space is truncated by transforming selected atom centred PAOs. A set of PAOs spanning the entire system is constructed using the projector 91

Pˆ = 1 −

occ X

|nihn|,

(2)

n

where the sum runs over the active, inactive and core orbital spaces, including any that have been frozen. Applying this projector to the atomic orbital basis generates a new set of basis functions that are orthogonal to the occupied spaces, but not to one another. Each PAO is centred on the same atom as its parent atomic orbital, allowing the PAOs to be easily partitioned based on the location of the atoms. We therefore classify the atoms

8


Page 9 of 61 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


as either host atoms or environment atoms depending on whether they are in the embedded region of the molecule or its environment. When the total IBO partial charge exerted by the embedded orbitals on a given atom exceeds 0.2 atomic units, the atom is classified as a host atom, otherwise it is an environment atom. PAOs located on environment atoms are discarded, and those on host atoms are transformed by diagonalising their overlap matrix, S = U> λU.

(3)

The diagonal elements of λ are the eigenvalues of S, and the columns of U its eigenfunctions. Functions with eigenvalues smaller than a threshold value of 1.0 × 10−8 are considered redundant and removed. The resulting matrix is used to transform the remaining PAOs into a smaller set of orthogonal virtual orbitals - called the virtual domain - located ‘above’ the embedded IBOs. Local PAO methods use the same approach to construct virtual domains for pairs of correlated orbitals; 91,112–117 in the present work, however, only a single domain is constructed for the total embedding region. The virtual domain may be extended to achieve higher accuracy, as is done in existing local correlation methods; 17–23,32,118–120 however we do not do so here, due to the extra cost of a larger virtual space. Nevertheless, an IBO in the centre of the embedding region will experience an effective domain extension due to the surrounding host atoms. Similarly an IBO on the fringe of the embedding region experiences a domain extension towards the embedded region, and we therefore expect the error associated with truncating the virtual space (the virtual domain error) to be small. This error is measured numerically in section 3, and is indeed found to be small. Finally, the effective Fock matrix is block diagonalised with respect to the embedded inactive, active and virtual orbital spaces. The resulting pseudo-canonical molecular orbitals speed up convergence of the correlated methods that are subsequently executed by the GeCCo program. The regional CC method of reference 69 uses a very similar orbital selection procedure to 9



Page 10 of 61

that described here; however, in that case the embedded region is treated using fully local single reference CC within Molpro.

2.4

Subtractive Multi-layer Embedding

The embedding procedure described in section 2.3 only allows for embedding within the surrounding CASSCF reference molecular orbitals. Any influence of dynamic correlation within the environment, and more crucially, between the environment and the embedded region is absent. These effects can be included using a simple subtractive embedding approach, similar to that of the ONIOM scheme 46 and that proposed in reference 68. We define the subtractive embedding energy as

emb emb E sub (A, B) = Elow (A + B) − Elow (A) + Ehigh (A)

(4)

where E emb is the energy calculated using the embedding procedure described in section 2.3. Arguments A and B denote the embedded and environment orbital subspaces respectively, and subscripts denote the accuracy of the methods applied to each region. All energies on the right hand side of equation 4 are calculated using the same CASSCF reference molecular orbitals, and therefore the subtractive embedding is acting on the dynamic and static correlation energies only. Several nested layers of correlated methods may be easily combined through repeated application of equation 4. In principle the size of the embedded region should be increased until the energy difference under study is converged. In practice the size of the embedding region is limited by the high cost of icMRCC calculations, though in future we hope to see the development of much more efficient variants of MRCC theory that are able to overcome this limitation. Nevertheless, at present we can study the effect of higher-order correlation even if convergence at the MRCC level is not complete, by monitoring the convergence of MRPT embedded within CASSCF.

10


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


2.5

Density Fitted Integral Transformation

icMRCC is implemented in GeCCo using a molecular orbital based algorithm, and hence requires a full integral transformation. 121–123 Although the number of embedded orbitals is small relative to the full set, this can still be prohibitively expensive for large systems with thousands of basis functions. Straightforward application of density fitting procedures circumvents this problem, and is used with the iron and nickel complexes in section 4. The cost of the transformation is dominated by the two-electron integrals, in particular those involving the truncated virtual space; however, density fitting allows them to be decomposed into a contraction over three-index objects,

(pq|rs) =

X

¯ ¯ (B|pq)( B|rs),

¯ B

with the three-index objects constructed via,

¯ (B|rs) =

X

(G−1 )BA ¯

A

X µ

Cµr

X

Cνs (A|µν).

ν

The labels µ and ν refer to atomic basis functions, A and B to auxiliary atom centred fitting functions, and Cµr are the embedded molecular orbital coefficients. The two electron three centre integrals are defined as, Z (A|µν) =

−1 dr 1 dr 2 χµ (r 1 )χν (r 1 )r12 χA (r 2 )

where r12 = |r 1 − r 2 |. The matrix G comes from the Cholesky decomposition of the Coulomb interaction over the auxiliary functions. Existing integral blocking and screening machinery within Molpro was used to implement the fitting procedure, adapted from the previous implementations in references 17,124 and 125.

11



3

Numerical testing of the Embedding Scheme

Here we assess the accuracy of the embedding procedure for several test cases. All calculations were performed on a single processor, using a development version of Molpro 107,108 and the general contraction code, GeCCo. IBO localization was performed using symmetric orbital orthogonalisation, and a localization exponent of 4. Molecular point group symmetry was exploited in the canonical calculations. All NEVPT2 and CASPT2 calculations were performed with fully internally contracted canonical algorithms, implemented in GeCCo, except those requiring a level shift, which were performed with Molpro. The PNO-CASPT2 calculations on the large complexes in section 4 were also performed using Molpro. Numerical values used to make the plots are given in the supporting information, along with all geometries, and lists of the host atoms for each system.

3.1

Alkane homolytic bond breaking

Breaking of the central carbon-carbon bond in butane and hexane cannot be properly described with low-level single reference methods, which predict dissociation into charged fragments and not the radicals of the true ground state. 8,9 We use these systems as benchmark cases for embedded MRCC. Reference molecular orbitals were calculated using the state-specific CASSCF method, 126–129 with an active space of two electrons in the central carbon-carbon σ bonding and antibonding orbitals. To ensure that these two radical orbitals remained in the active space we did the following: initially the CASSCF calculation at each geometry was performed using C2h point group symmetry, in which the highest energy orbitals of the Ag and Bu irreducible representations correspond to the radical orbitals. This wave function was then used as a starting guess for a second non symmetry-adapted calculation. Furthermore, the bond length was decreased from dissociation towards equilibrium, and the wave function of the previous geometry used as a starting guess for the next.

12


Page 12 of 61

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


Region A

Region B

Region C

Region D

Figure 1: Charge densities of the embedding regions in hexane homolytic bond breaking. Equivalent regions are also used with butane. Density plots rendered using the program Avogadro. 134,135 The two halves of the molecule were moved in a direction parallel to the central C-C bond, with the rest of the structure frozen at the equilibrium geometry. 130 The correlation consistent polarized valence triple-ζ basis (cc-pVTZ) of Dunning and co-workers 131–133 was used throughout. Four embedding regions were considered, shown in figure 1. Region A contains the active space only, region B the active space and the nearest two C-C bonds, region C the active space and the four nearest C-H bonds, and finally region D is the union of regions B and C. Table 2 details the number of orbitals in each space for each embedding region. Bond breaking potential surfaces for the alkanes were calculated relative to their energies at a fragment distance of 2000 Å, using icMRCCSD embedded in CASSCF and NEVPT2. Figure 2 shows these surfaces with their asymptotic limits aligned to zero. The error in each surface relative to icMRCCSD on the full molecule is also presented, with the largest unsigned errors given in table 1 for clarity. Both alkanes exhibit similar behaviour due to the localised nature of their radicals, which is particularly well suited to an embedding approach. For hexane, embedding with region D exhibits consistent errors across the range of distances, of at most 2.7 kJ mol−1 for icMRCCSD in CASSCF, −4.2 kJ mol−1 for icMRCCSD in NEVPT2 and −1.1 kJ mol−1 for icMRCCSD in CASPT2. For these systems embedding within only CASSCF performs well, though this is likely coincidental. 13



Butane

Hexane

−100 −200 −300

CASSCF NEVPT2 CASPT2 icMRCCSD

−400

CASSCF NEVPT2 CASPT2 icMRCCSD

−500 30.0 20.0 10.0 0.0 −10.0 −20.0 −30.0 30.0 20.0 10.0 0.0

CASSCF Region A Region B Region C Region D icMRCCSD

CASSCF Region A Region B Region C Region D icMRCCSD

NEVPT2 Region A Region B Region C Region D icMRCCSD

NEVPT2 Region A Region B Region C Region D icMRCCSD

CASPT2 Region A Region B Region C Region D icMRCCSD

CASPT2 Region A Region B Region C Region D icMRCCSD

−10.0 −20.0 −30.0 30.0 20.0 10.0

2.0

1.5

1.0

0.5

0.0

−30.0

2.0 0.0

−20.0

1.5

−10.0

1.0

0.0

0.5

Error in Ed / kJ mol−1


Ed / kJ mol−1

0


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

Page 14 of 61

(rC−C − req) / Å

Figure 2: Homolytic C-C bond breaking in butane (left column) and hexane (right column). The first row shows the dissociation surface. The remaining rows show the error when embedding icMRCCSD in CASSCF, NEVPT2 and CASPT2 respectively. All errors are relative to icMRCCSD/cc-pVTZ on the full molecule. The largest unsigned errors in the lower six panels can be found in table 1.

14


Page 15 of 61 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 1: Largest unsigned errors in dissociation surfaces of the alkanes butane and hexane (see figure 2) in kJ mol−1 , when embedding icMRCCSD in CASSCF, NEVPT2, and CASPT2 (defined as the maximum absolute error in Ed (r − req ) relative to icMRCCSD/ccpVTZ on the full molecule, over the interval r − req = [0.0, 2.0]). Environment

icMRCCSD Region None

A

B

C

D

Butane/cc-pVTZ CASSCF NEVPT2 CASPT2

67.3 10.9 4.9

28.5 22.0 5.0

28.1 16.7 5.2

15.1 10.1 4.3

1.5 3.6 1.3

Hexane/cc-pVTZ CASSCF NEVPT2 CASPT2

67.1 11.1 5.5

28.2 22.7 4.9

27.8 17.3 5.0

15.1 10.7 4.2

2.7 4.2 1.1

icMRCCSD embedded within NEVPT2 using region D leads to more consistent errors across the potential curve than embedding within CASSCF or CASPT2, and the error vanishes more quickly at stretched geometries. Region D consists of only 176 virtual orbitals, a reduction of 27% and 50% for butane and hexane respectively, and only 6 inactive orbitals (see table 2 for further details). The smaller virtual space in particular trans4 scaling behaviour of icMRlates into significant computational savings due to the Mvirt

CCSD. Other embedding regions exhibit larger errors that vary as the distance changes, suggesting that region D is the minimum set of intrinsically important orbitals that must be treated at the multireference level. Benchmark calculations are prohibitively expensive for larger alkanes; nevertheless, the consistency of the embedding error for the molecules presented here leads us to expect similar results for larger systems. In such cases, embedding region D (or a small extension of it) would allow for significant savings in computational effort. Finally, we comment of the contrasting behaviour of the non-embedded CASPT2 and NEVPT2 results. CASPT2 is known to underestimate bond energies by around 3-5 kJ mol−1 , 136 consistent with our calculations. Our non-embedded NEVPT2 result, however, overesti-

15



Page 16 of 61

Table 2: Number of molecular orbitals for each embedding region of the alkanes and NHCMe. The percentage of the full doubly-occupied and virtual spaces are given in brackets. The inactive, active and virtual orbitals are treated at the higher level. Region

Environment

Inactive (%)

Active

Butane/cc-pVTZ A B C D Full

16 14 12 10 4

0 (0) 2 (13) 4 (25) 6 (38) 12 (75)

2 2 2 2 2

60 (25) 120 (50) 116 (45) 176 (73) 242 (100)

Hexane/cc-pVTZ A B C D Full

24 22 20 18 6

0 (0) 2 (8) 4 (17) 6 (25) 18 (75)

2 2 2 2 2

60 (17) 120 (34) 116 (33) 176 (50) 350 (100)

NHC-Me/cc-pVTZ A B C D Full

25 21 19 17 7

0 (0) 4 (16) 6 (24) 8 (32) 18 (72)

2 2 2 2 2

90 (31) 90 (31) 150 (51) 150 (51) 295 (100)

16


Virtual (%)

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


mates the bond strength of hexane by 11.1 kJ mol−1 . We attribute the poorer NEVPT2 result to the small size of active space used here. At distances in the regime of chemical bonding, charge transfer configurations in the reference wave function become important. The ionic nature of these configurations distorts the wider σ bonding of the molecule, a process called dynamic σ polarisation. 137–140 A larger active space is needed to properly describe this effect. We have performed further calculations using the orbitals of region B as the active space, that result in CASPT2 and NEVPT2 bond energy errors of 5.9 and −1.0 kJ mol−1 respectively (measured from equilibrium, relative to icMRCCSD calculations with this same larger active space). The NEVPT2 result is therefore significantly improved when using a zeroth order wave function that can capture the dynamic σ polarisation. Schapiro and co-workers 141 suggest that CASPT2 does not suffer from a similar error because it treats all three orbital spaces equivalently with a single particle Hamiltonian. NEVPT2 on the other hand, treats the active space distinctly by including bielectronic effects in the Dyall Hamiltonian, making it more sensitive to changes in the zeroth order wave function.

3.2

A Small N-Heterocyclic Carbene

We now consider the vertical spin splitting of a persistent N-heterocyclic carbene, 142–145 1,3-di(methyl)imidazol-2-ylidene (referred to as NHC-Me from here on) shown in figure 3. This molecule has a five membered ring with the carbene centre between two nitrogen atoms, each of which hosts a methyl group. The highest occupied molecular orbital (HOMO) is the sp2 hybridised lone pair on the middle carbon in the ring, and the lowest unoccupied molecular orbital (LUMO) is the orthogonal p orbital of the same atom. The ground state is a singlet with the filled sp2 orbital stabilised by the neighbouring nitrogen atoms, that withdraw electron density through the σ bonds, and donate electron density to the empty p orbital via the π system. Reference orbitals were generated using state-specific CASSCF theory with an active space 17



Region A

Region B

Region C

Region D

Error in ∆ETS / kJ mol−1

Figure 3: Charge density plots of the correlated embedding regions of the carbene NHCMe.

Virtual domain error / kJ mol−1

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

Page 18 of 61

10.0 5.0 0.0 −5.0 −10.0 −15.0 −20.0 −25.0 −30.0 −35.0 −40.0 −45.0 −50.0 0.8 0.6 0.4 0.2 0.0 −0.2 −0.4 −0.6 −0.8 None

CASSCF → NEVPT2 CASSCF → CASPT2 CASSCF → icMRCCSD NEVPT2 → icMRCCSD CASPT2 → icMRCCSD

A

B C Correlated region

D

Full

Figure 4: NHC-Me triplet-singlet splitting energy (∆ETS = ET − ES ) with increasing embedding region, relative to icMRCCSD/cc-pVTZ on the whole molecule (for which ∆ETS = 432.1 kJ mol−1 ). The lower panel shows the virtual domain error for each region. consisting of two electrons in the HOMO and LUMO, using the cc-pVTZ basis set. Embedded correlated calculations were then performed using four nested embedding regions, shown in figure 3. Region A contains only the active space, region B further includes the two adjoining σ and π orbitals surrounding the carbene centre, region C incorporates the next two nearest C-N σ bonds in the ring, and finally region D contains the entire ring structure. Details of the number of orbitals in each space are given in table 2. The error in the vertical triplet-singlet splitting energy relative to icMRCCSD on the whole molecule is shown in figure 4. When embedding in CASSCF, convergence of the splitting 18


Page 19 of 61 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


energy with respect to increasing the embedding region is not smooth. Embedding icMRCCSD in NEVPT2 or CASPT2, however, results in much smoother convergence, with errors in region D of only 1.1 and −0.5 kJ mol−1 respectively. This is because the perturbative methods provide a sufficiently accurate description of the influence of the environment on the splitting energy, unlike CASSCF, and so artefacts due to the embedding are much smaller. The splitting energy is almost entirely converged when region D is used, and we therefore expect that larger residues attached to the ring will not reduce the accuracy of the approach (though this may not be the case for residues conjugated to the ring system). Region D represents a sizeable reduction of the inactive and virtual orbital spaces, of around 50% for this example (see table 2 for more details) Larger systems with additional residues will therefore also enjoy significant reduction in computational cost. We further calculate the virtual domain error by comparing the embedded calculations to their equivalents with the full canonical virtual space. Figure 4 shows the size of the domain error for each of the embedding regions, which is generally small, and not the limiting factor on the accuracy of the embedding. The magnitude of the virtual domain error depends on the importance of the PAOs on the next nearest set of possible host atoms, referred to as boundary PAOs from here onwards. Delocalised embedded orbitals (such as in region B) would interact more strongly with the boundary PAOs, leading to a larger domain error in their absence. Compact embedded orbitals (such as those of region A) would interact less strongly with the boundary PAOs, and so their absence leads to a smaller domain error. Furthermore, the magnitude of all embedding errors is attenuated as the embedding region boundary moves further from the multireference centre (as with regions C and D) .

19



[Fe(H2 O)6 ]2+

Region A

Region B

Figure 5: Hexa-aqua Iron (II) and charge density plots of its embedding regions.

3.3

Hexa-aqua Iron [Fe(H2 O)6 ]2+

Systems containing transition metal centres are frequently multireference in nature, with several low lying configurations arising from the d orbitals of the metal centre (see table 3 of the supporting information for explicit lists of significant configurations in the reference wave functions of the systems studied here). We consider the simple test case of hydrated iron(II), 49,146 shown in figure 5. In the octahedral ligand field environment, the lowest energy singlet, triplet and quintet states have term symbols 1 A1 , 3 T1 and 5 T2 , the energetic ordering of which depends upon the Fe-O bond distance. We examine the vertical singletquintet splitting as the ligand distance is increased. The 3 T1 state is not considered here because of convergence difficulties with the embedded calculations, due to near degeneracy of two low-energy triplet states, and the unavailability of point-group symmetry adapted calculations when embedding. Reference molecular orbitals were calculated using CASSCF theory, with an active space of 6 electrons in the 5 Fe 3d orbitals, and state averaging only over the degenerate states within a given term. All calculations were performed with the cc-pVTZ basis set. The geometries of the ligand water molecules were fixed as the Fe-O distance increased, with O-H bond lengths of 0.97204 Å and H-O-H angles of 108.706 degrees, following reference 49. The equilibrium Fe-O distance (subject to the aforementioned constraints) is approximately 2.025 Å for the singlet state and 2.153 Å for the quintet, calculated by fitting a quartic polynomial to the total icMRCCSD/cc-pVTZ energies used to make figure 6 (see the supporting information for full details). Two embedding regions were considered,

20


Page 20 of 61

300.0 200.0 100.0

NEVPT2 CASPT2 icMRCCSD

0.0 50.0

Error in ∆ESQ / kJ mol−1

40.0 30.0 20.0 10.0 0.0 −10.0 NEVPT2 NEVPT2 → icMRCCSD R−A NEVPT2 → icMRCCSD R−B CASPT2 CASPT2 → icMRCCSD R−A CASPT2 → icMRCCSD R−B

−20.0 −30.0 −40.0 −50.0 1.9

Virtual domain error / kJ mol−1

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


∆ESQ / kJ mol−1

Page 21 of 61

2.0

2.1 2.2 2.3 Fe−O distance / Å

2.4

0.30 0.20 0.10 0.00 −0.10 −0.20 −0.30

Fe−O = 2.0 Å

Fe−O = 2.2 Å

A B Correlated region

A B Correlated region

Figure 6: The singlet-quintet splitting (∆ESQ = ES − EQ ) of [Fe(H2 O)6 ]2+ is shown in the upper panel, and its associated error relative to icMRCCSD/cc-pVTZ on the whole system shown in the middle panel. The lower panels show the virtual domain error, at selected Fe-O bond lengths.

21



Page 22 of 61

Table 3: Number of molecular orbitals for each region of [Fe(H2 O)6 ]2+ and the nickel oxide cluster model. The percentage of the full doubly-occupied and virtual spaces are given in brackets. The inactive, active and virtual orbitals are treated at the higher level. Region

Environment

Inactive (%)

Active

Virtual (%)

[Fe(H2 O)6 ]2+ /cc-pVTZ A 39 B 33 full 15

0 (0) 6 (15) 24 (62)

5 5 5

63 (17) 243 (65) 372 (100)

Nickel-oxide cluster model/def2-SVP A 79 B 69 C 59 Full 29

0 (0) 10 (13) 20 (25) 50 (63)

4 4 4 4

60 (45) 73 (55) 133 (100) 133 (100)

shown in figure 5. Region A contains only the active space, whilst region B also includes the six Fe-O sigma bonds. Table 3 details the specific number of molecular orbitals involved in each space. Figure 6 shows the vertical singlet-quintet splitting calculated using icMRCCSD embedded in NEVPT2 and CASPT2, along with the error relative to icMRCCSD on the entire system. Embedding region B within NEVPT2 performs well with errors of −4.1 and −2.1 kJ mol−1 in the vicinity of the singlet and triplet equilibrium bond lengths respectively. Embedding in CASPT2 fares less well however, with errors of 8.0 and 6.3 kJ mol−1 at the same bond lengths. This is because the CASPT2 method gives an unbalanced description of open shell and closed shell states, leading to an overestimate of the spin-splitting energy. 98,147 A level shift correction of 0.1 was applied to the non-embedded CASPT2 calculation, to enable convergence of the energy. The virtual domain error is also shown in the lower panels of figure 6 for two selected bond lengths, and is typically less than ±0.3 kJ mol−1 for this system. This error is larger for region B than for region A, despite the former having a greater proportion of the full virtual space. This occurs because the boundary PAOs are closer to region B (around 1.0 Å away) than they are to region A (around 2.0 Å away), and so their interaction with the embedded orbitals would be stronger, and the domain error is larger in their absence.

22


Page 23 of 61 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


3.4

Crystalline Nickel Oxide

Crystalline solids of many transition metal compounds exhibit magnetically ordered structure when unpaired electrons on nearby sites interact. This coupling can occur between neighbouring sites, called the exchange interaction, and also between next-nearest neighbours via intervening orbitals, called the super-exchange interaction. 148,149 The archetypal example of such a super-exchange interaction is seen in nickel oxide, a crystalline solid composed of Ni2+ and O2− ions in a rock-salt structure. Naked Ni2+ cations have a 3d8 3 F electronic configuration, which is split by the octahedral crystal field into the lowest energy 3 A2g state, followed by 3 T2g and 3 T1g states. 150 Super-exchange interaction between unpaired electrons on the Nickel sites further splits the lowest-energy state, into a nearly degenerate trio with singlet, triplet and quintet spin respectively. Calculations on these states require a multireference method. Previous theoretical studies 151–157 suggest that dynamic correlation arising from the bridging oxygen, as well as the response of the nearest shell of oxygen ions to charge-transfer configurations, are both important. The Heisenberg Hamiltonian is a widely used model that describes the interaction of unpaired spins between different sites. In our case we only consider the pairwise interaction between two selected nickel sites, and so the model reduces to

ˆ = −2J S ˆ1 · S ˆ2 . H ˆ is the vector of electron spin operaThe super-exchange coupling constant is J, and S tors on each of the Ni sites. The multiplet splittings of this Hamiltonian are related to J, 149 with a triplet-singlet energy gap of −2J, a quintet-triplet gap of −4J, and a quintetsinglet gap of −6J. Comparison to our calculated splitting therefore gives an ab-inito value for J. We determine J using all three splittings and report the mean average, with individual values given in the supporting information. Here we follow the approach of previous authors and employ an embedded cluster model.

23



Cluster model including pseudopotentials

Region A

Region B

Region C

Figure 7: Cluster model for solid nickel oxide with Ni2+ pseudopotentials shown as grey balls, along with charge-density plots of the three embedding regions. Two Ni2+ ions and one shell of 11 O2− ions are treated explicitly with all electrons. This cluster is embedded in a cavity of equal and opposite charge within a cubic array of 12154 point charges, that has a side length of about 46 Å. The Madelung potential is reproduced in the centre of the crystal by using the Evjen method. 158 To prevent ‘leakage’ of electron density from the explicit oxygen anions into the surrounding cloud of point charges, Ni2+ pseudo potentials are superimposed onto the first shell of nickel point charges. 159 An experimental Ni-O distance of 2.0842 Å was used throughout. 151,160 Reference molecular orbitals were calculated using the state-specific CASSCF method, with an active space composed of 4 electrons in the dx2 −y2 and dz2 orbitals of both nickel sites. The singlet and triplet states have strong multireference character, with each consisting of four configurations with almost equal weight (see table 3 of the supporting information for further details). Three nested embedding regions were investigated. Region A contains the active space only, region B further includes the remaining 6 d orbitals from the nickel ions, as well as the 4 valence orbitals of the bridging oxygen, and finally region

24


Page 24 of 61

Page 25 of 61

0.0 −10.0 −20.0 J / cm−1

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


−30.0 −40.0 −50.0

CASSCF → NEVPT2 NEVPT2 → icMRCCSD icMRCCSD → icMRCCSD(T) Experimental

−60.0 −70.0 −80.0 none

A

B C Correlated region

full

Figure 8: The super-exchange coupling constant for crystalline nickel oxide with increasing sized embedding region. The experimental result is taken from reference 161. C also includes 10 Ni-O σ-bonds from the next surrounding shell of oxygen atoms. Charge density plots of these regions are shown in figure 7, and table 3 details the size of their orbital spaces. Calculations on this system cannot utilise symmetry because of the point charge lattice, and we are therefore only able to compare to reference icMRCCSD calculations using the double-ζ split-valence polarization (def2-SVP) Karlsruhe basis set. 162 Convergence of the embedding procedure with respect to region size is demonstrated in figure 8. Clearly the bridging oxygen orbitals included in region B are critical for describing the interaction between the two nickel sites, in contrast to the surrounding O-Ni σ-bonds that do not appear to contribute significantly. Ab-initio calculations on the full cluster model (without the orbital embedding described in this paper) are known to correctly reproduce the ordering of the spin states, though frequently underestimate the magnitude of the J value. Our numerical value for J calculated using icMRCCSD on the whole cluster is −35.9 cm−1 , compared to the experimental value of −79 ± 1 cm−1 . 161,163 Table 4 details our values for J along with several taken from the literature, and a discussion of this discrepancy follows. Our CASSCF result is in general agreement with previous findings, 151,156 at −14.6 cm−1 . NEVPT2 gives −28.7 cm−1 , however our CASPT2 result (using the same basis) gives 25



Page 26 of 61

Table 4: Comparison of various bulk nickel oxide J values (in cm−1 ) from the present work and the literature. All calculations are for the full cluster model, without further orbital embedding. The cluster model itself, and the splitting energies used to calculate J, vary subtly between works (see relevant references for full details). Method CASSCF SA-CASSCFa SA-CASSCFb CASSCF+SOC NEVPT2 NEVPT2+SOCc CASPT2d CASPT2e,d CASPT2a,e CASPT2a,e,f icMRCCSD icMRCCSD+SOCc MCCEPAg,b Experiment

J

Reference

−14.6 −13.3 −16.0 −20.1 −28.7 −26.6 −37.4 −38.0 −47.2 −54.4 −35.9 −34.1 −29.2 −79 ± 1

This work 156 151 This work This work This work This work This work 156 156 This work This work 151 161 and 163

a

Uses a more flexible ANO basis. Uses approximately aug-V5ZPP basis. c Diagonal elements of the CASSCF SOC matrix replaced with higher level calculations. d Required a level shift of 0.1 to converge. e Ni 3s and 3p ‘semi-core’ orbitals included in dynamic correlation. f Uses larger active space, equivalent to our region B. g Dynamic correlation applied to 22 orbitals only, equivalent to our region B. b

−37.4 cm−1 , a result similar to even icMRCCSD. The two MRPT methods clearly have pathologically different behaviour for this system; in this case we would expect the CASPT2 to underestimate the splittings (and so the J value), because the ground state is a singlet. It is not clear if the performance of CASPT2 is legitimate or anomalous in this case. Fink and Staemmler calculated J using MCCEPA, an approximate form of MRCC (that they applied only to the equivalent of region B), and report a J value of −29.2 cm−1 , comparable to our own icMRCCSD result. 151 Graaf, Nieuwpoort and co-workers achieve a J value of −47.2 cm−1 using CASPT2 and a more flexible atomic natural orbital (ANO) basis set, and −54.4 cm−1 when using a larger active space equivalent to our region B. 156

26


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


This suggests that a larger basis set can account for perhaps 10 cm−1 at the correlated level. Extending the active space to include the bridging orbitals of region B may also account for another 7 cm−1 . Contributions from triple excitations might also be significant, as can be seen in embedding region B, where icMRCCSD(T) gives a J value of −38.6 cm−1 . Perhaps 3 cm−1 or more may be accounted for in the full model by the effect of higher excitations. Unfortunately we are unable to investigate these possibilities at the MRCC level, due to the prohibitive cost. We have also investigated the effect of spin-orbit coupling (SOC) on this system (see appendix A for full details). The CASSCF+SOC J value was −20.1 cm−1 , but the NEVPT2+SOC and icMRCCSD+SOC results were −26.6 and −34.1 respectively, and we therefore conclude that SOC is not important for this system. In summary, whilst the computations in this section demonstrate convergence of the J value from our embedded approach towards the icMRCCSD benchmark, the final estimate for J is still in error by a factor of 2 when compared to experiment. Cooperative effects arising from three or more nickel sites may account for the missing interactions, giving a better value for J, however, the present cluster model cannot account for this. A more extensive investigation will be the subject of future work.

4 4.1

Application to large Transition Metal Complexes Iron Complex FeC72 N2 H100

To demonstrate the applicability of this approach to much larger multireference systems, we calculate the first two excitation energies of the FeC72 N2 H100 complex. 31,32,164,165 The structure is shown in figure 9. The geometry and basis were taken from reference 31 (175 atoms, valence triple-ζ polarization (def2-TZVP) Karlsruhe basis set 162 with 2939 contracted functions). The central iron atom lies in a conjugated plane including two N p orbitals, and the π-system of the central ring of each ligand. The remaining four rings stand 27



Schematic

Region A

Region B

Region C

Region D

Region E

Figure 9: Schematic of the iron complex FeC27 N2 H100 , and charge density plots of its five embedding regions.

28


Page 28 of 61

Page 29 of 61

360.0

icMRCCSD → icMRCCSD(T) CASPT2 → icMRCCSD CASSCF → CASPT2

340.0 Excitation Energy / kJ mol−1

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


320.0

2nd Excitation

300.0 280.0 W2

260.0

X2

240.0 220.0

1st Excitation

200.0 180.0

W1

X1

160.0 None

A

B C D Correlated region

E

Full

Figure 10: Excitation energies of the iron complex. The 1st excitation is the quintet-triplet splitting, and the 2nd the quintet-singlet splitting. The PNO-CASPT2 method of reference 32 was used for the full system calculations. The labels W1 and W2 highlight our best embedded icMRCCSD results for the 1st and 2nd excitation energies respectively, and X1 and X2 similarly for icMRCCSD(T). orthogonal to this plane and almost face one another; however, one ring from each ligand is able to use its π system to coordinate with the iron atom in a σ like interaction. Although the complex possesses Ci point group symmetry, the non-embedded calculations were performed without symmetry adaption to allow for the use of density fitted integral transformations in the HF, 166 CASSCF, 167 and PNO-CASPT2 32 methods. The auxiliary basis sets recommended for def2-TZVPP were used for the density fitting. 168 To reduce the cost of state-specific CASSCF calculations, we follow the procedure of reference 32 for generating the reference molecular orbitals. Firstly, a restricted open-shell Hartree–Fock (ROHF) calculation is performed on the high spin quintet state (with four Fe d orbitals singly occupied), after which the doubly occupied space (excluding the core) is localised using IBOs. Nine orbitals local to the central Fe atom and its immediate N neighbours (corresponding to region C, discussed below) were subsequently optimised using CASSCF with an active space of six elections in the five Fe d orbitals. The remaining

29



environment orbitals were frozen at the ROHF/IBO level. We found it necessary to also freeze the two Fe-N sigma bonds when optimising the quintet state, to prevent the Fe d orbitals from falling out of the active space. The SCF energy is not affected in this particular case, because the quintet state is totally dominated by the ROHF configuration (with coefficient 1.000000), although such an orbital re-ordering would clearly affect a multireference correlated calculation. The singlet and triplet states of this complex have strong multireference character, consisting of 7 and 4 configurations respectively with coefficients of at least magnitude 0.1. Five nested embedded regions were investigated: region A contains the active space, region B further includes two additional Fe-N sigma bonds, region C two additional N p orbitals, and region D includes 4 additional π orbitals from the ligand central rings. Finally region E also includes two additional π orbitals from the non-central ligand rings, which coordinate directly to the Fe atom. All embedded regions are displayed in figure 9. Following the embedding procedure outlined in section 2.3 leads to prohibitively large virtual domains for regions D and E, that would require in excess of 250GB of RAM for the icMRCCSD calculations. We therefore reduced the number of host atoms involved in constructing these virtual spaces, and so decreased the number of virtual orbitals from 416 to 292 and from 601 to 354 for regions D and E respectively. (See the supporting information for full lists of host atoms for each region). Although this compromise may increase the virtual domain error in these cases, the effect will be small as the energy differences are converged with respect to embedding region size. Furthermore we expect cancellation of embedding artefacts to occur, as was the case with the carbene studied in section 3.2. An implementation of the local PNO-CASPT2 method within Molpro 32 was used to perform perturbative calculations on the whole system, with a level shift of 0.1 necessary to converge the solution. To render the calculation feasible, pair energies of less than 1.0 × 10−6 were approximated using the multipole approximation (Molpro keyword THRDIST). PNO and density-fitting domains were extended to include two shells of bonded atoms

30


Page 30 of 61

Page 31 of 61 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 5: Number of molecular orbitals for each region of the iron and nickel transition metal complexes. The percentage of the full doubly-occupied and virtual spaces are given in brackets. The inactive, active and virtual orbitals are treated at the higher level. Region

Environment

Inactive (%)

Active

Virtual (%)

FeC27 N2 H100 /def2-TZVP A 283 B 281 C 279 D 275 E 273 Full 83

0 (0) 2 (1) 4 (1) 8 (3) 10 (4) 200 (71)

5 5 5 5 5 5

44 (2) 106 (4) 106 (4) 292 (11) 354 (13) 2651 (100)

[NiC90 N20 H124 ]2+ /def2-TZVP A 409 B 403 C 399 Full 119

0 (0) 6 (1) 10 (2) 290 (71)

5 5 5 5

44 (1) 230 (6) 354 (9) 3761 (100)

Table 6: Comparison of excitation energies (in kJ mol−1 ) for the iron complex FeC72 N2 H100 . The final two rows refer to embedded icMRCC calculations labelled as W1, W2 and X1, X2 respectively in figure 10. The DLPNO-NEVPT2 result is based on SACASSCF orbitals, and the other results on the CASSCF orbitals detailed in the text. Method

1st Excitation

SA-CASSCF CASSCF DLPNO-NEVPT2 PNO-CASPT2 icMRCCSD in PNO-CASPT2

2nd Excitation

204.9 218.6 178.7 204.1

346.1 303.0

179.9 (W1) icMRCCSD(T) in icMRCCSD in PNO-CASPT2 176.6 (X1)

Reference 31 32 31 32

272.1 (W2)

-

262.4 (X2)

-

(Molpro keyword IEXT = 2), and a PNO occupation number threshold of 1.0 × 10−8 was used (Molpro keyword THRPNO_OCC). The lowest energy levels of this molecule are quintet, triplet and singlet spin states in ascending energy order. Figure 10 presents the first two vertical excitation energies of the complex, using CASPT2 embedded in CASSCF, icMRCCSD embedded in CASPT2, and finally icMRCCSD(T) embedded in icMRCCSD. In the latter, three layers of embedding 31



are used, with icMRCCSD(T) embedded in the largest possible icMRCCSD calculation (Region E), and the remainder of the molecule treated at the PNO-CASPT2 level. Both excitations show stable convergence of the energy with respect to the embedded region size, confirming the central assumption that the spin states of the metal atom are predominantly influenced by its immediate environment. The addition of certain orbitals to the embedding region affects the excitation energies differently, in particular, going from region B to D lowers the 2nd excitation energy but not the 1st. This happens because the extra orbitals in these regions enjoy significant overlap with the Fe dyz orbital, which is significantly populated in the open shell states but not in the singlet. The stabilisation subsequently cancels in the 1st excitation energy, but not in the 2nd. (We define the N-Fe-N vector as the z axis, and the plane containing the Fe atom and the ligand central rings as the xz plane). The perturbative triples correction in embedding region C has a larger effect on the second excited state for similar reasons. Table 6 compares our best embedded icMRCCSD results (labelled as W1 and W2 for the 1st and 2nd excitation energies respectively) and icMRCCSD(T) results (similarly labelled X1 and X2) with the linear scaling MRPT methods (the DLPNO-NEVPT2 result is taken from reference 31, and is based on CASSCF orbitals state-averaged (SA) over the quintet and triplet states). A detailed comparison between the CASSCF and SA-CASSCF energies, as well as between the PNO-CASPT2 and DLPNO-NEVPT2 energies can be found in reference 32, and is not reproduced here. Dynamic correlation energy is significant in this system, with the embedded MRCC calculations W1 and W2 changing the CASSCF excitation energies by -38.7 and -74 kJ mol−1 respectively. Furthermore the PNO-CASPT2 results are around 25-30 kJ mol−1 higher than W1 and W2, consistent with CASPT2 theory systematically underestimating the correlation energy of states with more paired electrons. Embedded icMRCCSD(T) calculations X1 and X2 reduce the icMRCCSD excitation energies by a further -3.3 and -9.7 kJ mol−1 respectively, however these energies are presumably not converged with respect

32


Page 32 of 61

Page 33 of 61 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


to the size of the embedding region. Perturbative triples calculations on larger regions unfortunately required too much RAM. The DLPNO-NEVPT2 1st excitation energy is very close to the embedded icMRCCSD result W1, albeit using different reference orbitals, though some orbital relaxation will have taken place with both methods due to the effect of singly excited determinants. To investigate this further, we have also performed an embedded calculation using NEVPT2 in region E and (not state averaged) CASSCF on the environment. This gives a 1st excitation energy of 178.4 (very close to both the DLPNO-NEVPT2 and embedded icMRCC results, suggesting that convergence has been reached with this region) and a 2nd excitation energy of 290.3 kJ mol−1 (nearly 20 kJ mol−1 higher than the embedded icMRCCSD result). These calculations suggest that although NEVPT2 theory may perform well for the 1st excitation energy, the same may not be true of the 2nd excitation energy. Here we also discuss the rationale behind choosing the embedding regions. For most systems studied in this work the choice of embedding regions is highly constrained by the point group symmetry of the molecule. The large iron complex, however, allows for many more possibilities when choosing regions. Initially, nested regions were constructed by adding groups of similar orbitals (both π and σ) with increasing distance from the multireference centre (the transition metal in this case). It was subsequently found that addition of the four next-nearest σ-bonds to region C (two N-H and two N-C bonds) had only a small impact on the excitation energies (the biggest change being +2.3 kJ mol−1 in the 2nd excitation energy). These σ bonds were subsequently omitted from regions D and E, allowing computational effort to be redeployed to the extended π system, which is seen to be important due to the strong planarity at the centre of the complex. For routine applications, an automatized approach to construct the embedding regions is certainly desirable and will be the subject of future work.

33



4.2

Nickel Complex [NiC90 N20 H120 ]2+

We now investigate a large nickel complex [NiC90 N20 H120 ]2+ taken from reference 169. The central nickel(II) metal atom is accompanied by two tri-dentate ligands (abbreviated as ‘tdta’ in reference 169) forming a distorted octahedral structure, with a triplet ground state. The structure is shown schematically in figure 11. This complex is a catalyst formed from a ‘click’ reaction, and has been previously studied using DLPNO-NEVPT2 to assign spectra and investigate spectroscopic properties. 31 The same geometry and def2-TZVP basis was used as in that work. Density fitting was used for the integral transformation, and HF, CASSCF and PNO-CASPT2 calculations, with the auxiliary basis sets recommended for def2-TZVPP used for the Ni atom, and those recommended for cc-pVDZ on all other atoms. 170 The reference molecular orbitals were calculated using an approach similar to that of the Fe complex in section 4.1. A triplet ROHF calculation was performed first and its doubly occupied orbitals IBO localised (excluding the core orbitals), before re-optimising the active space and those IBOs spatially nearby with state-specific CASSCF. We refer to these reference orbitals as semi-relaxed. An active space of 8 electrons in 5 Ni d orbitals was used, with 24 IBOs also optimised (those exerting a partial charge of at least 0.2 atomic units on the Ni or its six adjoining N atoms). The remaining inactive orbitals were frozen at the ROHF/IBO level. Similarly to the Fe complex, we found that the six Ni-N sigma bonds also had to be frozen in the high spin case, to prevent the Ni d orbitals from dropping out of the active space. The high-spin state was overwhelmingly dominated by the ROHF configuration, with a CI coefficient of 0.9998965. For comparison we also performed correlated calculations using a second set of unrelaxed reference orbitals, in which only the active space was optimised with CASSCF and all the inactive orbitals frozen at the ROHF/IBO level. In both cases, the singlet-state reference wave functions contained 5 configurations with coefficients of magnitude greater than 0.1 (see table 3 of the supporting information for full details). 34


Page 34 of 61

Page 35 of 61 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


Schematic

Region A

Region B

Region C

Figure 11: Structure and charge density plots of the embedding regions used in the nickel complex [NiC90 N20 H124 ]2+ .

35


∆EST / kJ mol−1

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

∆EST / kJ mol−1


Page 36 of 61

225.0 Unrelaxed: Active space only optimised 220.0 215.0 icMRCCSD → icMRCCSD(T) CASPT2 → icMRCCSD 210.0 CASSCF → CASPT2 205.0 200.0 195.0 190.0 185.0 Z1 Y1 180.0 225.0 220.0 Semi−relaxed: Active space + 24 IBOs optimised 215.0 icMRCCSD → icMRCCSD(T) CASPT2 → icMRCCSD 210.0 CASSCF → CASPT2 205.0 200.0 195.0 190.0 185.0 Z2 Y2 180.0 None A B C Full Correlated region

Figure 12: Singlet-triplet splitting ∆EST = ES − ET for the nickel complex. The reference orbitals are unrelaxed in the top panel, and semi-relaxed in the bottom panel. The PNOCASPT2 algorithm of reference 32 is used for the full calculations. Three embedding regions were considered. Region A contains the active space only, region B further includes the six surrounding Ni-N sigma bonds, and region C another four π orbitals on the ligand rings nearest the metal centre. Table 5 details the number orbitals in each region, which are also shown in figure 11. Perturbative calculations on the whole molecule were performed using the local PNO-CASPT2 program in Molpro with the same parameters as for the iron complex, though without any level shift. The vertical singlet-triplet splitting energy with increasing correlated region is shown in figure 12, for CASPT2 embedded in CASSCF, icMRCCSD embedded in CASPT2, and the triple layer case of icMRCCSD(T) embedded in icMRCCSD embedded in PNO-CASPT2. Our best icMRCCSD splitting energies are labelled Y1 and Y2 for the unrelaxed and semirelaxed reference orbitals respectively, and our best icMRCCSD(T) energies similarly labelled Z1 and Z2. Table 7 compares the embedded MRCC splitting energies to PNO-CASPT2 and DLPNO-NEVPT2, the latter of which is taken from reference 31 and which uses fully

36


Page 37 of 61 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 7: Comparison of excitation energies (in kJ mol−1 ) for the nickel complex [NiC90 N20 H120 ]2+ . The final two rows refer to embedded icMRCC calculations labelled in figure 12. The DLPNO-NEVPT2 result is based on fully relaxed SA-CASSCF orbitals, and the others on the CASSCF orbitals detailed in the text. Method

Unrelaxeda

Semi-relaxedb

Reference

207.0c 216.6 189.5c 184.4

31 32 31 32

184.4 (Y2)

-

182.6 (Z2)

-

SA-CASSCF CASSCF 220.4 DLPNO-NEVPT2 PNO-CASPT2 199.1 icMRCCSD in PNO-CASPT2 186.6 (Y1) icMRCCSD(T) in icMRCCSD in PNO-CASPT2 185.6 (Z1) a

Only the active space is optimised in the CASSCF. The active space and 24 IBOs are optimised in the CASSCF. c All reference orbitals are relaxed in SA-CASSCF. b

relaxed SA-CASSCF reference orbitals (see reference 169 for full details). For the case of semi-relaxed reference orbitals, the PNO-CASPT2, DLPNO-NEVPT2 and embedded MRCC results (Y2 and Z2) are all within a range of 6.9 kJ mol−1 , indicating that dynamic correlation is not very strong in this system, and that a perturbative treatment is adequate. The unrelaxed reference orbitals on the other hand produce more varied results, with PNO-CASPT2 giving a splitting 12.5 kJ mol−1 higher than Y1 and 13.5 kJ mol−1 higher than Z1. We attribute this to the explicit treatment of single excitations in icMRCC, which is better able to relax the reference orbitals than its perturbative counterpart, and can therefore compensate somewhat for less well optimised reference calculations.

5

Conclusions

We have presented an embedded icMRCC method, in which accurate icMRCC calculations are applied to the important target region of a molecule, and MRPT applied to the

37



Page 38 of 61

chemical environment. Using MRCC allows strongly correlated systems to be described properly, with both static and dynamic correlation energy included. Furthermore use of the internally contracted variant confers the desirable properties of size extensivity, invariance of the energy with respect to rotation of the active orbitals, and affordability when using moderate sized active spaces of five or more orbitals. The embedding scheme used here is based on partitioning of the orbital spaces of a CASSCF reference calculation, to construct a truncated subset of semi-canonical molecular orbitals in which the icMRCC theory is applied. Truncation of the virtual space is achieved by transforming PAOs resident on a set of ‘host atoms’ within the embedding region. The domain error introduced because of this was shown numerically to be less than 0.5 kJ mol−1 for benchmark systems. Simple subtractive embedding of the correlation energy was then used to apply MRPT to the environment region, and allow for multilayer embedding. The procedure is implemented in a development version of Molpro, interfaced to the general contraction code GeCCo. Benchmark calculations on small multireference systems have shown convergence of calculated energy differences with increasing embedded region size. The error due to the lower level treatment of the environment was found to be typically less than 5 kJ mol−1 (relative to icMRCCSD/cc-pVTZ on the whole system, when embedding within a MRPT environment); however, this error depends on the quality of the description of the environment. When using the subtractive embedding scheme, the embedded energies fall between the full MRPT and full MRCC results, provided that the embedding region is sufficiently large. Small embedding regions may lead to anomalies that make the result worse than MRPT (see region A of the large nickel complex, and regions A and B of the alkanes when embedding in NEVPT2). It is therefore necessary to confirm convergence of the calculated energy difference with respect to embedded region size. We have also presented multilayer embedding calculations using icMRCCSD(T) within icMRCCSD within PNO-CASPT2 for two large transition metal complexes with a triple-ζ

38


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


basis - FeC27 N2 H100 with 2939 basis functions and [NiC90 N20 H124 ]2+ with 4175 basis functions. Although we cannot measure the true embedding error for these large systems, energy differences do show convincing convergence, implying that a perturbative description of the wider environment is adequate. These calculations demonstrate that strongly correlated systems of interesting size can be treated with the approach outlined here, at attainable cost, and to higher accuracy than MRPT alone.

Supporting Information Available Supporting Information Available: All the energies used to make the plots, Cartesian geometries of all the molecules studied, lists of host atoms used for the embedded calculations, J values calculated from individual spin splittings for the nickel oxide cluster model, significant configurations and their coefficients in the CASSCF reference wave functions of each system, and details of the quartic polynomial functions fitted to the [Fe(H2 O)6 ]2+ singlet and quintet states. This material is available free of charge via the Internet at http://pubs.acs.org.

Acknowledgement D.J.C. thanks Dr. Qianli Ma for helpful discussions. H.-J.W. acknowledges support from the European Research Council (ERC), Advanced Grant 320723 (ASES).

A

Details of the Spin-orbit Coupling Calculations for Nickel Oxide

The effect of SOC was calculated at the CASSCF level using the perturbative Breit–Pauli SOC operator. 171 The full electronic plus SOC Hamiltonian was diagonalised in the basis

39



of electronic states. SA-CASSCF reference orbitals were used, with a larger active space consisting of the 10 d orbitals from both nickel atoms. The lowest set of singlet, triplet and quintet states cannot interact via the SOC operator due to the point group symmetry of the cluster model. We therefore also include the next six singlet, triplet and quintet states, to give a total of 21 interacting states, which allows the lowest set to interact indirectly. J values were then calculated using the lowest energy singlet, the middle state of the formerly degenerate lowest energy triplet term, and similarly for the quintet term. The spin orbit interaction lowered the J value from −14.4 to −20.1 cm−1 at the CASSCF level (with both calculations using the same aforementioned reference orbitals). Dynamic correlation effects can be approximately incorporated by substituting the diagonal elements of the electronic plus SOC Hamiltonian matrix, with previously calculated NEVPT2 and icMRCCSD energies. We do this only for the lowest singlet, triplet and quintet states, using the reference orbitals and smaller active space discussed in the main body of the text. The energies of the higher states were then shifted, so that their excitation energies relative to the ground state were unchanged. Following this procedure with NEVPT2 energies raises the J value from −28.7 to −26.6 cm−1 , and with icMRCCSD also raises the J value from −36.0 to −34.1 cm−1 . We conclude therefore that the spin-orbit interaction is small for this system and cannot account for the discrepancy with the experimentally measured value of J.

References (1) Cramer, C. J.; Włoch, M.; Piecuch, P.; Puzzarini, C.; Gagliardi, L. Theoretical models on the Cu2 O2 torture track: Mechanistic implications for oxytyrosinase and small-molecule analogues. J. Phys. Chem. A 2006, 110, 1991–2004. (2) Pierloot, K.; Zhao, H.; Vancoillie, S. Copper corroles: the question of noninnocence. Inorg. Chem. 2010, 49, 10316–10329. 40


Page 40 of 61

Page 41 of 61 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


(3) Vancoillie, S.; Zhao, H.; Radon, M.; Pierloot, K. Performance of CASPT2 and DFT for relative spin-state energetics of heme models. J. Chem. Theory Comput. 2010, 6, 576–582. (4) Rohrmüller, M.; Hoffmann, A.; Thierfelder, C.; Herres-Pawlis, S.; Schmidt, W. G. The Cu2 O2 torture track for a real-life system:[Cu2 (btmgp)2 O2 ]2+ oxo and peroxo species in density functional calculations. J. Comput. Chem. 2015, 36, 1672–1685. (5) Singh, S. K.; Rajaraman, G. Deciphering the origin of giant magnetic anisotropy and fast quantum tunnelling in Rhenium (IV) single-molecule magnets. Nat. Commun. 2016, 7 . (6) Neese, F.; Pantazis, D. A. What is not required to make a single molecule magnet. Faraday Discuss. 2011, 148, 229–238. (7) Chowdhury, S. R.; Mishra, S. Heavy Ligand Atom Induced Large Magnetic Anisotropy in Mn (II) Complexes. Phys. Chem. Chem. Phys. 2017, (8) Yang, K. R.; Jalan, A.; Green, W. H.; Truhlar, D. G. Which ab initio wave function methods are adequate for quantitative calculations of the energies of biradicals? The performance of coupled-cluster and multi-reference methods along a single-bond dissociation coordinate. J. Chem. Theory Comput. 2012, 9, 418–431. (9) Robinson, D. Splitting multiple bonds: A comparison of methodologies on the accuracy of bond dissociation energies. J. Comput. Chem. 2013, 34, 2625–2634. (10) Kedziora, G. S.; Barr, S. A.; Berry, R.; Moller, J. C.; Breitzman, T. D. Bond breaking in stretched molecules: multi-reference methods versus density functional theory. Theor. Chem. Acc. 2016, 135, 79. (11) Zhang, D.; Liu, C. Electronic Structures of Anti-Ferromagnetic Tetraradicals: Ab Initio and Semi-Empirical Studies. J. Chem. Theory Comput. 2016, 12, 1714–1727. 41



(12) Roca-Sanjuán, D.; Aquilante, F.; Lindh, R. Multiconfiguration second-order perturbation theory approach to strong electron correlation in chemistry and photochemistry. WIREs Comput. Mol. Sci. 2012, 2, 585–603. (13) Marazzi, M.; Gattuso, H.; Monari, A. Nile blue and Nile red optical properties predicted by TD-DFT and CASPT2 methods: static and dynamic solvent effects. Theor. Chem. Acc. 2016, 135, 57. (14) Liu, L.; Liu, J.; Martinez, T. J. Dynamical correlation effects on photoisomerization: Ab initio multiple spawning dynamics with MS-CASPT2 for a model transprotonated schiff base. J. Phys. Chem. B 2016, 120, 1940–1949. (15) Chung, L. W.; Hayashi, S.; Lundberg, M.; Nakatsu, T.; Kato, H.; Morokuma, K. Mechanism of efficient firefly bioluminescence via adiabatic transition state and seam of sloped conical intersection. J. Am. Chem. Soc. 2008, 130, 12880–12881. (16) Chen, S.-F.; Liu, Y.-J.; Navizet, I.; Ferré, N.; Fang, W.-H.; Lindh, R. Systematic theoretical investigation on the light emitter of firefly. J. Chem. Theory Comput. 2011, 7, 798–803. (17) Werner, H.-J.; Schütz, M. An efficient local coupled cluster method for accurate thermochemistry of large systems. J. Chem. Phys. 2011, 135, 144116. (18) Yang, J.; Chan, G. K.-L.; Manby, F. R.; Schütz, M.; Werner, H.-J. The orbitalspecific-virtual local coupled cluster singles and doubles method. J. Chem. Phys. 2012, 136, 144105. (19) Riplinger, C.; Neese, F. An efficient and near linear scaling pair natural orbital based local coupled cluster method. J. Chem. Phys. 2013, 138, 034106. (20) Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E. F.; Neese, F. Sparse maps - A systematic infrastructure for reduced-scaling electronic structure methods. II. Linear 42


Page 42 of 61

Page 43 of 61 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


scaling domain based pair natural orbital coupled cluster theory. J. Chem. Phys. 2016, 144, 024109. (21) Saitow, M.; Becker, U.; Riplinger, C.; Valeev, E. F.; Neese, F. A new near-linear scaling, efficient and accurate, open-shell domain-based local pair natural orbital coupled cluster singles and doubles theory. J. Chem. Phys. 2017, 146, 164105. (22) Schwilk, M.; Ma, Q.; Köppl, C.; Werner, H.-J. Scalable electron correlation methods. 3. Efficient and accurate parallel local coupled cluster with pair natural orbitals (PNO-LCCSD). J. Chem. Theory Comput. 2017, 13, 3650–3675. (23) Ma, Q.; Schwilk, M.; Köppl, C.; Werner, H.-J. Scalable Electron Correlation Methods. 4. Parallel Explicitly Correlated Local Coupled Cluster with Pair Natural Orbitals (PNO-LCCSD-F12). J. Chem. Theory Comput. 2017, (24) Edmiston, C.; Krauss, M. Pseudonatural Orbitals as a Basis for the Superposition of Configurations. I. He2+ . J. Chem. Phys. 1966, 45, 1833–1839. (25) Meyer, W. Ionization energies of water from PNO-CI calculations. Int. J. Quantum Chem. 1971, 5, 341–348. (26) Meyer, W. PNO–CI Studies of electron correlation effects. I. Configuration expansion by means of nonorthogonal orbitals, and application to the ground state and ionized states of methane. J. Chem. Phys. 1973, 58, 1017–1035. (27) Werner, H.-J.; Meyer, W. PNO-CI and PNO-CEPA studies of electron correlation effects: V. Static dipole polarizabilities of small molecules. Mol. Phys. 1976, 31, 855–872. (28) Jungen, M.; Ahlrichs, R. Ab initio calculations on small hydrides including electron correlation. Theor. Chim. Acta 1970, 17, 339–347.

43



(29) Neese, F.; Hansen, A.; Liakos, D. G. Efficient and accurate approximations to the local coupled cluster singles doubles method using a truncated pair natural orbital basis. J. Chem. Phys. 2009, 131, 064103. (30) Neese, F.; Wennmohs, F.; Hansen, A. Efficient and accurate local approximations to coupled-electron pair approaches: An attempt to revive the pair natural orbital method. J. Chem. Phys. 2009, 130, 114108. (31) Guo, Y.; Sivalingam, K.; Valeev, E. F.; Neese, F. Sparse Maps - A systematic infrastructure for reduced-scaling electronic structure methods. III. Linear-scaling multireference domain-based pair natural orbital N-electron valence perturbation theory. J. Chem. Phys. 2016, 144, 094111. (32) Menezes, F.; Kats, D.; Werner, H.-J. Local complete active space second-order perturbation theory using pair natural orbitals (PNO-CASPT2). J. Chem. Phys. 2016, 145, 124115. (33) Chaudhuri, R. K.; Freed, K. F.; Hose, G.; Piecuch, P.; Kowalski, K.; Włoch, M.; Chattopadhyay, S.; Mukherjee, D.; Rolik, Z.; Szabados, Á.; Tóth, G.; Surján, P. R. Comparison of low-order multireference many-body perturbation theories. J. Chem. Phys. 2005, 122, 134105. (34) Lang, J.; Švaňa, M.; Demel, O.; Brabec, J.; Kedžuch, S.; Noga, J.; Kowalski, K.; Pittner, J. A MRCC study of the isomerisation of cyclopropane. Mol. Phys. 2017, 1–12. (35) Demel, O.; Pittner, J.; Neese, F. A Local Pair Natural Orbital-Based Multireference Mukherjee’s Coupled Cluster Method. J. Chem. Theory Comput. 2015, 11, 3104– 3114. (36) Köhn, A.; Hanauer, M.; Mück, L. A.; Jagau, T.-C.; Gauss, J. State-specific multireference coupled-cluster theory. WIREs Comput. Mol. Sci. 2013, 3, 176–197. 44


Page 44 of 61

Page 45 of 61 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


(37) Mahapatra, U. S.; Datta, B.; Mukherjee, D. A size-consistent state-specific multireference coupled cluster theory: Formal developments and molecular applications. J. Chem. Phys. 1999, 110, 6171–6188. (38) Evangelista, F. A.; Allen, W. D.; Schaefer III, H. F. High-order excitations in stateuniversal and state-specific multireference coupled cluster theories: Model systems. J. Chem. Phys. 2006, 125, 154113. (39) Gomes, A. S. P.; Jacob, C. R. Quantum-chemical embedding methods for treating local electronic excitations in complex chemical systems. Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 2012, 108, 222–277. (40) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227–249. (41) Senn, H. M.; Thiel, W. QM/MM methods for biomolecular systems. Angew. Chem., Int. Ed. Engl. 2009, 48, 1198–1229. (42) Anderson, J. A.; Hopkins, B. W.; Chapman, J. L.; Tschumper, G. S. A systematic assessment of density functionals and ONIOM schemes for the study of hydrogen bonding between water and the side chains of serine, threonine, asparagine, and glutamine. J. Mol. Struct.: THEOCHEM 2006, 771, 65–71. (43) Roca-Sanjuán, D.; Olaso-González, G.; Coto, P. B.; Merchán, M.; SerranoAndrés, L. Modeling hole transfer in DNA. II. Molecular basis of charge transport in the DNA chain. Theor. Chem. Acc. 2010, 126, 177–183. (44) Manathunga, M.; Yang, X.; Luk, H. L.; Gozem, S.; Frutos, L. M.; Valentini, A.; Ferrè, N.; Olivucci, M. Probing the photodynamics of rhodopsins with reduced retinal chromophores. J. Chem. Theory Comput. 2016, 12, 839–850.

45



(45) Maseras, F.; Morokuma, K. IMOMM: A new integrated ab initio + molecular mechanics geometry optimization scheme of equilibrium structures and transition states. J. Comput. Chem. 1995, 16, 1170–1179. (46) Chung, L. W.; Sameera, W.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z.; Liu, F.; Li, H.-B.; Ding, L.; Morokuma, K. The ONIOM method and its applications. Chem. Rev. 2015, 115, 5678– 5796. (47) Vreven, T.; Byun, K. S.; Komáromi, I.; Dapprich, S.; Montgomery Jr, J. A.; Morokuma, K.; Frisch, M. J. Combining quantum mechanics methods with molecular mechanics methods in ONIOM. J. Chem. Theory Comput. 2006, 2, 815–826. (48) Wesolowski, T. A.; Warshel, A. Frozen density functional approach for ab initio calculations of solvated molecules. J. Phys. Chem. 1993, 97, 8050–8053. (49) Goodpaster, J. D.; Barnes, T. A.; Manby, F. R.; Miller III, T. F. Density functional theory embedding for correlated wavefunctions: Improved methods for open-shell systems and transition metal complexes. J. Chem. Phys. 2012, 137, 224113. (50) Hégely, B.; Nagy, P. R.; Ferenczy, G. G.; Kállay, M. Exact density functional and wave function embedding schemes based on orbital localization. J. Chem. Phys. 2016, 145, 064107. (51) Libisch, F.; Huang, C.; Carter, E. A. Embedded correlated wavefunction schemes: theory and applications. Acc. Chem. Res. 2014, 47, 2768–2775. (52) Jacob, C. R.; Neugebauer, J. Subsystem density-functional theory. WIREs Comput. Mol. Sci. 2014, 4, 325–362. (53) Manby, F. R.; Stella, M.; Goodpaster, J. D.; Miller III, T. F. A simple, exact

46


Page 46 of 61

Page 47 of 61 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


density-functional-theory embedding scheme. J. Chem. Theory Comput. 2012, 8, 2564–2568. (54) Goodpaster, J. D.; Barnes, T. A.; Manby, F. R.; Miller III, T. F. Accurate and systematically improvable density functional theory embedding for correlated wavefunctions. J. Chem. Phys. 2014, 140, 18A507. (55) Stella, M.; Bennie, S. J.; Manby, F. R. Computational study of adsorption of cobalt on benzene and coronene. Mol. Phys. 2015, 113, 1858–1864. (56) Bennie, S. J.; Stella, M.; Miller III, T. F.; Manby, F. R. Accelerating wavefunction in density-functional-theory embedding by truncating the active basis set. J. Chem. Phys. 2015, 143, 024105. (57) Barnes, T. A.; Goodpaster, J. D.; Manby, F. R.; Miller III, T. F. Accurate basis set truncation for wavefunction embedding. J. Chem. Phys. 2013, 139, 024103. (58) Bennie, S. J.; van der Kamp, M. W.; Pennifold, R. C.; Stella, M.; Manby, F. R.; Mulholland, A. J. A projector-embedding approach for multiscale coupled-cluster calculations applied to citrate synthase. J. Chem. Theory Comput. 2016, 12, 2689– 2697. (59) Knizia, G.; Chan, G. K.-L. Density matrix embedding: A simple alternative to dynamical mean-field theory. Phys. Rev. Lett. 2012, 109, 186404. (60) Knizia, G.; Chan, G. K.-L. Density matrix embedding: A strong-coupling quantum embedding theory. J. Chem. Theory Comput. 2013, 9, 1428–1432. (61) Wouters, S.; Jiménez-Hoyos, C. A.; Sun, Q.; Chan, G. K.-L. A practical guide to density matrix embedding theory in quantum chemistry. J. Chem. Theory Comput. 2016, 12, 2706–2719.

47



(62) Booth, G. H.; Chan, G. K.-L. Spectral functions of strongly correlated extended systems via an exact quantum embedding. Phys. Rev. B 2015, 91, 155107. (63) Bulik, I. W.; Chen, W.; Scuseria, G. E. Electron correlation in solids via density embedding theory. J. Chem. Phys. 2014, 141, 054113. (64) Zheng, B.-X.; Chan, G. K.-L. Ground-state phase diagram of the square lattice Hubbard model from density matrix embedding theory. Phys. Rev. B 2016, 93, 035126. (65) Gutdeutsch, U.; Birkenheuer, U.; Krüger, S.; Rösch, N. On cluster embedding schemes based on orbital space partitioning. J. Chem. Phys. 1997, 106, 6020–6030. (66) Whitten, J.; Yang, H. Theoretical studies of surface reactions on metals. Int. J. Quantum Chem. 1995, 56, 41–47. (67) Buenker, R.; Liebermann, H.; Kokh, D.; Izgorodina, E.; Whitten, J. Configuration interaction study of the excited states of CO adsorbed on a Pt97 cluster. Chem. Phys. 2003, 291, 115–124. (68) Henderson, T. M. Embedding wave function theory in density functional theory. J. Chem. Phys. 2006, 125, 014105. (69) Mata, R. A.; Werner, H.-J.; Schütz, M. Correlation regions within a localized molecular orbital approach. J. Chem. Phys. 2008, 128, 144106. (70) Li, W.; Piecuch, P. Multilevel Extension of the Cluster-in-Molecule Local Correlation Methodology: Merging Coupled-Cluster and Møller- Plesset Perturbation Theories. J. Phys. Chem. A 2010, 114, 6721–6727. (71) Rolik, Z.; Kállay, M. A general-order local coupled-cluster method based on the cluster-in-molecule approach. J. Chem. Phys. 2011, 135, 104111.

48


Page 48 of 61

Page 49 of 61 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


(72) Myhre, R. H.; Sánches de Merás, A. M.; Koch, H. The extended CC2 model ECC2. Mol. Phys. 2013, 111, 1109–1118. (73) Myhre, R. H.; Sánchez de Merás, A. M.; Koch, H. Multi-level coupled cluster theory. J. Chem. Phys. 2014, 141, 224105. (74) Myhre, R. H.; Koch, H. The multilevel CC3 coupled cluster model. J. Chem. Phys. 2016, 145, 044111. (75) Myhre, R. H.; Coriani, S.; Koch, H. Near-Edge X-ray Absorption Fine Structure within Multilevel Coupled Cluster Theory. J. Chem. Theory Comput. 2016, 12, 2633–2643. (76) Høyvik, I.-M.; Myhre, R. H.; Koch, H. Correlated natural transition orbitals for core excitation energies in multilevel coupled cluster models. J. Chem. Phys. 2017, 146, 144109. (77) Stoll, H. On the correlation energy of graphite. J. Chem. Phys. 1992, 97, 8449– 8454. (78) Paulus, B. The method of increments - a wavefunction-based ab initio correlation method for solids. Phys. Rev. 2006, 428, 1–52. (79) Friedrich, J.; Dolg, M. Fully automated incremental evaluation of MP2 and CCSD (T) energies: Application to water clusters. J. Chem. Theory Comput. 2009, 5, 287– 294. (80) Friedrich, J.; Hanrath, M.; Dolg, M. Fully automated implementation of the incremental scheme: Application to CCSD energies for hydrocarbons and transition metal compounds. J. Chem. Phys. 2007, 126, 154110.

49



(81) Friedrich, J.; Hanrath, M.; Dolg, M. Error analysis of incremental electron correlation calculations and applications to clusters and potential energy surfaces. Chem. Phys. 2007, 338, 33–43. (82) Friedrich, J.; Dolg, M. Implementation and performance of a domain-specific basis set incremental approach for correlation energies: Applications to hydrocarbons and a glycine oligomer. J. Chem. Phys. 2008, 129, 244105. (83) Friedrich, J.; Dolg, M.; Gansäuer, A.; Geich-Gimbel, D.; Lauterbach, T. A combined theoretical and experimental study of efficient and fast titanocene-catalyzed 3-exo cyclizations. J. Am. Chem. Soc. 2005, 127, 7071–7077. (84) Banerjee, A.; Simons, J. The coupled-cluster method with a multiconfiguration reference state. Int. J. Quantum Chem. 1981, 19, 207–216. (85) Evangelista, F. A.; Gauss, J. An orbital-invariant internally contracted multireference coupled cluster approach. J. Chem. Phys. 2011, 134, 114102. (86) Hanauer, M.; Köhn, A. Pilot applications of internally contracted multireference coupled cluster theory, and how to choose the cluster operator properly. J. Chem. Phys. 2011, 134, 204111. (87) Hanauer, M.; Köhn, A. Communication: Restoring full size extensivity in internally contracted multireference coupled cluster theory. J. Chem. Phys. 2012, (88) Hanauer, M. Internally contracted multireference coupled-cluster methods. Ph.D. thesis, Mainz, Univ., Diss., 2013; see http://nbn-resolving.org/urn:nbn:de:hebis:7734238. (89) Aoto, Y. A.; Köhn, A. Internally contracted multireference coupled-cluster theory in a multistate framework. J. Chem. Phys. 2016, 144, 074103.

50


Page 50 of 61

Page 51 of 61 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


(90) Roos, B. O. The Complete Active Space Self-Consistent Field Method and its Applications in Electronic Structure Calculations. Adv. Chem. Phys. 2007, 69, 399–445. (91) Pulay, P. Localizability of dynamic electron correlation. Chem. Phys. Lett. 1983, 100, 151–154. (92) Hanauer, M.; Köhn, A. Perturbative treatment of triple excitations in internally contracted multireference coupled cluster theory. J. Chem. Phys. 2012, 136, 204107. (93) Aoto, Y. A.; Köhn, A. Revisiting the F + HCl → HF + Cl reaction using a multireference coupled-cluster method. Phys. Chem. Chem. Phys. 2016, 18, 30241– 30253. (94) Aoto, Y. A.; de Lima Batista, A. P.; Köhn, A.; de Oliveira-Filho, A. G. How to arrive at accurate benchmark values for transition metal compounds: Computation or experiment? J. Chem. Theory Comput. 2017, (95) Samanta, P. K.; Mukherjee, D.; Hanauer, M.; Köhn, A. Excited states with internally contracted multireference coupled-cluster linear response theory. J. Chem. Phys. 2014, 140, 134108. (96) Werner, H.-J. Third-order multireference perturbation theory The CASPT3 method. Mol. Phys. 1996, 89, 645–661. (97) Celani, P.; Werner, H.-J. Multireference perturbation theory for large restricted and selected active space reference wave functions. J. Chem. Phys. 2000, 112, 5546– 5557. (98) Ghigo, G.; Roos, B. O.; Malmqvist, P.-Å. A modified definition of the zeroth-order Hamiltonian in multiconfigurational perturbation theory (CASPT2). Chem. Phys. Lett. 2004, 396, 142–149.

51



(99) Roos, B. O.; Andersson, K. Multiconfigurational perturbation theory with level shift - the Cr2 potential revisited. Chem. Phys. Lett. 1995, 245, 215–223. (100) Roos, B. O.; Andersson, K.; Fülscher, M. P.; Serrano-Andrés, L.; Pierloot, K.; Merchán, M.; Molina, V. Applications of level shift corrected perturbation theory in electronic spectroscopy. J. Mol. Struct.: THEOCHEM 1996, 388, 257–276. (101) Forsberg, N.; Malmqvist, P.-Å. Multiconfiguration perturbation theory with imaginary level shift. Chem. Phys. Lett. 1997, 274, 196–204. (102) Angeli, C.; Cimiraglia, R.; Evangelisti, S.; Leininger, T.; Malrieu, J.-P. Introduction of n-electron valence states for multireference perturbation theory. J. Chem. Phys. 2001, 114, 10252–10264. (103) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. N-electron valence state perturbation theory: a fast implementation of the strongly contracted variant. Chem. Phys. Lett. 2001, 350, 297–305. (104) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. n-electron valence state perturbation theory: A spinless formulation and an efficient implementation of the strongly contracted and of the partially contracted variants. J. Chem. Phys. 2002, 117, 9138– 9153. (105) Neese, F. The ORCA program system. WIREs Comput. Mol. Sci. 2012, 2, 73–78. (106) Menezes, F. A local complete active space 2nd-order perturbation theory method. Ph.D. thesis, Univ. Stuttgart, 2016; see http://dx.doi.org/10.18419/opus-8991. (107) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: a general-purpose quantum chemistry program package. WIREs Comput. Mol. Sci. 2012, 2, 242–253.

52


Page 52 of 61

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


(108) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Györffy, W.; Kats, D.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; Shamasundar, K. R.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hesselmann, A.; Hetzer, G.; Hrenar, T.; Jansen, G.; Köppl, C.; Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; O’Neill, D. P.; Palmieri, P.; Peng, D.; Pflüger, K.; Pitzer, R.; Reiher, M.; Shiozaki, T.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M. MOLPRO, version 2017.0, a package of ab initio programs; 2017; see http://www.molpro.net. (109) Knizia, G. Intrinsic atomic orbitals: An unbiased bridge between quantum theory and chemical concepts. J. Chem. Theory Comput. 2013, 9, 4834–4843. (110) Lykos, P. G.; Parr, R. G. On the Pi-Electron Approximation and Its Possible Refinement. J. Chem. Phys. 1956, 24, 1166–1173. (111) Hosteny, R.; Dunning Jr, T.; Gilman, R.; Pipano, A.; Shavitt, I. Ab initio study of the π-electron states of trans-butadiene. J. Chem. Phys. 1975, 62, 4764–4779. (112) Sæbø, S.; Pulay, P. Local configuration interaction: An efficient approach for larger molecules. Chem. Phys. Lett. 1985, 113, 13–18. (113) Sæbø, S.; Pulay, P. Fourth-order Møller–Plessett perturbation theory in the local correlation treatment. I. Method. J. Chem. Phys. 1987, 86, 914–922. (114) Sæbø, S.; Pulay, P. The local correlation treatment. II. Implementation and tests. J. Chem. Phys. 1988, 88, 1884–1890. (115) Saebø, S.; Pulay, P. Local treatment of electron correlation. Annu. Rev. Phys. Chem. 1993, 44, 213–236.

53



(116) Hampel, C.; Werner, H.-J. Local treatment of electron correlation in coupled cluster theory. J. Chem. Phys. 1996, 104, 6286–6297. (117) Schütz, M.; Werner, H.-J. Low-order scaling local electron correlation methods. IV. Linear scaling local coupled-cluster (LCCSD). J. Chem. Phys. 2001, 114, 661–681. (118) Werner, H.-J.; Pflüger, K. On the selection of domains and orbital pairs in local correlation treatments. Annu. Rep. Comput. Chem. 2006, 2, 53–80. (119) Werner, H.-J.; Knizia, G.; Krause, C.; Schwilk, M.; Dornbach, M. Scalable electron correlation methods I.: PNO-LMP2 with linear scaling in the molecular size and near-inverse-linear scaling in the number of processors. J. Chem. Theory Comput. 2015, 11, 484–507. (120) Ma, Q.; Werner, H.-J. Scalable Electron Correlation Methods. 2. Parallel PNOLMP2-F12 with Near Linear Scaling in the Molecular Size. J. Chem. Theory Comput. 2015, 11, 5291–5304. (121) Lindh, R.; Ryu, U.; Liu, B. The reduced multiplication scheme of the Rys quadrature and new recurrence relations for auxiliary function based two-electron integral evaluation. J. Chem. Phys. 1991, 95, 5889–5897. (122) Knowles, P.; Handy, N. A new determinant-based full configuration interaction method. Chem. Phys. Lett. 1984, 111, 315–321. (123) Knowles, P. J.; Handy, N. C. A determinant based full configuration interaction program. Comput. Phys. Commun. 1989, 54, 75–83. (124) Werner, H.-J.; Manby, F. R.; Knowles, P. J. Fast linear scaling second-order MøllerPlesset perturbation theory (MP2) using local and density fitting approximations. J. Chem. Phys. 2003, 118, 8149–8160.

54


Page 54 of 61

Page 55 of 61 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


(125) Schütz, M.; Manby, F. R. Linear scaling local coupled cluster theory with density fitting. Part I: 4-external integrals. Phys. Chem. Chem. Phys. 2003, 5, 3349–3358. (126) Werner, H.-J.; Knowles, P. J. A second order multiconfiguration SCF procedure with optimum convergence. J. Chem. Phys. 1985, 82, 5053–5063. (127) Knowles, P. J.; Werner, H.-J. An efficient second-order MC SCF method for long configuration expansions. Chem. Phys. Lett. 1985, 115, 259–267. (128) Werner, H.-J.; Meyer, W. A quadratically convergent multiconfiguration–selfconsistent field method with simultaneous optimization of orbitals and CI coefficients. J. Chem. Phys. 1980, 73, 2342–2356. (129) Werner, H.-J.; Meyer, W. A quadratically convergent MCSCF method for the simultaneous optimization of several states. J. Chem. Phys. 1981, 74, 5794–5801. (130) Johnson III, R. D. NIST Computational Chemistry Comparison and Benchmark Database, Release 18 ; NIST Standard Reference Database Number 101: http://cccbdb.nist.gov/, 2016. (131) Dunning Jr, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (132) Woon, D. E.; Dunning Jr, T. H. Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. J. Chem. Phys. 1993, 98, 1358–1371. (133) Balabanov, N. B.; Peterson, K. A. Systematically convergent basis sets for transition metals. I. All-electron correlation consistent basis sets for the 3 d elements Sc–Zn. J. Chem. Phys. 2005, 123, 064107. (134) Avogadro: an open-source molecular builder and visualization tool. Version 1.20. see http://avogadro.cc. 55



(135) Hanwell, M. D.; Curtis, D. E.; Lonie, D. C.; Vandermeersch, T.; Zurek, E.; Hutchison, G. R. Avogadro: an advanced semantic chemical editor, visualization, and analysis platform. J. Cheminf. 2012, 4, 17. (136) Andersson, K.; Roos, B. O. Multiconfigurational second-order perturbation theory: A test of geometries and binding energies. Int. J. Quantum Chem. 1993, 45, 591– 607. (137) Malmqvist, P.-Å.; Roos, B. O. Inclusion of dynamic σ-π polarization in π-electron ab initio calculations. Theor. Chim. Acta 1992, 83, 191–199. (138) Davidson, E. R. The spatial extent of the v state of ethylene and its relation to dynamic correlation in the cope rearrangement. J. Phys. Chem. 1996, 100, 6161–6166. (139) Angeli, C. On the nature of the π → π* ionic excited states: The V state of ethene as a prototype. J. Comput. Chem. 2009, 30, 1319–1333. (140) Gozem, S.; Huntress, M.; Schapiro, I.; Lindh, R.; Granovsky, A. A.; Angeli, C.; Olivucci, M. Dynamic electron correlation effects on the ground state potential energy surface of a retinal chromophore model. J. Chem. Theory Comput. 2012, 8, 4069. (141) Schapiro, I.; Sivalingam, K.; Neese, F. Assessment of n-electron valence state perturbation theory for vertical excitation energies. J. Chem. Theory Comput. 2013, 9, 3567–3580. (142) Herrmann, W. A.; Koecher, C. N-Heterocyclic Carbenes. Angew. Chem., Int. Ed. Engl. 1997, 36, 2162–2187. (143) Hopkinson, M. N.; Richter, C.; Schedler, M.; Glorius, F. An overview of Nheterocyclic carbenes. Nature 2014, 510, 485.

56


Page 56 of 61

Page 57 of 61 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


(144) Díez-González, S.; Marion, N.; Nolan, S. P. N-heterocyclic carbenes in late transition metal catalysis. Chem. Rev. 2009, 109, 3612–3676. (145) Goedecke, C.; Leibold, M.; Siemeling, U.; Frenking, G. When does carbonylation of carbenes yield ketenes? A theoretical study with implications for synthesis. J. Am. Chem. Soc. 2011, 133, 3557–3569. (146) Miliordos, E.; Xantheas, S. S. An accurate and efficient computational protocol for obtaining the complete basis set limits of the binding energies of water clusters at the MP2 and CCSD(T) levels of theory: Application to (H2 O)m , m= 2-6, 8, 11, 16, and 17. J. Chem. Phys. 2015, 142, 234303. (147) Andersson, K.; Malmqvist, P.; Roos, B. O. Second-order perturbation theory with a complete active space self-consistent field reference function. J. Chem. Phys. 1992, 96, 1218–1226. (148) Anderson, P. Antiferromagnetism. Theory of superexchange interaction. Phys. Rev. 1950, 79, 350. (149) Illas, F.; Moreira, I. P.; De Graaf, C.; Barone, V. Magnetic coupling in biradicals, binuclear complexes and wide-gap insulators: a survey of ab initio wave function and density functional theory approaches. Theor. Chem. Acc. 2000, 104, 265–272. (150) Maki, G. Ligand Field Theory of Ni (II) Complexes. I. Electronic Energies and Singlet Ground-State Conditions of Ni (II) Complexes of Different Symmetries. J. Chem. Phys. 1958, 28, 651–662. (151) Fink, K.; Staemmler, V. A modified CAS-CI approach for an efficient calculation of magnetic exchange coupling constants. Mol. Phys. 2013, 111, 2594–2605. (152) de PR Moreira, I.; Illas, F.; Martin, R. L. Effect of Fock exchange on the electronic structure and magnetic coupling in NiO. Phys. Rev. B 2002, 65, 155102. 57



(153) Wang, C.; Fink, K.; Staemmler, V. A quantum chemical ab initio study of the superexchange coupling in binuclear oxygen-bridged Ni(II) complexes. Chemical Physics 1995, 192, 25 – 35. (154) de Graaf, C.; Sousa, C.; de PR Moreira, I.; Illas, F. Multiconfigurational perturbation theory: An efficient tool to predict magnetic coupling parameters in biradicals, molecular complexes, and ionic insulators. J. Phys. Chem. A 2001, 105, 11371– 11378. (155) de Grad, C.; Broer, R.; Nieuwpoort, W. Electron correlation effects on the dd excitations in NiO. Chem. Phys. 1996, 208, 35–43. (156) De Graaf, C.; Broer, R.; Nieuwpoort, W. Comparison of the superexchange interaction in NiO and in a NiO [100] surface. Chem. Phys. Lett. 1997, 271, 372–376. (157) De Graaf, C.; Illas, F.; Broer, R.; Nieuwpoort, W. On the magnetic coupling in NiO. J. Chem. Phys. 1997, 106, 3287–3291. (158) Evjen, H. On the stability of certain heteropolar crystals. Phys. Rev. 1932, 39, 675. (159) Doll, K.; Dolg, M.; Fulde, P.; Stoll, H. Quantum chemical approach to cohesive properties of NiO. Phys. Rev. B 1997, 55, 10282. (160) Wyckoff, R. W. G. Crystal Structures, 2nd ed.; Interscience Publishers, New York, 1965. (161) Hutchings, M. T.; Samuelsen, E. Measurement of spin-wave dispersion in NiO by inelastic neutron scattering and its relation to magnetic properties. Phys. Rev. B 1972, 6, 3447. (162) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. 58


Page 58 of 61

Page 59 of 61 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


(163) Massey, M.; Chen, N.; Allen, J.; Merlin, R. Pressure dependence of two-magnon Raman scattering in NiO. Phys. Rev. B 1990, 42, 8776. (164) Zadrozny, J. M.; Atanasov, M.; Bryan, A. M.; Lin, C.-Y.; Rekken, B. D.; Power, P. P.; Neese, F.; Long, J. R. Slow magnetization dynamics in a series of twocoordinate iron (II) complexes. Chem. Sci. 2013, 4, 125–138. (165) Atanasov, M.; Zadrozny, J. M.; Long, J. R.; Neese, F. A theoretical analysis of chemical bonding, vibronic coupling, and magnetic anisotropy in linear iron (II) complexes with single-molecule magnet behavior. Chem. Sci. 2013, 4, 139–156. (166) Polly, R.; Werner, H.-J.; Manby, F. R.; Knowles, P. J. Fast Hartree–Fock theory using local density fitting approximations. Mol. Phys. 2004, 102, 2311–2321. (167) Győrffy, W.; Shiozaki, T.; Knizia, G.; Werner, H.-J. Analytical energy gradients for second-order multireference perturbation theory using density fitting. J. Chem. Phys. 2013, 138, 104104. (168) Weigend, F. Hartree–Fock exchange fitting basis sets for H to Rn. J. Comput. Chem. 2008, 29, 167–175. (169) Schweinfurth, D.; Krzystek, J.; Schapiro, I.; Demeshko, S.; Klein, J.; Telser, J.; Ozarowski, A.; Su, C.-Y.; Meyer, F.; Atanasov, M.; Frank, N.; Biprajit, S. Electronic Structures of Octahedral Ni (II) Complexes with “Click” Derived Triazole Ligands: A Combined Structural, Magnetometric, Spectroscopic, and Theoretical Study. Inorg. Chem. 2013, 52, 6880–6892. (170) Weigend, F. A fully direct RI-HF algorithm: Implementation, optimised auxiliary basis sets, demonstration of accuracy and efficiency. Phys. Chem. Chem. Phys. 2002, 4, 4285–4291.

59



(171) Glass, R.; Hibbert, A. Relativistic effects in many electron atoms. Comput. Phys. Commun. 1978, 16, 19–34.

60


Page 60 of 61

Page 61 of 61 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


For Table of Contents Only

61


Embedded Multireference Coupled Cluster Theory - Journal of

Recommend Documents