Spin-Adapted Formulation and Implementation of Density Cumulant

Article pubs.acs.org/JCTC

Spin-Adapted Formulation and Implementation of Density Cumulant Functional Theory with Density-Fitting Approximation: Application to Transition Metal Compounds Xiao Wang,† Alexander Yu. Sokolov,*,‡ Justin M. Turney,† and Henry F. Schaefer*,† †

Center for Computational Quantum Chemistry, University of Georgia, Athens, Georgia 30602, United States Department of Chemistry, Princeton University, Princeton, New Jersey 08544, United States

‡

S Supporting Information *

ABSTRACT: Density cumulant functional theory (DCT) has recently emerged as an attractive ab initio approach for the treatment of electron correlation. In its orbital-optimized formulation (ODC-12) [J. Chem. Phys. 139, 204110 (2013)], DCT has been shown to provide reliable results for a variety of challenging chemical systems. Among the attractive properties of DCT are its size-consistency and size-extensivity, as well as the efficient computation of the molecular properties and analytic gradients. In this work, we present a new formulation and implementation of DCT that takes advantage of spin adaptation and the density-fitting approximation (DF-ODC-12). Our new spin-adapted DF-ODC-12 implementation is more efficient than the previous ODC-12 implementation with up to a ∼12-fold speed-up. We demonstrate the capabilities of DFODC-12 with a study of transition metal compounds, which require high levels of electron correlation treatment. For transition metal carbonyl complexes [Fe(CO)5, Cr(CO)6] and the ferrocene molecule [Fe(Cp)2], the DF-ODC-12 equilibrium parameters and bond dissociation energies extrapolated to the complete basis set limit are in very good agreement with reference data derived from experiment. particle density cumulant.22−25 While DCT is formally related to methods that are based on the variational optimization of the two-particle reduced density matrix,24,26−30 it also has relationships with popular wave function-based approaches. In particular, some of the DCT formulations (such as the orbital-optimized ODC-12 method)19,20 are related to linearized CC theory with double excitations, while the higher-level DCT variants can be considered finite-order approximations to the unitary and variational CC theory.21 Not surprisingly, DCT shares many appealing features with CC theory, in particular, its size-extensivity and size-consistency. In contrast to the traditional CC methods5−7 [such as CCSD or CCSD(T)], DCT is formulated using a Hermitian and stationary energy functional, which makes the computation of molecular properties very efficient.20 DCT methods have been shown to reliably describe molecules with open-shell electronic states, as well as symmetry-breaking and multireference effects,20,21,31,32 making it particularly attractive for computational studies of transition metal compounds. However, the performance of ODC-12 for transition metal compounds has never been tested. Herein, we present an efficient formulation and implementation of the ODC-12 method. Although the ODC-12 computa-

1. INTRODUCTION A reliable theoretical description of transition metal compounds presents many challenges. While the electronic structure of many transition metal systems near geometric equilibrium can be qualitatively described using Hartree−Fock or density functional theory (DFT),1 to obtain quantitatively accurate results, dynamic electron correlation effects need to be captured. Single-reference ab initio methods were developed to accomplish this task, but their application to transition metal compounds in combination with large one-electron basis sets is computationally demanding.2−4 In particular, coupled cluster theory5−7 (CC) has been demonstrated to achieve reliable results for a number of transition metal systems, provided that costly triple (and higher-order) excitations are included in the expansion of the electronic wave function.8−10 Apart from concerns about efficiency, the applicability of single-reference CC theory has been limited due to the complicated evaluation of molecular properties,11−14 as well as sensitivity to multireference and symmetry-breaking effects, which are particularly common in transition metal compounds.15 In this regard, the development of accurate yet efficient theoretical methods for the description of transition metal compounds is still an active area of research. We have recently developed and implemented density cumulant functional theory16−21 (DCT), an electronic structure approach that describes electron correlation using the two© 2016 American Chemical Society

Received: June 8, 2016 Published: September 8, 2016 4833

DOI: 10.1021/acs.jctc.6b00589 J. Chem. Theory Comput. 2016, 12, 4833−4842

Article

Journal of Chemical Theory and Computation via the following expression:22−25

tional cost is similar to that of the linearized coupled cluster theory with double excitations, our original implementation was limited to modest systems (300−400 molecular orbitals) due to the additional costs associated with the transformation and storage of the two-electron integral tensors at every iteration. In this work, we demonstrate that the efficiency of the DCT algorithm can be greatly improved using the density-fitting approximation for two-electron integrals,33−38 which significantly reduces the cost of the integral transformation. In addition, we present a spin-adapted formulation of DCT that achieves further computational savings for closed-shell molecules. To demonstrate the capabilities of our new ODC12 implementation, we present computations of the equilibrium structures and bond dissociation energies for a number of transition metal complexes, including the computation of the ferrocene complex with a quadruple-ζ basis set (970 basis functions). We begin our discussion with a brief overview of DCT (section 2.1) and its implementation for the ODC-12 method (section 2.2). We then describe how the density-fitting approximation can be used to reduce the computational cost of the ODC-12 implementation (section 2.3). In section 2.4, we present the spin-adapted formulation of DCT for closedshell molecules. We evaluate the performance of our new implementation in section 3 and present computational details in section 4. Finally, we present the ODC-12 results for equilibrium properties of transition metal compounds in section 5 and conclude our discussion in section 6.

λrspq = γrspq − γrpγsq + γrqγsp

Equations 6 and 7 can be used to compute γpq, λpq rs , and the electronic energy [eq 4] directly from a many-electron wave function |Ψ⟩. In DCT we take an alternative approach, which avoids the costly computation of |Ψ⟩ by minimizing the energy functional [eq 4] with respect to the density cumulant. To accomplish this, the one-particle reduced density matrix is exactly decomposed into two contributions:16

γpq = κpq + τpq κqp

⎛ ∂E ⎞ =0 ⎟ ⎜ ⎝ ∂λ2 ⎠ N − rep

(1)

⎛ ∂E ⎞ ⎜ ⎟ ⎝ ∂X ⎠

where and ap are the Fermion creation and annihilation operators, respectively, while hqp and vrspq are the usual one- and antisymmetrized two-electron integrals:

rs g pq = ⟨ψp(1)ψq(2)|

1 |ψ (1)ψs(2)⟩ r12 r

where

f qp

(3)

(4)

E=

are the elements of the generalized Fock operator

qs r f pq = hpq + vpr γs

(5)

The two-particle density cumulant is related to the familiar oneand two-particle reduced density matrices γqp = ⟨Ψ|ap†aq|Ψ⟩, γrspq = ⟨Ψ|ap†aq†asar |Ψ⟩

=0 (10)

X=0

where elements of X parametrize the unitary orbital trans† formation U = eX−X . A similar parametrization is used in the wave function-based orbital-optimized approaches, such as the orbital-optimized coupled-cluster and perturbation theories.44−47 In the current work, we will restrict our discussion to the ODC-12 variant of DCT,19,20 which formulates the stationarity conditions [eqs 9 and 10] using the Nrepresentability constraints derived from second-order perturbation theory.48 2.2. Implementation of the ODC-12 Method. In this section, we briefly describe the implementation of the ODC-12 method. Using approximate N-representability conditions, the ODC-12 energy can be written in the following form:16,17,20

(2)

Indices p, q, r, s run over the entire spin−orbital basis ψp. In eq 1, we use the Einstein convention for the summation over repeated indices, and we will continue to do so henceforth. In DCT, we express electronic energy as a known functional of the one-particle reduced density matrix (γqp) and the twoparticle density cumulant (λrspq):16 1 1 rs pq E = ⟨Ψ|Ĥ |Ψ⟩ = (hpq + f pq )γqp + vpq λrs 2 4

(9)

where the derivative is constrained to variations of the density cumulant (λ2), which ensure that the resulting density matrices are N-representable, that is, can be derived from a physical Nelectron wave function.40−43 An important property of DCT is that the errors in the cumulant N-representability only affect the small cumulant-dependent (correlation) part of the electronic energy, which allows for approximate N-representability constraints to be employed. Although the energy [eq 4] is invariant to orbital transformations if exact N-representability constraints are used, for an approximate set of N-representability conditions the energy will depend on the choice of orbitals. The optimal set of orbitals is chosen by minimizing the energy functional with respect to all nonredundant orbital rotations20

a†p

rs rs sr hpq = ⟨ψp(1)|h|̂ ψq(1)⟩, vpq = g pq − g pq

(8)

τqp

where and are the mean-field and correlation components of γqp, respectively. Conveniently, in eq 8 only τqp depends on the density cumulant, while κqp is independent of λpq rs and can be obtained by specifying a set of molecular orbitals.20 The electronic energy can thus be obtained by minimizing the energy functional with respect to cumulant parameters:

2. THEORY 2.1. Overview of Density Cumulant Functional Theory. We begin our discussion with a brief overview of density cumulant functional theory (DCT).16−21 First, we define the electronic Hamiltonian in the form39 1 rs † † Ĥ = hpqap†aq + vpq ap aq asar 4

(7)

1 q 1 (hp + f pq )γqp + vijabλabij 2 2 1 ab kl ij 1 cd ab ij + λkl vij λab + λij vcd λab − λikacvjckbλabij 8 8

(11)

(κii

where indices i, j, k, l and a, b, c, d denote the occupied = 1) and virtual (κaa = 0) spin−orbitals, respectively. In eq 11, we have expressed the ODC-12 energy as a function of the occupied-virtual density cumulant elements λijab. The remaining

(6) 4834


Article

Journal of Chemical Theory and Computation ij elements of λpq rs do not appear explicitly and are related to λab through the N-representability conditions. To evaluate the energy, we ensure that it [eq 11] is stationary with respect to variations of the λijab parameters [eq 9] and orbital rotations [eq 10]:19,20

R ijab =

(d) Evaluate the orbital residual Rai [eq 13]. Check convergence of the energy, orbitals, and cumulant. If ΔE = |Enew − Eold| and the root-mean-square of Rai and Rab ij are less than a predefined threshold (ϵconv), then the ODC-12 procedure is complete. Otherwise, proceed to the next step. (e) Update the molecular orbitals:

∂E c i ij = P −(ab)[Fã λcbij] − P −(ij)[Fk̃ λabkj] + vab ∂λijab

Xia =

1 1 cd ij + vklijλabkl + vab λcd − P −(ij)P −(ab)[vkbjcλacik ] = 0 2 2

(12)

denotes a modified Fock matrix [see eq 16 below]

P −(ij)f (i , j) ≡ f (i , j) − f (j , i)

(14)

is the antisymmetric permutation operator. The two sets of coupled residual equations Rai and Rab ij are solved iteratively until self-consistency is achieved. The ODC-12 computational algorithm can be summarized as follows: 1. Obtain an initial guess for the cumulant parameters λijab and molecular orbitals. We typically begin the ODC-12 computation using the converged Hartree−Fock orbitals; however, other types of orbitals can be used [e.g., Kohn− Sham orbitals or second-order Møller−Plesset (MP2) natural orbitals]. Similarly, for the density cumulant, MP2 amplitudes are usually chosen as the guess (λijab = tijab). We have observed that the converged ODC-12 energy is rather insensitive to initial parameters as long as the ground-state wave function is qualitatively well represented by a single Slater determinant. 2. Begin iterations: (a) From the density cumulant, compute the correlation contribution τqp to γqp 1 1 τi j = − λikabλabjk − τikτkj , τab = λacijλijbc + τacτcb 2 2

rs g pq ≈

Y bpq =

j

1 + τi ′ + τj ′

b

Fã =

,

=

+

(pq|X ) = JXY =

1 − τa ′ − τb ′

a

∬ ψp(1)ψq(1) r1 χX (2) dr1 dr2

∬ χX (1) r1 χY (2) dr1 dr2 12

(21)

(22)

In the above equations, the capital indices X, Y label the auxiliary atom-centered basis set χX fitted to approximately satisfy eq 19 in the least-squares sense.33,34 This decomposition reduces the storage of the two-electron integrals from 6(N 4) to 6(NauxN 2) (Naux is the number of auxiliary basis functions), which can greatly reduce the disk usage and the input/output (I/O) operation count, while maintaining high accuracy of the approximation. For this reason, density-fitting has become a widely used tool in quantum chemistry. 51−57 In our implementation of DCT with density-fitting (DF-DCT), the antisymmetrized two-electron integrals vrspq in eqs 11 to 13 are computed from the bYpq coefficients directly in memory for each symmetry sector, which significantly reduces the I/O cost of 1 cd ij the ladder ( 2 vab λcd ) tensor contraction. Additional computational savings are achieved for the two-electron integral transformation, which is performed in each DCT iteration.

f ab′′ Y aa ′(Y †)bb ′

j

(20)

12

R ijab i

∑ (pq|X )[J−1/2 ]XY

where the coefficients are defined in terms of the two- and three-index integral tensors:

(16)

(λijab)old

Naux

bYpq

where Yii′ and Yaa′ are the eigenvector matrices obtained by diagonalizing τji and τba, respectively, with the corresponding eigenvalues τi′ and τa′.19 (b) Evaluate the cumulant residual Rab ij [eq 12] and update the density cumulant elements: (λijab)new

(19)

X

and update the generalized Fock matrix according to eq 5. The modified Fock elements can be computed as Fĩ =

(18)

Naux

∑ bprX bqsX X

(15)

fi ′j ′ Y ii ′(Y †) jj ′

†

U = e X−X ,

,

(f) Transform the integrals using the new set of coefficients C(new) and return to the start of the loop. We employ the direct inversion in the iterative subspace (DIIS) method49,50 to accelerate the iterative solution of the residual equations. The ODC-12 algorithm outlined above has 6(O2V 4) computational scaling (“O” and “V” are the number of occupied and virtual spin−orbitals, respectively), which 1 cd ij originates due to the “ladder” ( 2 vab λcd ) term in the cumulant ab residual Rij . The solution of the equations for the orbitals has a lower 6(O3V 3) scaling, but requires the costly 6(N 5) full integral transformation (N = O + V). As we will describe in the next section, these computational costs can be significantly reduced by employing the density-fitting approximation. 2.3. DCT with Density-Fitting. In the conventional density-fitting approach (DF),33−38 the two-electron integrals are approximated in the following form:

(13)

where and

a − Fã

C(new) = C(old)U

1 cd jk ab 1 ∂E b a j i vib λcd λjk − vajklλklbcλbcij = a = f i γb − fa γj + 4 4 ∂Xi 1 jk ab 1 bc ij + vib λjk − vaj λbc − vjikbλbcjlλklac + vabjcλcdik λjkbd = 0 2 2

R ia

F̃qp

R ia i Fĩ

b

Fĩ + Fj̃ − Fã − Fb̃

(17)

(c) From eq 11 determine the ODC-12 energy. 4835


Article

Journal of Chemical Theory and Computation

Figure 1. Total wall time of single-point energy computations using the UODC-12, RODC-12, DF-UODC-12, and DF-RODC-12 methods for two sets of (a) small- and (b) medium-sized systems selected from the S22 database59 (cc-pVTZ basis set). The wall time contributions from individual steps of the ODC-12 algorithm are shown as colored boxes (see section 2.2 for details). The “orbital update” refers to the orbital residual eq 13, the 1 cd ij λcd term in the cumulant residual eq 12, and the “cumulant update” refers to the remaining terms of the “V4O2 ladder term” corresponds to the 2 vab cumulant residual. All computations were performed using just one processor (Intel Xeon E5−2637 CPU @ 3.50 GHz).

spin (σ = α, β) components. The DCT equations can be significantly simplified if we write the electronic Hamiltonian in the spin-summed (spin-adapted) form:

Using the density-fitting approximation, the cost of the full integral transformation is reduced dramatically [from 6(N 5) to 6(NauxN3)], since only the transformation of the three-index tensors bYpq is required. As in conventional correlated and orbital-optimized densityfitted methods, our DF-DCT implementation uses two types of auxiliary basis sets. For the mean-field (κqp-dependent) terms of eqs 11 to 13, the two-electron integrals (gklij ) are evaluated using JKFIT basis sets that are optimized for fitting Coulomb and exchange integrals [the occupied−occupied ψi(1)ψj(1) orbital products].38 For the remaining (cumulant-dependent) terms, the grspq integrals are fitted using RI auxiliary basis sets.58 As a result, in the limit of a zero density cumulant (λrspq = 0), the DFDCT energy yields exactly the energy of the density-fitted selfconsistent field method (DF-SCF).38 2.4. Spin-Adapted Formulation of DCT for ClosedShell Systems. In sections 2.1 and 2.2, we presented a formulation of DCT in a general spin−orbital basis, where each spin−orbital ψp is expressed as a product of the spatial (ϕP) and

1 RS PQ Ĥ = hPQ EQP + gPQ ERS 2 EQP =

∑ aP†σ aQσ , ERSPQ

(23)

=

σ

∑ aP†σ aQ†ρaSρaRσ (24)

σ ,ρ

where capital indices P, Q, R, ... denote spatial orbitals. Using eq 23, we can express the DCT electronic energy in terms of the spin-free reduced density matrices ΓPQ and ΓPQ RS E = hPQ ΓQP +

ΓPQ

1 RS PQ g Γ RS 2 PQ

⟨Ψ|EPQ|Ψ⟩

(25)

ΓPQ RS

⟨Ψ|EPQ RS |Ψ⟩.

where = and = We are now in a position to define the spin-free density 25 cumulant ΛPQ RS . As discussed by Kutzelnigg, for a general wave PQ function |Ψ⟩, the definition of ΛRS depends on the spinquantum numbers S and MS. If we assume that |Ψ⟩ is a closed4836


Article

Journal of Chemical Theory and Computation shell electronic wave function (S = 0), ΛPQ RS can be written in a simple form: PQ P Q ΛPQ RS = Γ RS − Γ RΓ S +

1 Q P Γ R ΓS 2

Table 1. Nonbonded Interaction and Radical Stabilization Energies Computed Using the ODC-12 and DF-ODC-12 Methods (kcal mol−1)a

(26)

ODC-12

PQ Unlike the spin−orbital λpq rs , the spin-free ΛRS is no longer antisymmetric with respect to the exchange of two indices in the bra or ket. Nevertheless, ΛPQ RS still satisfies several important properties of λrspq. In particular, its partial trace contains information about the correlation contribution to ΓPQ:

⎛1 2 ⎞Q Q ⎜ Γ − Γ⎟ 1 1⎠ ≡ DP 2 P

∑ ΛQR PR = ⎝ R

H2O···NH3 HCN···HCN HF···HF NH3···NH3 CH2O···CH2O C2H4···C2H4 C2H4···Ar C2H2···C2H2

(27)

1

where DQP is the deviation of 2 Γ1 from idempotency. Equation 27 suggests that it is possible to decompose ΓPQ as follows: ΓQP = KPQ + TPQ

·CH2NO2 ·CH2CN ·CH2NH+3 ·CH2OH ·CH2CCH ·CH2BH2 ·CHO ·CHCH2

(28)

1 K 2

where is the idempotent contribution. The correlation component TQP can be obtained from ΛRS PQ according to eq 27. These spin-adapted equations can be used to formulate the ODC-12 energy and the N-representability and stationarity conditions. We refer readers to the Appendix (section 7) for the explicit form of these equations.

DF-ODC-12 error

Nonbonded Interaction Energiesb −7.0 −7.0 −5.1 −3.3 −4.8 −1.7 −0.6 1.0 Radical Stabilization Energiesc 42.1 −8.7 4.6 −8.7 −13.2 73.1 −15.5 5.9

1.2 5.6 2.1 3.0 1.5 1.2 3.0 3.0

× × × × × × × ×

10−3 10−3 10−3 10−4 10−3 10−3 10−4 10−4

2.9 3.0 2.0 3.1 1.7 1.0 4.0 8.8

× × × × × × × ×

10−3 10−3 10−3 10−3 10−3 10−2 10−3 10−3

a

For DF-ODC-12, only the unsigned error relative to the ODC-12 energy is shown. The nonbonded interaction energies were computed using the aug-cc-pVTZ basis set, while for the radical stabilization energies the cc-pVTZ basis set was used. The corresponding JKFIT and RI auxiliary basis sets were used for the DF-ODC-12 method. Molecular geometries were taken from refs 60 and 61. bThe nonbonded interaction energy is defined as the enthalpy of the association reaction: X + Y → X ··· Y. cThe radical stabilization energy is defined as the enthalpy of the hydrogen-transfer reaction: ·CH3 + RH → CH4 + ·R

3. ANALYSIS OF COMPUTATIONAL EFFICIENCY AND DENSITY-FITTING ERRORS In this section we analyze the efficiency of our new ODC-12 implementation. We refer to the original (nonspin-adapted) ODC-12 implementation as the unrestricted ODC-12 method (UODC-12) and the spin-adapted (restricted) implementation as RODC-12. Figure 1 shows the relative computational efficiency of UODC-12 and RODC-12, as well as their formulations with the density-fitting approximation (DFUODC-12 and DF-RODC-12), by comparing the computational wall time for a number of closed-shell systems selected from the S22 database.59 For all systems considered, spin-adaptation of the ODC-12 method results in about a 3-fold reduction in computational time, and using the density-fitting approximation gives rise to additional computational savings. Figure 1 demonstrates that density-fitting greatly reduces the time necessary to perform the two-electron integral transformation, as well as the evaluation of the 6(O2V 4) ladder tensor contraction (see section 2.3 for details). As a result, both DF-UODC-12 and DF-RODC-12 are about 2−5× more efficient than the conventional counterparts for the molecular systems under consideration. Combining the effect of density-fitting and spin-adaptation in the DF-RODC12 method results in an overall 6- to 13-fold reduction in computational cost for closed-shell systems, relative to our original UODC-12 implementation. We now turn our attention to the analysis of errors introduced by the density-fitting approximation in DF-ODC12. Table 1 shows the errors of DF-ODC-12 relative to the conventional ODC-12 method for the interaction energies of closed-shell nonbonded complexes and radical stabilization energies selected from the A24 and RSE42 benchmark databases.60,61 In all systems the effect of the density-fitting approximation on the computed relative energies is negligibly small, with the largest errors of ∼10−2 kcal mol−1. These results are consistent with the previous analyses of the density-fitting

approximation, where errors of a similar magnitude have been found.51,52,56,58,62−65

4. COMPUTATIONAL DETAILS The ODC-12 method with spin-adaptation (denoted as RODC-12) and density-fitting (DF-RODC-12 and DFUODC-12) have been implemented in the developer’s version of the PSI4 program.66 The equilibrium structures of all transition metal compounds were fully optimized using our ODC-12 analytic gradient code in PSI4. All electrons were correlated in the ODC-12 computations. We employed the def2-XZVPP (X = T, Q) basis sets of Weigend and Ahlrichs,67 unless noted otherwise. The corresponding auxiliary def2XZVPP-JKFIT and def2-XZVPP-RI basis sets67−69 were used for the DF-ODC-12 method, as described above in section 2.3. The dissociation energies of all compounds were computed as the difference between the total energy of the dissociated products and the energy of the reactants. To compute the dissociation energies of the 3d-transition metal complexes (section 5.2), the molecular energies were extrapolated to the complete basis set (CBS) limit using the following equations: Eref = A + Be−α

X

Ecorr = A + BX −β

(29) (30)

where Eref is the SCF (DF-SCF) energy and Ecorr is the ODC12 (DF-ODC-12) correlation energy (Ecorr = EODC-12 − Eref) computed using the def2-XZVPP basis sets (X = T, Q). The optimal values of the exponential parameters α = 7.88 and β = 4837


Article


Table 2. Equilibrium Bond Lengths (re, Å) and Dissociation Energies (De, kcal mol−1) Obtained Using the CCSD and DF-ODC12 Methods (def2-QZVPP basis set)a re

De

molecule

state

DF-ODC-12

CCSD

ref

DF-ODC-12

CCSD

ref

MnCl ZnCl FeCl CrCl ZnS ZnH CuCl

7 +

Σ 2 Σ 6 Δ 6 + Σ 1 Σ 2 + Σ 1 Σ

2.233 2.145 2.178 2.198 2.070 1.582 2.076

2.232 2.142 2.176 2.193 2.056 1.578 2.076

2.243 2.100 2.179 2.194 2.100 1.590 2.050

81.6 51.9 81.4 89.2 31.8 24.3 86.6

82.6 52.4 82.5 87.0 29.0 25.2 85.7

81.5 54.3 79.4 90.9 34.9 21.6 88.5

ZnO NiCl TiCl CuH VO VCl MnS CrO CoCl VH FeH CrH

1

Σ Π 4 Φ 1 + Σ 4 − Σ 5 Δ 6 + Σ 5 Π 3 Φ 5 Δ 4 Δ 6 + Σ

1.730 2.083 2.264 1.468 1.574 2.225 2.078 1.660 2.092 1.683 1.568 1.651

1.709 2.083 2.265 1.471 1.573 2.227 2.072 1.601 2.104 1.682 1.565 1.643

1.800 2.073 2.265 1.463 1.589 2.215 2.070 1.621 2.087 1.730 1.568b 1.656

30.8 88.1 99.5 63.5 132.8 98.9 61.6 91.0 89.5 59.7 44.1 53.4

26.1 87.4 97.0 62.6 133.9 99.4 59.3 84.7 87.6 60.3 44.7 51.4

38.1 90.1 101.7 62.6 152.1 103.2 71.1 104.5 81.5 51.8 44.8b 46.8

0.018

0.020

5.0

6.0

MAE

2

a

For the DF-ODC-12 method, the corresponding JKFIT and RI auxiliary basis sets were used. The two sets of molecules correspond to the 3dMLBE-SR7 and 3dMLBE-MR13 datasets of ref 4. The reference re and De values for all systems but FeH were taken from ref 4. The reference De values do not include the spin-orbit coupling corrections reported in ref 4. bReference values are taken from ref 8.

for ΔDe. In particular, the ΔDe values of both methods are substantially smaller for 3dMLBE-SR7, relative to 3dMLBEMR13, which may be explained by the higher multireference character in the latter. For the 3dMLBE-SR7 data set, DFODC-12 exhibits a small ΔDe value of 2.0 kcal mol−1, while CCSD has the larger ΔDe = 3.2 kcal mol−1. For the 3dMLBEMR13 set, ΔDe values of both methods increase by more than a factor of 2, with ΔDe = 6.7 and 7.7 kcal mol−1 for DF-ODC-12 and CCSD, respectively. Out of 12 diatomics in the 3dMLBEMR13 data set, both methods exhibit large ΔDe errors of ≥5 kcal mol−1 for seven molecules (ZnO, VO, MnS, CrO, CoCl, VH, and CrH). In these challenging systems, DF-ODC-12 predicts smaller ΔDe errors in five systems (ZnO, VO, MnS, CrO, VH), with a substantial 6.4 kcal mol−1 improvement over CCSD for CrO. 5.2. Complexes of 3d-Transition Metals. Finally, we assess the accuracy of the ODC-12 method for predicting the equilibrium properties of 3d-transition metal complexes. Table 3 presents the optimized structural parameters and dissociation energies (De) for the iron carbonyl, chromium carbonyl, and ferrocene complexes [Fe(CO)5, Cr(CO)6, and Fe(Cp)2], computed using CCSD(T), DF-ODC-12, as well as the popular density functional method B3LYP. These complexes have been the subject of many computational studies,2,10,71−85 partly due to their high point-group symmetry, as well as the availability of the gas-phase experimental results. For Fe(CO)5 and Cr(CO)6, the molecular geometries were fully optimized using the def2-TZVPP basis set. For Fe(Cp)2, the B3LYP and DF-ODC-12 geometries (Figure 2) were optimized using the def2-TZVP basis, while the CCSD(T) optimized geometry was taken from ref 82. For all three methods, the equilibrium bond

2.97 have been suggested by Neese and Valeev based on the extrapolation of the CCSD(T) energies to the CBS limit.70 We have verified that the same parameters yield accurate CBSextrapolated dissociation energies for the ODC-12 method with a benchmark set of diatomic molecules and have adopted these parameters for the ODC-12 basis set extrapolation using eqs 29 and 30 (see Supporting Information for details).

5. RESULTS AND DISCUSSION 5.1. Transition Metal Diatomics. In this section, we investigate the performance of ODC-12 for diatomic molecules containing transition metals. Table 2 presents the equilibrium bond lengths (re) and dissociation energies (De) obtained using the CCSD and DF-ODC-12 methods for a set of 19 diatomics, which comprise the 3dMLBE20 benchmark database reported by Truhlar and co-workers (we exclude CoH from the original data set due to the SCF convergence problems in CCSD).4 As suggested by ref 4, the 19 molecules are divided into two subsets: the 3dMLBE-SR7 set with predominantly singlereference character of the molecular wave functions and the 3dMLBE-MR13 set with significant multireference character. We used the re and De values from ref 4 as our benchmark for all systems except for FeH for which reference values are taken from recent highly accurate ab initio computations.8 The mean absolute errors (MAE) in computed re and De (Δre and ΔDe) averaged over all 19 molecules are shown in Table 2. Both CCSD and DF-ODC-12 exhibit similar MAEs, with the DF-ODC-12 errors smaller than CCSD by 0.002 Å and 1.0 kcal mol−1 for re and De, respectively. When comparing MAEs of the two methods for the 3dMLBE-SR7 and 3dMLBEMR13 subsets individually, noticeable differences are observed 4838


Article


43.7 kcal mol−1, Table 3) yields somewhat better agreement with experiment (De = 40.7 kcal mol−1) than DF-ODC-12 with errors of 3 kcal mol−1. DF-ODC-12 (De = 35.1 kcal mol−1) underestimates the available experimental values by 5.6 kcal mol−1. The B3LYP method (De = 28.9 kcal mol−1) significantly underestimates either experimental dissociation energy by more than 10 kcal mol−1. The Cr(CO)6 complex has an octahedral Oh structure with equivalent Cr−C and C−O bonds (Figure 2b). The obtained CCSD(T) and DF-ODC-12 distances of the Cr−C bonds are in a very good agreement with experiment, while B3LYP significantly overestimates the Cr−C bond length (by ∼0.01 Å). For the carbonyl ligand bond dissociation energies [Cr(CO)6 → Cr(CO)5 + CO], the DF-ODC-12 method (De = 40.9 kcal mol−1) shows the best agreement with the available experimental data (De = 37.9, 41.3 kcal mol−1). The CCSD(T) method overestimates the experimental result by ∼4 kcal mol−1, while the density functional theory B3LYP further underestimates the observed dissociation energy (De = 35.7 kcal mol−1). The ferrocene complex’s [Fe(Cp)2] gas-phase equilibrium structure has D5h symmetry (eclipsed conformation) as shown in Figure 2c. Although the Fe(Cp)2 ground-state wave function has predominantly single-reference character, the accurate determination of its equilibrium geometry and heterolytic dissociation energy [Fe(Cp)2 → Fe2+ + 2Cp−] has been very challenging for theoretical methods.10,81−84 The computational difficulty also arises from the relatively large size of the system (970 basis functions for ferrocene with the def2-QZVPP basis set). Out of the three methods compared in Table 3, CCSD(T) and DF-ODC-12 exhibit the best agreement with experiment for the equilibrium structural parameters, while B3LYP noticeably overestimates the Fe−C bond length (by 0.02− 0.03 Å). For the heterolytic ligand dissociation energies, the results of the three methods are different by a larger extent. In particular, DF-ODC-12 exhibits an excellent agreement with experiment, with an error of just 1 kcal mol−1. Both CCSD(T) and B3LYP errors are much larger: the former method overestimates the experimental De value by 15 kcal mol−1, while the latter method underestimates it by 16 kcal mol−1. We stress, however, that the CBS-extrapolated De energies of all three methods may be affected by the lack of the additional core−valence correlation basis functions in the def2-XZVPP basis sets. 97 As we demonstrated in the Supporting Information, the inclusion of the core−valence correlation basis functions in transition metal diatomics can amount up to 2 kcal mol−1 difference in De at the CBS limit, which is significantly smaller than the relative difference in De between the three methods, as well as the experimental uncertainty. To account for these effects, the development of density-fitting basis sets with additional core−valence correlation functions is highly desirable.

Table 3. Equilibrium Structural Parameters (re, Å) and Dissociation Energies (De, kcal mol−1) Obtained Using the B3LYP, CCSD(T), and DF-ODC-12 Methodsa B3LYP

CCSD(T)

DF-ODC-12

r(Fe−C)ax

1.827

Fe(CO)5 1.809e 1.819

r(Fe−C)eq

1.821

1.812e

1.812

r(C−O)ax

1.137

1.145

e

1.134

r(C−O)eq

1.140

1.147e

1.138

Deb

28.9

43.9e

r(Cr−C) r(C−O) Dec

1.924 1.139 35.7

r(Fe−C)

2.079

35.1 Cr(CO)6 1.913e 1.916 1.137 1.146e 44.0e 40.9 Fe(Cp)2 2.047f 2.061

r(C−C)

1.423

1.427f

1.423

r(C−H)

1.078

1.079f

1.079

∠C5−H Ded

0.87 620

0.52f 651e,f

1.43 637

expt

expt ref

1.811, 1.807, 1.810 1.803, 1.827, 1.842 1.117, 1.152, 1.142 1.133, 1.152, 1.149 40.7

86, 87, 101 86, 87, 101 86, 87, 101 86, 87, 101 88

1.918, 1.916 1,141, 1.140 37.9, 41.3

91, 93 91, 93 89, 92

2.054, 2.031, 2.064 1.435, 1.415, 1.440 1.080, 1.073, 1.104 3.7, 1.7 636

10, 90, 94 10, 90, 94 10, 90, 94 90, 95 96

a

The Fe(CO)5 and Cr(CO)6 geometries were optimized using the def2-TZVPP basis set, while the def2-TZVP basis was used for Fe(Cp)2. The computed De values were extrapolated to the complete basis set limit (see section 4 for details). The experimental De values were corrected by the zero-point vibrational, thermal, and relativistic contributions (see Supporting Information for details). bDe is the enthalpy of the reaction: Fe(CO)5 (1A′1) → Fe(CO)4 (3B2) + CO (1Σ+). cDe is the enthalpy of the reaction: Cr(CO)6 (1A1 g) → Cr(CO)5 (1A1) + CO (1Σ+). dDe is the enthalpy of the reaction: Fe(Cp)2 (1A1′) → Fe2+ (5D) + 2Cp− (1A1′ ). eThe details of the CCSD(T) computations are reported in the Supporting Information. f The CCSD(T) optimized geometry for Fe(Cp)2 was taken from ref 82 and was used to compute De at the CCSD(T)/CBS level of theory.

dissociation energies (De) were computed at the optimized geometries using the def2-XZVPP (X = T, Q) basis sets and then extrapolated to the CBS limit, as described in section 4. The reference De energies presented in Table 3 were derived from experimental dissociation energies86−96 by subtracting out the computed zero-point vibrational, thermal, and relativistic energy contributions, in order to directly compare with the results of ab initio computations. Details regarding this procedure can be found in the Supporting Information. The D3h gas-phase structure of Fe(CO)5 (Figure 2a) has been the subject of many computational and experimental studies.2,71−78,86−89 Despite this, the Fe(CO)5 experimental structural parameters reported in the literature have a large degree of uncertainty, while the results of theoretical computations strongly depend on the method employed. As one can see from Table 3, different experimental studies do not agree about the relative lengths of the axial and equatorial Fe− C bonds. Out of the three methods employed in our study, B3LYP and DF-ODC-12 predict the shorter equatorial Fe−C bond, while CCSD(T) predicts the axial bond to be slightly shorter. Overall, all three methods predict structural parameters that are within the span of the experimental conclusions. When the dissociation energies for one carbonyl ligand are compared, [Fe(CO)5 → Fe(CO)4 + CO], the CCSD(T) method (De =

6. CONCLUSIONS We presented an efficient spin-adapted implementation of density cumulant functional theory (DCT) in its ODC-12 formulation. Taking advantage of the density-fitting approximation for the two-electron integrals, our new implementation (DF-ODC-12) achieves 6- to 13-fold reduction in computational time, relative to the previous implementation. As we demonstrated, these computational savings mainly originate from the lower disk usage and input/output operation count in 4839


Article


Figure 2. Equilibrium structures of the (a) Fe(CO)5, (b) Cr(CO)6, and (c) ferrocene complexes computed using the CCSD(T) and DF-ODC-12 methods, as well as those obtained in experiments.10,86,87,90,91,93,94,101 Only the smallest and the largest experimental values are shown. The Fe(CO)5 and Cr(CO)6 geometries were optimized using the def2-TZVPP basis set, while the def2-TZVP basis was used for Fe(Cp)2.

■

APPENDIX: SPIN-ADAPTED ODC-12 ENERGY AND RESIDUAL EQUATIONS The spin-free second-order cumulant N-representability conditions take the following form:

the two-electron integral transformation and the density cumulant residual equations. To demonstrate the performance of our new DF-ODC-12 implementation, we carried out a benchmark study of the ODC-12 method for a number of transition metal compounds, such as the diatomic molecules containing transition metals, as well as the iron carbonyl [Fe(CO)5], chromium carbonyl [Cr(CO)6], and ferrocene [Fe(Cp)2] complexes. For diatomic molecules, the performance of DF-ODC-12 was found to be similar to that of CCSD, with some improvement over CCSD in computed equilibrium structural parameters and bond dissociation energies. For transition metal complexes, the DFODC-12 optimized geometries and ligand dissociation energies were shown to be in very good agreement with the available experimental results. In particular, for the Cr(CO)6 and Fe(Cp)2 complexes, DF-ODC-12 outperforms CCSD(T), predicting the ligand dissociation energies that are much closer to the reference data derived from experiment. Our results demonstrate that the DF-ODC-12 method can be used to study the equilibrium properties of the small and medium-sized transition metal systems with large basis sets in a reliable and efficient way. Further efficiency improvements can be achieved by incorporating the local approximations57,98−100 and taking advantage of distributed-memory parallelization techniques. An attractive feature of the DF-ODC-12 method is the availability of the one- and two-particle reduced density matrices at the end of the energy computation, which can be efficiently used to compute ground-state equilibrium properties of transition metal complexes, such as optimized geometries and vibrational spectra. In this regard, the generalization of DCT to allow for a description of systems in electronically excited states is of great interest, and will be explored in the future.

ΛIJKL =

1 IJ AB 1 IJ AB ΛABΛKL + ΛABΛLK 3 6

(31)

AB ΛCD =

1 AB IJ 1 IJ ΛIJ ΛCD + ΛAB JI ΛCD 3 6

(32)

1 1 AC IK 1 AC IK 1 AC IK AI IK ΛBJ = − ΛAC ΛJK ΛCB − ΛKJ ΛBC − ΛKJ ΛCB JK ΛBC − 3 6 6 3 (33)

ΛAI JB =

1 AC IK ΛJK ΛBC 2

(34)

where indices I, J, K, L and A, B, C, D label the occupied (KII = 2) and virtual (KAA = 0) spatial orbitals. Using eqs 31 to 34, the spin-adapted ODC-12 energy functional can be written as 1 P 1 IJ AB KL IJ AB (hQ + fQP )ΓQP + gAB ΛIJ + ΛABΛKLgIJ̃ 2 12 1 AB IJ CD 1 BJ ̃ AC IK AC (Λ KJ ΛCB + Λ̃ JK ΛIK + ΛIJ ΛCDgAB ̃ − gAI BC) 12 6 1 JB AC IK + gAI ΛJK ΛBC , (35) 2 RS where we defined the Λ̃RS PQ and g̃PQ intermediates and the generalized Fock matrix f QP : E=

4840

⎛ QS 1 QS ⎞⎟ S f PQ = hPQ + ⎜gPR − gRP ΓR ⎝ 2 ⎠

(36)

RS RS Λ̃ PQ = 2ΛRS PQ + ΛQP

(37) DOI: 10.1021/acs.jctc.6b00589 J. Chem. Theory Comput. 2016, 12, 4833−4842

Journal of Chemical Theory and Computation RS RS gPQ + gQP ̃ RS = 2gPQ

■

ACKNOWLEDGMENTS This research was supported by the U.S. National Science Foundation, Grant No. CHE-1361178. A.Y.S. thanks Professor Werner Kutzelnigg and Dr. Andrew Simmonett for insightful discussions.

(38)

Minimizing the energy functional [eq 35] with respect to the variation of ΛAB IJ and orbital rotations, we obtain the spinadapted cumulant and orbital residual equations: 1 1 IJ KL C IJ I KJ IJ P+(IJ , AB)[FÃ Λ̃CB − FK̃ Λ̃ AB] + 2gAB + gKL ̃ ΛAB 3 3 1 CD IJ 1 IC KJ CJ ̃ KI CI ̃ KJ + gAB ΛCB − gAK Λ BC − gAK ΛCB] ̃ ΛCD + P+(IJ , AB)[3gAK 3 3

RIJAB =

1 CD JK ̃ AB 1 KL ̃ BC IJ g ΛCDΛ JK − gAJ Λ KLΛBC 6 IB 6 1 KB JL AC 1 CJ IK BD JK AB BC IJ + gIB ΛJK − gAJ ΛBC + gIJ ΛBCΛKL − gAB ΛCDΛJK 2 2 1 1 JC ̃ IK BD JL AC JL AC IK − gIJBK (Λ̃ BCΛKL + Λ̃CBΛLK ) + gAB (ΛCDΛJK + Λ̃ DCΛBD KJ ) 6 6

RIA = f IB Γ BAf IA Γ IJ +

■

(39)

(41) (42) (43)

f IJ′′ YII ′(Y †) JJ ′

J

2 + TI ′ + TJ ′

(44)

f AB′′ YAA ′(Y †)BB ′

B

FÃ =

2 − TA ′ − TB ′

YII′

(45)

(YAA′)

where and TI′ (TA′) are the eigenvectors and the corresponding eigenvalues of the TJI (TBA) matrix, which can be evaluated by solving the nonlinear equations:

TIJ = DIJ −

1 K J TI TK 2

TAB = −DAB +

(46)

1 C B T A TC 2 DJI

(47)

DBA

In eqs 46 and 47, and are the occupied- and virtualorbital form of the partial trace of Λ2, respectively, and can be explicitly written as

■

1 ̃ AB DIJ = − ΛIK ABΛ JK 6

(48)

1 BC DAB = − ΛIJACΛ̃ IJ 6

(49)

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00589. Details of the basis set extrapolation, corrections to the experimental ligand dissociation energies, as well as the details of the CCSD(T) computations for Fe(CO)5, Cr(CO)6, and Fe(Cp)2 (PDF)

■

REFERENCES

(1) Cramer, C. J.; Truhlar, D. G. Phys. Chem. Chem. Phys. 2009, 11, 10757−10816. (2) DeYonker, N. J.; Williams, T. G.; Imel, A. E.; Cundari, T. R.; Wilson, A. K. J. Chem. Phys. 2009, 131, 024106. (3) Jiang, W.; DeYonker, N. J.; Determan, J. J.; Wilson, A. K. J. Phys. Chem. A 2012, 116, 870−885. (4) Xu, X.; Zhang, W.; Tang, M.; Truhlar, D. G. J. Chem. Theory Comput. 2015, 11, 2036−2052. (5) Crawford, T. D.; Schaefer, H. F. Rev. Comp. Chem. 2000, 14, 33− 136. (6) Bartlett, R. J.; Musiał, M. Rev. Mod. Phys. 2007, 79, 291−352. (7) Shavitt, I.; Bartlett, R. J. Many-Body Methods in Chemistry and Physics; Cambridge University Press: Cambridge, UK, 2009. (8) DeYonker, N. J.; Allen, W. D. J. Chem. Phys. 2012, 137, 234303. (9) DeYonker, N. J.; Halfen, D. T.; Allen, W. D.; Ziurys, L. M. J. Chem. Phys. 2014, 141, 204302. (10) Coriani, S.; Haaland, A.; Helgaker, T.; Jørgensen, P. ChemPhysChem 2006, 7, 245−249. (11) Scheiner, A. C.; Scuseria, G. E.; Rice, J. E.; Lee, T. J.; Schaefer, H. F. J. Chem. Phys. 1987, 87, 5361−5373. (12) Salter, E. A.; Trucks, G. W.; Bartlett, R. J. J. Chem. Phys. 1989, 90, 1752−1766. (13) Gauss, J.; Stanton, J. F.; Bartlett, R. J. J. Chem. Phys. 1991, 95, 2623−2638. (14) Gauss, J.; Lauderdale, W. J.; Stanton, J. F.; Watts, J. D.; Bartlett, R. J. Chem. Phys. Lett. 1991, 182, 207−215. (15) Jiang, W.; DeYonker, N. J.; Wilson, A. K. J. Chem. Theory Comput. 2012, 8, 460−468. (16) Kutzelnigg, W. J. Chem. Phys. 2006, 125, 171101. (17) Simmonett, A. C.; Wilke, J. J.; Schaefer, H. F.; Kutzelnigg, W. J. Chem. Phys. 2010, 133, 174122. (18) Sokolov, A. Y.; Wilke, J. J.; Simmonett, A. C.; Schaefer, H. F. J. Chem. Phys. 2012, 137, 054105. (19) Sokolov, A. Y.; Simmonett, A. C.; Schaefer, H. F. J. Chem. Phys. 2013, 138, 024107. (20) Sokolov, A. Y.; Schaefer, H. F. J. Chem. Phys. 2013, 139, 204110. (21) Sokolov, A. Y.; Schaefer, H. F.; Kutzelnigg, W. J. Chem. Phys. 2014, 141, 074111. (22) Kutzelnigg, W.; Mukherjee, D. J. Chem. Phys. 1997, 107, 432− 449. (23) Mazziotti, D. A. Chem. Phys. Lett. 1998, 289, 419−427. (24) Mazziotti, D. A. Phys. Rev. A: At., Mol., Opt. Phys. 1998, 57, 4219−4234. (25) Kutzelnigg, W.; Mukherjee, D. J. Chem. Phys. 1999, 110, 2800− 2809. (26) Colmenero, F.; Valdemoro, C. Phys. Rev. A: At., Mol., Opt. Phys. 1993, 47, 979−985. (27) Nakatsuji, H.; Yasuda, K. Phys. Rev. Lett. 1996, 76, 1039−1042. (28) Mazziotti, D. A. Phys. Rev. Lett. 2006, 97, 143002. (29) Kollmar, C. J. Chem. Phys. 2006, 125, 084108. (30) DePrince, A. E.; Mazziotti, D. A. Phys. Rev. A: At., Mol., Opt. Phys. 2007, 76, 042501. (31) Mullinax, J. W.; Sokolov, A. Y.; Schaefer, H. F. J. Chem. Theory Comput. 2015, 11, 2487−2495. (32) Copan, A. V.; Sokolov, A. Y.; Schaefer, H. F. J. Chem. Theory Comput. 2014, 10, 2389−2398. (33) Whitten, J. L. J. Chem. Phys. 1973, 58, 4496−4501. (34) Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. J. Chem. Phys. 1979, 71, 3396−3402. (35) Vahtras, O.; Almlöf, J.; Feyereisen, M. W. Chem. Phys. Lett. 1993, 213, 514−518.

(40)

where P+(IJ, AB)[f(I, J, A, B)] ≡ f(I, J, A, B) + f(J, I, B, A) is the symmetric permutation operator that acts only on the terms in the square brackets. The modified Fock matrix elements F̃QP are defined as FĨ =

Article

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest. 4841


Article

Journal of Chemical Theory and Computation (36) Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Chem. Phys. Lett. 1993, 208, 359−363. (37) Rendell, A. P.; Lee, T. J. J. Chem. Phys. 1994, 101, 400−408. (38) Weigend, F. Phys. Chem. Chem. Phys. 2002, 4, 4285−4291. (39) Kutzelnigg, W. J. Chem. Phys. 1982, 77, 3081−3097. (40) Coleman, A. J. Rev. Mod. Phys. 1963, 35, 668−686. (41) Garrod, C.; Percus, J. K. J. Math. Phys. 1964, 5, 1756−1776. (42) Weinhold, F. J. Chem. Phys. 1967, 47, 2298−2311. (43) Davidson, E. R. Chem. Phys. Lett. 1995, 246, 209−213. (44) Bozkaya, U.; Sherrill, C. D. J. Chem. Phys. 2013, 139, 054104. (45) Bozkaya, U. J. Chem. Phys. 2011, 135, 224103. (46) Sherrill, C. D.; Krylov, A. I.; Byrd, E. F.; Head-Gordon, M. J. Chem. Phys. 1998, 109, 4171−4181. (47) Scuseria, G. E.; Schaefer, H. F. Chem. Phys. Lett. 1987, 142, 354−358. (48) Kutzelnigg, W.; Mukherjee, D. J. Chem. Phys. 2004, 120, 7350− 7368. (49) Pulay, P. Chem. Phys. Lett. 1980, 73, 393−398. (50) Pulay, P. J. Comput. Chem. 1982, 3, 556−560. (51) Manby, F. R. J. Chem. Phys. 2003, 119, 4607−4613. (52) Bozkaya, U. J. Chem. Theory Comput. 2014, 10, 2371−2378. (53) Bozkaya, U.; Sherrill, C. D. J. Chem. Phys. 2016, 144, 174103. (54) Bozkaya, U. Phys. Chem. Chem. Phys. 2016, 18, 11362−11373. (55) Bozkaya, U. J. Chem. Theory Comput. 2016, 12, 1179−1188. (56) Epifanovsky, E.; Zuev, D.; Feng, X.; Khistyaev, K.; Shao, Y.; Krylov, A. I. J. Chem. Phys. 2013, 139, 134105. (57) Riplinger, C.; Neese, F. J. Chem. Phys. 2013, 138, 034106. (58) Weigend, F.; Köhn, A.; Hättig, C. J. Chem. Phys. 2002, 116, 3175−3183. (59) Jurečka, P.; Šponer, J.; Č erný, J.; Hobza, P. Phys. Chem. Chem. Phys. 2006, 8, 1985−1993. (60) Ř ezác,̌ J.; Hobza, P. J. Chem. Theory Comput. 2013, 9, 2151− 2155. (61) Soydaş, E.; Bozkaya, U. J. Chem. Theory Comput. 2013, 9, 1452− 1460. (62) Ten-No, S.; Manby, F. R. J. Chem. Phys. 2003, 119, 5358−5363. (63) Schütz, M.; Manby, F. R. Phys. Chem. Chem. Phys. 2003, 5, 3349−3358. (64) Aquilante, F.; Lindh, R.; Pedersen, T. B. J. Chem. Phys. 2007, 127, 114107. (65) Manby, F. R.; Knowles, P. J.; Lloyd, A. W. J. Chem. Phys. 2001, 115, 9144−9148. (66) Turney, J. M.; Simmonett, A. C.; Parrish, R. M.; Hohenstein, E. G.; Evangelista, F. A.; Fermann, J. T.; Mintz, B. J.; Burns, L. A.; Wilke, J. J.; Abrams, M. L.; Russ, N. J.; Leininger, M. L.; Janssen, C. L.; Seidl, E. T.; Allen, W. D.; Schaefer, H. F.; King, R. A.; Valeev, E. F.; Sherrill, C. D.; Crawford, T. D. WIREs Comput. Mol. Sci. 2012, 2, 556−565. (67) Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (68) Hättig, C. Phys. Chem. Chem. Phys. 2005, 7, 59−66. (69) Weigend, F. J. Comput. Chem. 2008, 29, 167−175. (70) Neese, F.; Valeev, E. F. J. Chem. Theory Comput. 2011, 7, 33−43. (71) Jonas, V.; Thiel, W. J. Chem. Phys. 1995, 102, 8474−8484. (72) Grimme, S. J. Chem. Phys. 2003, 118, 9095−9102. (73) Decker, S. A.; Klobukowski, M. J. Am. Chem. Soc. 1998, 120, 9342−9355. (74) Furche, F.; Perdew, J. P. J. Chem. Phys. 2006, 124, 044103. (75) Burow, A. M.; Bates, J. E.; Furche, F.; Eshuis, H. J. Chem. Theory Comput. 2014, 10, 180−194. (76) Delley, B.; Wrinn, M.; Lüthi, H. P. J. Chem. Phys. 1994, 100, 5785−5791. (77) Ehlers, A. W.; Frenking, G. Organometallics 1995, 14, 423−426. (78) van Wüllen, C. J. Chem. Phys. 1996, 105, 5485−5493. (79) Barnes, L. A.; Liu, B.; Lindh, R. J. Chem. Phys. 1993, 98, 3978− 3989. (80) Li, J.; Schreckenbach, G.; Ziegler, T. J. Phys. Chem. 1994, 98, 4838−4841. (81) Klopper, W.; Lüthi, H. P. Chem. Phys. Lett. 1996, 262, 546−552.

(82) Harding, M. E.; Metzroth, T.; Gauss, J.; Auer, A. A. J. Chem. Theory Comput. 2008, 4, 64−74. (83) Vancoillie, S.; Zhao, H.; Tran, V. T.; Hendrickx, M. F. A.; Pierloot, K. J. Chem. Theory Comput. 2011, 7, 3961−3977. (84) Phung, Q. M.; Vancoillie, S.; Pierloot, K. J. Chem. Theory Comput. 2012, 8, 883−892. (85) Xu, Z.-F.; Xie, Y.; Feng, W.-L.; Schaefer, H. F. J. Phys. Chem. A 2003, 107, 2716−2729. (86) Braga, D.; Grepioni, F.; Orpen, A. G. Organometallics 1993, 12, 1481−1483. (87) Beagley, B.; Schmidling, D. G. J. Mol. Struct. 1974, 22, 466−468. (88) Shen, J.; Gao, Y.; Shi, Q.; Basolo, F. Inorg. Chem. 1989, 28, 4304−4306. (89) Lewis, K. E.; Golden, D. M.; Smith, G. P. J. Am. Chem. Soc. 1984, 106, 3905−3912. (90) Takusagawa, F.; Koetzle, T. F. Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1979, 35, 1074−1081. (91) Jost, A.; Rees, B.; Yelon, W. B. Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1975, 31, 2649−2658. (92) Graham, J. R.; Angelici, R. J. Inorg. Chem. 1967, 6, 2082−2085. (93) Rees, B.; Mitschler, A. J. Am. Chem. Soc. 1976, 98, 7918−7924. (94) Haaland, A.; Nilsson, J. E.; Olson, T.; Norin, T. Acta Chem. Scand. 1968, 22, 2653−2670. (95) Haaland, A.; Lusztyk, J.; Novak, D. P.; Brunvoll, J.; Starowieyski, K. B. J. Chem. Soc., Chem. Commun. 1974, 54−55. (96) Ryan, M. F.; Eyler, J. R.; Richardson, D. E. J. Am. Chem. Soc. 1992, 114, 8611−8619. (97) Balabanov, N. B.; Peterson, K. A. J. Chem. Phys. 2005, 123, 064107. (98) Schütz, M.; Werner, H.-J. J. Chem. Phys. 2001, 114, 661−681. (99) Werner, H.-J.; Manby, F. R.; Knowles, P. J. J. Chem. Phys. 2003, 118, 8149−8160. (100) Riplinger, C.; Sandhoefer, B.; Hansen, A.; Neese, F. J. Chem. Phys. 2013, 139, 134101. (101) McClelland, B. W.; Robiette, A. G.; Hedberg, L.; Hedberg, K. Inorg. Chem. 2001, 40, 1358−1362.

4842


Spin-Adapted Formulation and Implementation of Density Cumulant

Recommend Documents