Excited States of Large Open-Shell Molecules: An Efficient, General

Article pubs.acs.org/JPCA

Excited States of Large Open-Shell Molecules: An Efficient, General, and Spin-Adapted Approach Based on a Restricted Open-Shell Ground State Wave function Michael Roemelt and Frank Neese* Max Planck Institut für Chemische Energiekonversion, Stiftstrasse 34-36, D-45470 Mülheim an der Ruhr, Germany S Supporting Information *

ABSTRACT: A spin-adapted configuration interaction with singles method that is based on a restricted open-shell reference function (ROCIS) with general total spin S is presented. All excited configuration state functions (CSFs) are generated with the aid of a spin-free second quantization formalism that only leads to CSFs within the first order interacting space. By virtue of the CSF construction, the formalism involves higher than singly excited determinants but not higher than singly excited configurations. Matrix elements between CSFs are evaluated on the basis of commutator relationships using a symbolic algebra program. The final equations were, however, hand-coded in order to maximize performance. The method can be applied to fairly large systems with more than 100 atoms in reasonable wall-clock times and also parallelizes well. Test calculations demonstrate that the approach is far superior to UHF-based configuration interaction with single excitations but necessarily falls somewhat short of quantitative accuracy due to the lack of dynamic correlation contributions. In order to implicitly account for dynamic correlation in a crude way, the program optionally allows for the use of Kohn−Sham orbitals in combination with a modest downscaling of two-electron integrals (DFT/ROCIS). All two-electron integrals of Kohn−Sham orbitals that appear in the Hamiltonian matrix are reduced by a total of three scaling parameters that are suitable for a wide range of molecules. Test calculations on open-shell organic radicals as well as transition metal complexes demonstrate the wide applicability of the method and its ability to calculate the electronic spectra of large molecular systems.

■

INTRODUCTION The description of electronically excited states of molecules has been a major goal of quantum chemistry over the last decades.1−3 It is commonly assumed that configuration interaction with single excitations (CIS) can be regarded as a valid zeroth order approximation to excited states. However, for an accurate description of transition energies and transition moments, the inclusion of higher excitations relative to the Hartree−Fock (or multireference) ground state wave function and the excited state wave function either by configuration interaction (CI)4−7 or perturbation theory (PT)8−12 or coupled cluster (CC) linear response13−17 or equation of motion18,19 schemes is inevitable. These higher excitations serve the dual purpose of compensating for the bias of the calculation for the ground state as the orbitals have been optimized for this state and, second, bring in the differential dynamic correlation between the ground and excited states. Thus, a number of single-reference methods exist that provide different trade-offs between computational cost and accuracy. The success of all of those methods is based on the fact that the CIS wave function often provides a qualitatively correct description of the excited states under consideration. If this is not the case, the inclusion of dynamic correlation cannot compensate for the qualitative errors inherent in the CIS approximation. For example, excited states that are © 2013 American Chemical Society

dominated by double excitations from the Hartree−Fock (HF) reference cannot successfully be described in this way. Alternatively, time-dependent density functional theory (TD-DFT)1,20,21 provides a computationally affordable method for the calculation of transition energies and moments on the basis of a Kohn−Sham determinant. Within the widely used Tamm−Dancoff approximation (TDA)22 it can be regarded as an equivalent of CIS theory in the framework of DFT. However, unlike CIS, TD-DFT is formally exact (however, see Schirmer and Dreuw for arguments against the formal exactness of TDDFT23) and includes electron correlation through the exchange-correlation kernel (XC-kernel). In practice, the adiabatic approximation is almost invariably made in which the XC-kernel is taken to be time independent. Together with the many approximations that enter the specific functionals that are presently used, TD-DFT in its realization is a highly approximate method with many severe qualitative failures.24,25 It is, however, under most circumstances, still more accurate than CIS. A second way to apply DFT to the calculation of electronically excited states is realized in the DFT/SCI and DFT/MRCI methods that have been proposed by Grimme and co-workers.26,27 Received: December 21, 2012 Revised: March 14, 2013 Published: March 19, 2013 3069

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

In the present work, an efficient method is reported that is based on a single high-spin ROHF ground state determinant of spin S. The excited state expansion space is spanned by the spatially singly excited and properly spin-coupled CSFs of the same total spin, including the single member of the “trip-(2S+1)” excitation space that belongs to the first-order interacting space. This is achieved through the application of suitable spin-adapted excitation operators to the ground state determinant. Furthermore, we report a new parametrization scheme for the usage of Kohn− Sham orbitals in a ROCIS calculation in the spirit of Grimme’s DFT/MRCI.27 The resulting DFT/ROCIS approach allows for the implicit description of electron correlation while keeping the formal properties of a regular configuration interaction method. The implementation of ROCIS and DFT/ROCIS is of production level type and fairly general. It includes the ability to calculate absorption (ABS) and circular dichroism (CD) spectra. Further extension of the approach to magnetic properties and core-level excitations will be described elsewhere. The program is part of the ORCA suite of electronic structure programs41 and is fully parallelized.

They feature a conventional configuration interaction calculation based on Kohn−Sham molecular orbitals. Since twoelectron integrals over Kohn−Sham orbitals are generally too large, they are scaled down by an empirical parametrization scheme. Excellent results that compare well with highly correlated multireference ab initio methods for excitation energies of organic molecules have been obtained with the DFT/MRCI method.28 As has been known for a long time,29,30 and reviewed by Dreuw and Head-Gordon,1,31 there exist two possibilities for a CIS theory depending on the type of open-shell HF treatment of the ground state. If the electronic ground state is described by an unrestricted Slater-determinant (UHF), the corresponding treatment of excited states would be an unrestricted CIS (UCIS). Accordingly, if a restricted high-spin open-shell Hartree−Fock (ROHF) determinant is chosen as reference, the restricted open-shell configuration interaction with singles (ROCIS) method for excited states arises. Obviously, the ROHF method could be used to describe more complicated open-shell situations.32 However, tailored excited state treatments based on such more general open-shell situations do not appear to have been developed. In the present work, we will also tacitly assume the high-spin open-shell case. For ROCIS, the question arises as to what is being understood as a “single excitation”. If substitutions are made at the spin−orbital level, the manifold of excited configuration state functions (CSFs) contains only single determinants. This space is insufficient for the construction of eigenfunctions of the total spin, even if the ground state determinant is a pure spin eigenfunction. Alternatively, substitutions can be made at the orbital level followed by proper coupling of unpaired electrons to spin eigenfunctions. Such a method has, for example, been worked out by Zerner and coworkers in the framework of the semiempirical INDO/S method.29,30,33−35 Obviously, the space of excited determinants that is required to provide properly spin adapted excited CSFs includes double excitations at the spin−orbital level. Hence, the excited CSFs are not orthogonal to the ground state through the nonrelativistic (Born−Oppenheimer) Hamiltonian and hence after diagonalization there is a (small) correction to the ground state coming from those excited CSFs. In general, these CSFs can be thought of as representing triplet excitations out of the doubly occupied molecular orbitals (DOMOs) that are coupled to spin-flips in the singly occupied molecular orbitals (SOMOs). For example, if the total spin after spin-coupling is a doublet state, these excitations have been referred to as “trip-doublets”. Gouterman and co-workers have identified excited states that are dominated by such excitations early on in studies of porphyrins.36 In our own work, such states were met, for example, in metal−radical systems.37 Cory and Zerner have provided a highly informative review on the subject.38 In the ab initio field, later developments by Head-Gordon and co-workers have discussed various extensions of the ROCIS expansion space with doubly excited determinants relative to the ROHF ground state determinant. This has resulted in the ROCIS(D)12 and XCIS39 methods. One potential complication is that the number of linearly independent spin-couplings for a configuration of a given total spin S increases with the number of open shells. However, it is important to realize that only one member of the CSF set of a given spatial configuration belongs to the first-order interacting space (FOIS).40 Hence, through proper construction of the single (spatial) excitations, it should be possible to only include this member in the expansion space. It appears that such a method has not yet been developed despite the fact that it might serve as a good starting point for the more accurate description of open-shell excited states.

■

THEORY ROCIS. The general spin−orbital based CIS wave function is a linear combination of all singly excited determinants relative to the Hartree−Fock reference determinant: |ΨHF⟩ = |0⟩ = |ϕ1ϕ2···ϕi···ϕn| |ΨCIS⟩ =

∑ cia|Φia⟩ (1.1)

ia

Here, is the amplitude of the singly excited determinant |Φai ⟩, which is obtained by replacing spin−orbital ϕi by the formerly unoccupied spin−orbital ϕa. In a second quantization approach, this excited determinant is obtained by applying the usual Fermion excitation operator âai = â†a âi to the HF determinant: cai

|Φia⟩ = aî a|ΨHF⟩

(1.2)

For the ROCIS method, the reference function is the solution of the ROHF equations. The orbital set can be divided into three groups: (1) the group of doubly occupied molecular orbitals (DOMOs, labeled by i,j,...), (2) the group of singly occupied orbitals (SOMOs, labeled by t,u,v,w), (3) the group of empty orbitals (VMOs, labeled by a,b,...). In the following, orbitals from any subset are labeled with indices p,q,r,s. In order to define excited CSFs it is most convenient to use qβ spin traced excitation operators Epq = âqα pα+âpβ. This leads to a spin-free formalism in which the ROCIS excitation space is divided into three classes of (normalized) CSFs: |Φit ⟩ = Eit |0⟩

DOMO → SOMO

|Φat ⟩

SOMO → VMO

=

|Φia⟩ =

Eta|0⟩

1 a Ei |0⟩ DOMO → VMO 2

(1.3)

It has been shown that the effectiveness of this approach is superior to the UCIS approach but overall inferior to the performance of CIS for excited states of closed shell molecules.31 This was the motivation for the introduction of the XCIS approach which includes “trip-doublets.39 Still, the implementation of the XCIS approach was only able to treat spin-doublets. Zerner’s earlier program could also only handle spin-doublets.42 3070

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

In the present framework, two additional classes of double excitations are included in the excitation space in order to describe such excitations for general spin multiplicities: |Φatwi⟩

EwaEit |0⟩,

=

|Φat ti ⟩

1 a 2 a t = Ei |0⟩ − Et Ei |0⟩ 6 6

If it is assumed that the reference function is a Slater-determinant incorporating orthonormal orbitals, two additional rules apply:

Etu|0⟩ = δtu|0⟩

t≠w

⟨0|Eut = ⟨0|δtu

In cases where the reference function is more complicated (e.g., complete active space self-consistent field, CASSCF, reference functions), dealing with the “active” labels t, u, v, and so on leads to terms involving density matrices of nth order. Since the present work only focuses on the case of a single Hartree−Fock reference function, we will not go in any detail. Using eqs 1.7 and 1.8, it is possible to reduce any vacuum matrix element ⟨0|Epq...Ers|0⟩ to sums involving only Kronecker δ’s. The number of terms arising in the evaluation of operator strings rapidly increases with the number of strings in the operator. Hence, the rules in eqs 1.7 and 1.8 together with the strategy to evaluate matrix elements using them have been incorporated into a formula generation program (COMMUTATOR). The program can easily handle thousands of terms simultaneously in calculating the matrix elements of the BO Hamiltonian. Note that this approach does not make explicit reference to the form of the reference function (e.g., form of the spatial or spin parts) and is thus applicable to a wide range of reference functions. DFT/ROCIS. Pure CIS methods such as the presented one completely neglect dynamic electron correlation. In many cases, this leads to significant deviations of predicted excitation energies and intensities from the experiment (vide inf ra). Grimme and co-workers developed DFT/SCI and DFT/MRCI methods that introduce dynamic electron correlation to a configuration interaction calculation implicitly by using Kohn− Sham orbitals.26,27 The Kohn−Sham orbitals inherently include some dynamic correlation due to the exchange-correlation potential that occurs in the Kohn−Sham operator. A similar approach was chosen to improve the presented method and is hence referred to as DFT/ROCIS approach. The ROHF reference determinant is replaced by a restricted open-shell Kohn−Sham (ROKS) determinant. It incorporates two sets of orbitals: the DOMO’s and SOMO’s, respectively. They are eigenvectors of the two Kohn−Sham matrices FC(KS) and FO(KS):

(1.4)

The first class represents triplet excitations from DOMOs to VMOs in conjunction with an excitation within the SOMO space that reduces the number of singly occupied MOs in this space by one. The second class specifically describes triplet excitations coupled to a spin-flip in a SOMO. The latter are the “classical” trip-doublet excitations. Obviously, these classes have nonvanishing matrix elements ⟨Φatwi|HBO|0⟩ since they are not subject to Brillouin’s theorem. In the derivations below, we generally do not assume the Brillouin conditions to be satisfied, such that the formulas derived below are also applicable to cases where the reference function does not consist of HF orbitals. Accordingly, it is then possible for the HF reference function to enter the total ROCIS wave function of excited states and, on the other hand, for excited determinants to enter the ground state wave function. However, in our experience, these mixings are only on the order of a few percent (Supporting Information). Taken together, the ROCIS wave function is written as |ΨROCIS⟩ = c0|0⟩ +

∑ cit|Φit ⟩ + ∑ cta|Φat ⟩ + ∑ cia|Φia⟩ it

+

∑

ta

cwiat |Φatwi⟩

+

∑

itwa

ia

ctiat |Φat ti ⟩

ita

t≠w

(1.5)

As described below, transition energies for excited states of open-shell molecules are obtained as the solution to the eigenvalue problem of the Hamiltonian matrix in the basis spanned by the five excitation classes included in the ROCIS wave function and the reference determinant. Calculation of Matrix Elements. The spin-free second quantization form of the BO-Hamiltonian is (SQ) HBO =

∑ hpqEqp + p,q

1 2

∑

(1.8)

(pq|rs)(EqpEsr − δqrEsp)

p,q,r ,s

(1.6)

DOMOs

where hpq = ⟨p|ĥ|q⟩ are integrals over the one-electron part of the Hamiltonian and (ps|rq) = ∫ ϕp*(1)ϕr*(2)(1/r12)ϕs(1)ϕq(2) dr1 dr 2 are two-electron repulsion integrals.43 The spin-free replacement operators Epq are defined above. Provided that the one- and two-electron integrals are available, evaluation of matrix elements over the Hamiltonian reduces to the evaluation of strings of replacement operators between the reference function, denoted |0⟩. The evaluation of all such strings can be performed on the basis of commutation rules of the replacement operators together with some simple additional relations:

SOMOs

+

=

∑

SOMOs

+

∑

XC 2(J ii )pq − c HF(K ii)pq − c DFVpq [ρ]

i XC (J tt )pq − c HF(K tt )pq − c DFVpq [ρ]

t C(KS) = F pq −

1 XC {c HF(K tt )pq − c DFVpq [ρ]} 2 (1.9)

Here, ρ] = ⟨ϕp|(δEXC[ρ])/(δρ)|ϕq⟩ denotes elements of the exchange-correlation matrix that is functional specific.21 Both matrices are given in the hybrid density functional form, since hybrid density functionals will be used in all presented DFT/ROCIS calculations. The coefficients cHF and cDF denote VXC pq [

=0

[Eqp ; Esr ] = δqrEsp − δpsEqr

1 XC {c HF(K tt )pq − c DFVpq [ρ]} 2

DOMOs

= 2δpi|0⟩

⟨0|Epa

(J tt )pq −

∑

O(KS) F pq = hpq +

⟨0|Eip = ⟨0|2δpi Eap|0⟩

XC 2(J ii )pq − c HF(K ii)pq c DFVpq [ρ]

i

t

⟨0|0⟩ = 1 Epi|0⟩

∑

C(KS) F pq = hpq +

(1.7) 3071

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

We note in passing that the numbers obtained for the BhLYP functional48 are quite similar to the B3LYP45,46 parameters (c1 = 0.20, c2 = 0.30, and c3 =0.40). Only the c2 parameter differs considerably, which is due to the different admixtures of Hartree−Fock exchange in the functionals. Although the results obtained with BhLYP for transition metal L-edges are of comparable accuracy than the equivalent B3LYP calculations,44 the results for optical spectra are significantly inferior (see Supporting Information). As described above, the first parameter scales down Coulomb terms in all diagonal elements of the CI matrix, as does the equivalent parameter in Grimme’s original work (cDFT‑SCI = 1 0.317). Interestingly, the numerical value that we obtained differs from the original one. This difference is due to two reasons: (a) The two sets of parameters were optimized with respect to different training sets and electronic excitations and (b) in the original scheme, cDFT‑SCI is also used to scale down 1 Coulomb terms in off-diagonal elements of the CI matrix, whereas the off-diagonal matrix elements are treated by c3 in the present scheme. With regard to reason b, it appears plausible, that the numerical value of cDFT‑SCI has a value intermediate 1 between our present parameters c1 and c3. It has to be noted that, while in the original parametrization scheme off-diagonal exchange terms are left untouched, they are scaled by c3 in the present scheme. Moreover, off-diagonal terms that involve the Fock matrix (e.g., δijFCab−δabFCji ) are scaled by c3 too. This emphasizes the diagonal-dominated character of the DFT/ROCIS matrix. A further difference between the original and the present scheme is the handling of the diagonal exchange terms. In Grimme’s original work, they are modified by an additive term,26 while in the present scheme they are multiplied with c2. The scaling of the exchange terms results in more accurate relative energies of triplet and singlet excitations, which is essential for the calculation of transition metal L-edge spectra. In Grimme’s work on multireference density functional theory-based CI, the exchange terms (diagonal and off-diagonal) are treated in a similar fashion as in the present work.27 Test calculations on a training set of organic closed-shell molecules (see Supporting Information) demonstrate that the DFT/MRCI method is, not surprisingly, better suited to evaluate excitation energies of such molecules than DFT/ROCIS. This is mainly due to two reasons: (a) The empirical parameters that enter the DFT/MRCI method were optimized with respect to a similar test set of valence excited states of organic closed shell molecules, while the parameters of DFT/ROCIS were optimized with respect to transition metal L-edge spectra (vide supra). (b) The DFT/MRCI treatment is based on a more flexible wave function that incorporates higher excitations and hence also treats states with predominantly double excitation character correctly (see Supporting Information). Such double excitations are known to be of importance in extended π-systems of linear polyenes. However, due to its more restricted expansion space, DFT/ ROCIS is applicable to large molecules that are out of reach for DFT/MRCI and highly correlated ab initio methods. Furthermore, DFT/ROCIS is able to treat molecules of arbitrary spin states, whereas the DFT/MRCI parameters have only been fitted for singlet and triplet states.

the amount of HF exchange and pure density functional related exchange-correlation contribution, respectively. Thus for a HF calculation, cDF = 0. It was found that the two-electron integrals that occur in the CI matrix are too large when Kohn−Sham orbitals are used.26,27 According to the nature of the integrals and their position in the CI matrix, the integrals are thus scaled down by the parameters c1, c2, and c3. This is in accordance with the idea of a “screened” electron−electron interaction as would arise from a proper dynamic correlation treatment. Obviously one can only hope to capture these effects in an average way through a crude parametrization of the present form. The parameters c1 and c2 reduce the Coulomb- and exchange-like terms in the diagonal entries of the CI matrix, whereas c3 is used to scale down all off-diagonal elements. In the subspace spanned by the |Ψai ⟩ basis functions, the CI matrix elements then become C(KS) HiaDFT/ROCIS = Faa − FiiC(KS) − c1(ii|aa) + 2c 2(ia|ia) , ia C HiaDFT/ROCIS = c3{δijFab − δabF jiC − (ij|ab) + 2(ia|jb)} , jb

(1.10)

Formulas for the remaining diagonal elements are found in the Supporting Information. Note that the scaling of the off-diagonal elements as introduced in eq 1.10 implies not only the scaling of two-electron integrals but also of off-diagonal Fock-terms. This scheme has been chosen for the sake of internal consistency. In particular, these Fock terms contain two electron integrals over active labels. It depends on whether one bases a formula on FC(KS) or FO(KS) whether these two electron integrals occur as a correction to F outside the formula or inside F. Not scaling the Fock terms in the same way would have led to inconsistent treatment of these two-electron integrals. In general, the DFT/ROCIS approach accounts for the fact that virtual Kohn−Sham orbitals are better suited to represent electronically excited states than their HF counterparts. In order to not partially relax these Kohn−Sham orbitals back to HF orbitals, the Fock matrix terms in the matrix elements of the form ⟨0|Ĥ |Φqp⟩ are discarded. The presented ROCIS and DFT/ROCIS methods were motivated by the desire to provide the nonrelativistic basis for the calculation of transition metal L-edge X-ray absorption spectra. Consequently, the three parameters were optimized with respect to a test set of L-edges of first-row transition metal complexes in different oxidation and spin states. This work is the subject of another publication where the parametrization process will be described in detail.44 For the purposes of the present paper, it is sufficient to state that the parameters that we use have not been optimized for valence states of open-shell molecules. The values obtained are c1 = 0.18, c2 = 0.20, and c3 = 0.40 in conjunction with the B3LYP functional.45,46 The applications presented below support the assessment that a single set of parameters is reasonably successful for the prediction of both, core- and valence spectra. At this point it should be noted that the set of parameters was optimized utilizing the B3LYP functional45,46 in conjunction with the def2-TZVP(-f) basis set.47 A change of the functional and, to a much lesser extent, the basis set necessitates the use of different parameters to give accurate results. In fact, one could improve the achievable accuracy slightly by reoptimizing the parameters. We do, however, prefer to not do so since any of such reparameterizations does not revise the physical content of the calculations but only serves to improve the apparent numerical accuracy.

■

IMPLEMENTATION Direct CI. The ROCIS method sketched above has been implemented in a development version of the ORCA program package.41 Excited state energies and wave functions are obtained as eigenvalues and eigenvectors of the CI matrix in the

3072

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

notation. In this way the bookkeeping during the integral transformation is greatly simplified. The integrals (ij|kl) and (ik|jl) enter the expression of three different types of Fock matrices that are used during the evaluation of matrix elements of the BO-Hamiltonian: the closed-shell Fock matrix (FC), the open-shell Fock matrix (FO), and an intermediate Fock matrix (FI), which is only defined for convenience:

ROCIS multielectron-function space. They are determined in a direct CI fashion using a Davidson diagonalization routine.49 This requires the calculation of the σ-vector: σI = (H − SE I)cI

(1.11)

S denotes the overlap matrix of the ROCIS expansion space and cI is a vector. Once the CI problem has been solved by an appropriate vector cI, the σ-vector has only vanishing entries. The existing implementation of the Davidson algorithm in the ORCA program assumes a fully orthonormal CI basis,50 thus requiring the overlap matrix to equal unity. A closer look at the ROCIS expansion space given in eqs 1.3 and 1.4 reveals that this is generally not the case. If the reference determinant contains more than one SOMO, off-diagonal elements arise in the overlap matrix in the ⟨Φatti |Φbu uj ⟩ block. For each pair of (i,a) labels, the overlap of ROCIS expansion functions Φatti and Φau ui , ⟨Φatti |Φau ui ⟩ = (1/3) with t ≠ u, does not vanish (all other offdiagonal elements of S do vanish). Thus the expansion space has to be transformed such that the overlap matrix becomes the unit matrix (S = 1). This is achieved by using the well-known Löwdin orthonormalization procedure.43,51 Since the formulas for HcI were derived in the ROCIS basis given in eqs 1.3 and 1.4 and are also calculated in this way, two transformations are conducted during each CI iteration step. First, the trial vector c′I is transformed from the orthonormal basis to the nonorthonormal ROCIS basis. Then σI = HcI is calculated followed by back transformation to the orthonormal basis according to σI′ = S−(1/2)σI. The orthonormalization represents no specific problem as the overlap matrix has a sparse block structure: Stu(i , a)

=

au ⟨Φat ti |Φui ⟩

DOMOs

∑

SOMOs

2(J ii )pq − (K ii)pq +

∑

i

(J tt )pq

t

− (K tt )pq SOMOs

=

C F pq

∑

−

t

1 tt (K )pq 2

DOMOs I Fpq

∑

= hpq + C = F pq +

SOMOs ii

ii

2(J )pq − (K )pq +

i SOMOs

∑ t

∑

(J tt )pq

t

1 tt (K )pq 2 (1.14)

pq

pq

Note that (J )rs = (pq|rs) and (K )rs (pr|qs) are four-index integrals. Together with the six four-index integral classes, these Fock operators enable the program to calculate all desired entries of the σ-vector (vide inf ra). Working Equations. Including the reference determinant, the presented ROCIS space consists of six classes of multielectron functions. Thus, the σ-vector is naturally divided into six parts, each of which has to be calculated separately. For simplicity, the reference energy times the overlap is subtracted from the σ-vector. The resulting expressions for the σ-vector are then given by

(1.12)

σ0 =

Stu(i , a) − 1/2σI, iuua

+

2

which only scales as O(N ), where N is the number of basis functions used in the calculation. Integral Transformation and Fock Matrices. The orca_rocis program has been developed as part of the ORCA program package. The orbitals are generated from either a ROHF (ROKS) calculation or from the quasi-restricted orbitals of a UHF (UKS) calculation, which yields a nearly identical determinant.52 Two-electron integrals are transformed to the MO basis. For the transformation step, the resolution of identity (RI) approximation53−56 can be employed. In this case, the three classes (ij|K), (ia|K) and (ab|K) are generated, where K is an orthogonalized auxiliary basis function in the Coulomb metric.57 From these integrals. the four index integral classes (ij|kl), (ij|ka), (ij|ab), (ik|jl), (ik|ja), and (ia|jb) are produced and stored on disk using the parallel file handling and data compression tools available within the ORCA framework. Only for this purpose, do the integral indices i, j, k, and l include the SOMOs in order to minimize the number of different integral classes. For example, the integral class (ij|ab) has some members that would be referred to as (tu|ab) in the chosen

1 2

∑ FjuI cju + ∑ FubOcub + ju

(1.13)

u

(J tt )pq

t

DOMOs O F pq = hpq +

∑ Stu(i ,a) − 1/2c I,′ iuua ∑

∑

i

u

σI,′itta =

SOMOs

2(J ii )pq − (K ii)pq +

1 − (K tt )pq 2

Since the submatrices are not dependent on the (i,a) labels, a single (small) matrix contains all essential information about the nondiagonal part of S and hence about S itself. The transformation of the trial and the σ-vector only involves the Φatti part of the expansion space and is thus carried out as: c I, itta =

∑

C F pq = hpq +

ub

∑ F Cjbcjb jb

1 {∑ (K uu)jb cvjbv 6 juvb

∑ (K uv)jb cvjbu − juvb

∑ 2(K uu)jb cujbu} (1.15)

jub

σit = FtiIc0 −

∑ FjiI cjt + ∑ FtuI ciu + ∑ ((J it )ju − (K it )ju )cju j

u

ju

1 + ∑ (J it )ub cub + {∑ FtbI cib 2 ub b +

∑ (2(J it )jb − (K it )jb )cjb} + ∑ FvbOcvibt jb

+

vb

∑ (J

ut

)vb cvibu

−

∑ (K

uvb

+ +

3073

)jb cujbt

jub

1 {∑ −2FtbOctibt − 6 b

∑ 2(K jb

iu

it

)jb ctibt

+

∑ FtbI cuibu − ∑ 2(K uu)tb cuibu ub

∑ (K jub

ub it

)jb cuibu} (1.16)

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A σta = FatOc0 +

∑ (J ta )ju cju − ∑ FutOcua + ∑ FabOctb ju

+

u

ju

∑ (K ut)jv cvjau juv

+

−

cujau}

1 {− ∑ 2(K uu)jt ctjat + + 6 ju

−

jub

ju

∑ 2(K

ta

)jb ctjbt

−

∑ (K

−

jb

j

∑ 2(K

tt

ta

)jb cujbu}

−

∑

+

− =

+

∑

N u Fau ci

∑ (2(J

+

u

−

ia

)ju − (K

∑ (J

)ab cub

)ju )c ju

+ −

−

+

b

−

− +

juv

∑ (K uu)ab cvibv}

∑ (K tw)ua ciu − ∑ (K iw)ja cjt u

+ FtiIcwa +

+

u

uv

ub

∑ (j

bt )ab cwj

jb

− (K

+

−

ub

∑ ((K

iw

)ju cujat

au )ju cwj )

+ (J

au )ju cwj

−

j

ju

1 O at O aw cti − 2Ftw cwi + ∑ 2(J tw )ab ctibt {2Ftw 6 b 1 tw bw − ∑ 2(J )ab cwi } + {−∑ 2(K wt )ab ctibt 6 b b

∑ 2(J tu )ab ctibu − ∑ 2(K tu)ab cuibt

−

ub

∑ (K

iu

)vj cvjau

−

ub

∑ (K

uv

)ab cvibu

uvb

1 {−∑ 4FijOctjat + 6 j

∑ 4(K tt )ij ctjat − ∑ 2FijOcujau j

ju

∑ 4(K tt )ab ctibt − ∑ 2FabOcuibu + ∑ 2(K tt )ab cuibu} ub

∑ 2(K

tt

)ab cuiau

+

ub

∑ 4(K

tt

)uu cuiau

u

∑ 2(K tt )ab cviav} uvb

∑ (K uu)ij cvjav} +

∑ 4(K uu)ij ctjat ju

1 {−∑ 4(J ij )ab ctjbt 6 jb

∑ 2(J ij )ab cujbu} (1.20)

■

+

DFT/ROCIS As outlined above, DFT/ROCIS calculations do not only use a ROKS reference determinant and Kohn−Sham orbitals to build excited state CSF’s, they also necessitate the parametrization of various terms in the CI matrix. Thus DFT/ROCIS calculations feature two steps during the generation of the σ-vector.

∑ (K wt )ab cuibu − ∑ 2((K uu)tw cuiau − (K uu)tw cwiaw)} ub

)ab cuibt

The most expensive steps of the σ-vector construction scale as O(N4) and are identical to the steps encountered in a closed shell CIS calculation. It has to be stressed that the computational cost of a ROCIS calculation only slightly depends on the number of unpaired electrons, since the number of |Φatti ⟩ and |Φatwi⟩ expansion functions increases with their first and second power, respectively.

∑ (J it )jw cwjaw − ∑ (J it )jw cujau}

j

−

∑ 2(J

jub

1 aw + {∑ 2(K it )jw cwj 6 j

− 2 ∑ (K it )jw ctiat −

∑ 2(K uu)tv ctiav uv

tu

juv it

ju it

−

1 + {∑ 2(K uu)ij cujau − 6 ju

∑ (J tu )ab cwibu − ∑ (J wu )ab cuibt + ∑ (K wu)ab cuibt ji

∑ (K

)uv cviau

ub

j

∑ FuwO cuiat + ∑ FtuI cwiau − ∑ (K tw)uv cviau ub

−

+

∑ (K tw)ji cjb} − ∑ FjiI cwjat j

)ju ctjau

ju tt

1 + {∑ 4(K tt )uu ctiat − 6 u

u

u

∑ 2(J

−

b

∑ (J it )ab cwb − ∑ (J it )wu cua b

−

+

j

1 + {∑ (K tw)ab cib − 2 b

ju it

1 O bt + ∑ 2(K tt )ij cujau} + cti {∑ 4Fab 6 ju b

(1.18)

uvb

O t σwiat = (J it )wa + Faw ci +

∑ 2(K

)ju cujat

juv

∑ 2(K uu)ab cuibu + ∑ (K uu)ji cvjav ub

u it

ub

1 {∑ 2(K uu)ji cujau 6 ju

juv

−

+

uvb

ub

∑ 2FutOcuiat − ∑ 2FutOctiau + ∑ 2(K it )ju ctjau

uv

∑ (4(J ia )jb − 2(J ia )jb )cjb} + ∑ (K uv)ab cvibu ∑ (K uv)ji cvjau +

ju

ju

∑ 2FabC cib

jb

∑ 2(K tt )ab cib + ∑ (K uu)ij cja − ∑ (K uu)ab cib} u

∑ FuiOcua + (2(J ia )ju − (K ia)ju )cub u

+

b

1 + {∑ 2(K tt )ij c ja 2 j

b

ia

ju

1 + {∑ −2F jiCc ja + 2 j

∑ 2(K tt )iu cua − ∑ 2(J it )ab ctb u

iu

ub

) 2FaiCc0

∑ (K iu)ja cju − 2FitI cta

+

ju

FiuOcua

u

jub

(1.17)

2 σia

)ua ciu

u

∑ 2F Ojt ctjat + ∑ 2(K uu)it j

∑ FauI ciu − FatOcit + ∑ 2(K it )ja cjt u

∑ (2(J ta )jb − (K ta)jb )cjb} + ∑ FjuI ctjau + ∑ (J ju )ab ctjbu 1 + {∑ F Ojt cujau + 6 ju

∑ −(K uu)ia c0 − 2(K tt )ia c0 u

1 {∑ F jtOc ja 2 j

∑ ((J ta )ub − (K ta)ub )cub + jb

−

6σtiat =

b

ub

+

Article

u

(1.19) 3074

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

The first step is the evaluation of the σ-vector as shown in eqs 1.15−1.20 using the Fock matrices introduced in eq 1.14. All indices in these equations then refer to Kohn−Sham rather than Hartree−Fock orbitals. Subsequently, the entire σ-vector is multiplied with c3. In order to restore the original size of the contributions to σ that originate from the diagonal part of the CI matrix, the term (1 − c3)HIcI has to be added to each element σI. Accordingly, the σ-vector is divided in two parts according to

All three contributions contain integrals over the respective one-electron operator m = μ, M, Q and multielectron wave functions |I⟩ (initial state wave function) and |F⟩ (final state wave function). Evaluation of these integrals requires the α knowledge of the one-electron integrals mpq = ∫ ϕp*(r) α m̂ (r)ϕq(r)dτ (α = x,y,z) and the one-electron transition density ρ̂IF between the two states. In the second quantization, the transition density reads IF ρpq = ⟨ΨF |Epq|Ψ⟩ I

σI′ = c3 ∑ HIJcJ + (1 − c3)HIIcI J

and hence can be calculated employing the same strategy as for matrix elements over the BO-Hamiltonian. The transition density between two general ROCIS states then becomes

= c3 ∑ HIJcJ + HIIcI (1.21)

J≠I

CI‐basis

The diagonal elements of the CI matrix HII are calculated and stored prior to the σ-vector evaluation and can be used during each iteration of the Davidson diagonalization. In a second step the diagonal contribution HIIcI is altered according to the scheme presented in the Theory section (cf. eq 1.10). Since no summation over the compound label I is performed, this second step is computationally much less demanding than the formation of the σ-vector. The modification of the diagonal contributions involves two different actions. First, the elements of the Fock matrices from eq 1.14 have to be replaced by elements of the Kohn−Sham matrices from eq 1.9. Difference Fock matrices ΔFC = FC − FC(KS), which are generated and stored prior to the σ-vector evaluation simplify this task. The diagonal contribution HIIcI contains Fock matrix elements of the form FCpqcI. Subtraction of the same term using the difference Fock matrix results in the desired replacement of the Fock matrix.

IF ρpq =

Formulas for all 36 possible classes of can be found in the Supporting Information. Once the expression for a certain ⟨λ|Eqp|κ⟩ has been determined, symmetry arguments allow for its usage during the evaluation of the conjugated element ρIF qp as FI well as ρFI and ρ : pq qp CI − basis

ρqpIF =

feq =

1 2 3 α EFI 20

I

∑ i

ρpqFI =

∑

CλICκF⟨λ|Epq|κ ⟩

κ ,λ

∑

CλICκF⟨λ|Epq|κ ⟩

κ ,λ

ρqpFI =

∑ κ ,λ

CI‐basis

CλICκF⟨λ|Eqp|κ ⟩ =

∑

CλFCκI⟨λ|Epq|κ ⟩

κ ,λ

(1.27)

Note that eqs 1.27 are only valid for real valued CI basis functions and coefficients. Additionally the calculated transition energies and intensities enable the orca_mapspc (and, in principle, any other spectra simulation program) utility program to generate a crudely simulated absorption spectrum by convoluting the transitions with Gaussian-type line shape functions with a predefined full-width-athalf-maximum height (fwhm).

■

APPLICATIONS Computational Details. If not otherwise stated, all presented calculations were conducted using the following setup. All calculations were performed with the ORCA program package.41 Geometry optimizations used the BP86 functional61,62 together with the def2-TZVP(-f)47 basis set. Subsequent TDDFT and DFT/ROCIS calculations utilized the B3LYP functional45,46 and the same basis set. During the numerical integration step of all DFT calculations a dense, 320 points Lebedev integration grid (ORCA grid4) was employed. All two-electron integral transformations were carried out applying the RI approximation53−57 with the def2-TZV/J basis set.63 TD-DFT calculations were performed within the TDA approximation.22 DFT/ROCIS calculations employed the above introduced set of parameters. The convergence tolerance for ROCIS and DFT/ROCIS calculations was set to 1.00 × 10−6 eV for state energies and 1.00 × 10−6 for residual norms. Electronic Excitations of Organic Radicals. A series of electronic transition energies of medium-sized and large organic radicals was calculated with the presented ROCIS and DFT/ROCIS methods and compared to results obtained

(1.23)

2

1 ̂ (l (i) + 2s (̂ i)) F 2

x ,y,z

∑ |⟨I |Q ab|F ⟩|2 a,b

CI‐basis

CλFCκI⟨λ|Eqp|κ ⟩ =

κ ,λ

fed =

2 = α 2EFI 3

∑ CI‐basis

Here, ω and I are the angular frequency and time averaged intensity of the incident light beam, respectively. EIF is the energy difference between states |I⟩ and |F⟩, whereas fed, f md and feq denote contributions from the electric dipole (ed), the magnetic dipole (md), and the electric quadrupole (eq) terms. They are defined as59,60

fmd

(1.26)

⟨λ|Eqp|κ⟩

(1.22)

2πω (f + fmd + feq )δ(EIF − ℏω) I ed

2 EFI |⟨I |μ|F ⟩|2 3 2 = α 2EFI |⟨I |M |F ⟩|2 3

Cλ*FCκI⟨λ|Epq|κ ⟩

κ ,λ

A similar strategy is applied for the parametrization of the Coulomb and exchange terms contained in the diagonal contribution. After their original size is restored as described above (eq 1.24) they are reduced by subtracting (1 − c1) times the Coulomb and (1 − c2) times the exchange integral. Note that subtracting means that the (1 − c1) term always has the opposite sign than the original term from eqs 1.15−1.20. Absorption Spectra and Transition Densities. Within the framework of the minimal coupling scheme,58 the absorption intensity for a transition between states |I⟩ and|F⟩ is proportional to the absorption cross section σ, which is given by σ=

∑

CI‐basis

C(DFT) C C F pq CI = F pq cI − ΔF pq cI

(1.25)

(1.24) 3075

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

transition energies. Exclusion of vibronic fine structure in the calculated values leads to a constant overestimation dependent on the conformational shift of the excited state potential energy surface (PES) compared to the ground state PES. Certainly more important is the complete neglect of correlation effects and differential correlation effects in the ROCIS calculations. However, the trend is quite systematic, as can be quantified from the correlation of calculated and experimental data (Figure 1a). Good approximations to the experimental data can be obtained from the calculated values with the linear fit equation

from UCIS, TD-DFT, and experimentally determined 0→0 transition energies. The test set of organic radicals was taken from the study of Dierksen and Grimme.64 It comprises cyclic and polycyclic aromatic and rigid homoaromatic radicals (Chart 1). Chart 1. Lewis Structures of the Series of Organic Radicals That Were Investigated in This Study

ΔE 0 − 0(exp) ≈ 0.769 × ΔEvert(ROCIS) − 0.150 eV (1.28)

The corresponding scaled excitation energies listed in Table 2 match the experimental reference values within a maximum deviation of 0.35 eV. A linear correlation constant of R2 = 0.96 underlines the quality of the fit and hence the systematic nature of the error. The UCIS excitation energies are less accurate (rmsd = 1.24 eV) than the ROCIS excitation energies. As expected, UCIS also has a bias toward too high energies. However, the overestimation is much less systematic as for ROCIS. The linear correlation constant for the linear fit between experimental and calculated excitation energies significantly reduces to R2 = 0.67 (Figure 1b). It can thus be concluded that the spin-adapted ROCIS description of excited states of organic radicals is superior to an unrestricted CIS description due to the inclusion of properly spin adapted CSFs. This finding is in agreement with previous results.31 When the DFT/ROCIS method is employed, the absolute excitation energies of the present test set are reproduced more accurately than with pure ROCIS. The rmsd counts to 0.40 eV, which is significantly lower than the above obained 0.61 eV. In contrast to the ROCIS results, the DFT/ROCIS excitation energies are systematically lower than the experimental values. Despite the better rmsd of the predicted transition energies, the linear fit of experimental versus calculated transition energies (Figure 1c) exhibits a slightly lower correlation constant of R2 = 0.94. Accordingly, the maximum deviation of the scaled DFT/ ROCIS excitation energies in Table 2 is higher than for ROCIS. However, in view of the limited number of test cases, the difference of 0.02 in the R2 value is insignificant. The TDDFT results exhibit a similar systematic error as the pure ROCIS results but include electron correlation implicitly via the static exchange-correlation potential. As a consequence, the calculated excitation energies are too large but closer to the experimental values than the ab initio energies. A rmsd of 0.25 eV reflects the improvement of the predicted absolute transition energies over the ROCIS method. The systematic character of the predicted transition energies is demonstrated by an almost perfect R2 value of 0.99 for the linear fit of experimental versus predicted transition energies (Figure 1c). This excellent behavior of TD-DFT will not hold for larger, more delocalized radicals and will also not hold for open-shell transition metal complexes. [Fe(CO)5] and [Cr(CO)6]. Iron and chromium carbonyl (Chart 2) are typical examples of first row transition metal carbonyl complexes. They both obey the 18 electron rule and have a S = 0 electronic ground state. The correct description of the complicated bonding situation that features ligand-to-metal bonding and metal-to-ligand backbonding poses a challenge to any theoretical method.65−69 While the electronic transitions in the study of organic radicals are of π → π* and n → π* type,

The electronic transitions considered in the course of this investigation are of π → π* and n → π* nature. Table 1 summarizes the calculated and experimentally determined transition energies of the first dipole allowed Table 1. Comparison of Calculated Vertical Transition Energies (in eV) Obtained with ROCIS, DFT/ROCIS and TD-DFT (B3LYP) to Experimentally Determined 0→0 Transition Energies ROCIS UCIS

B3LYP DFT/ B3LYP ROCIS TD-DFT

compound

transition

exp

1 2 3 4 5 6 7 8 9 9 10 10 11

12B2→12A2 12B2→22A2 12B2→22A2 12A2u→12A3g 12B3g→12A2u 12B2→12A2 12A2→12B2 12B2g→12B1u 12Au→12B3g 12Au→12B2g 12A2→12B2 12A2→32A2 12B2→12A2

2.12 4.05 4.04 1.97 1.87 1.53 2.11 1.62 2.01 3.34 1.46 3.16 1.72

2.96 4.93 5.03 2.23 2.33 1.89 2.42 1.68 2.70 3.70 1.85 4.20 1.83

5.16 5.71 5.78 2.48 2.58 2.54 2.1 1.77 3.13 3.5 1.82 4.19 2.52

1.66 3.33 4.20 1.64 1.43 1.20 1.65 1.16 1.62 3.12 1.37 2.92 1.39

2.44 4.14 3.94 2.31 2.05 1.80 2.20 1.95 2.19 3.60 1.76 3.46 1.96

root mean square deviation mean abs error mean error max (+) error max (−) error

0.61 0.52 0.52 1.04

1.24 0.95 0.94 3.04 −0.01

0.40 0.36 −0.33 0.16 −0.76

0.25 0.23 0.22 0.33 −0.1

electronic excitations. Not surprisingly, the data demonstrate that the ROCIS method tends to overestimate the excitation energies of organic radicals. The average deviation from the experimental values is 0.37 eV (root-mean-square deviation (rmsd) = 0.61 eV). This can be partially attributed to the errors arising from comparing vertical transition energies and 0→0 3076

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

Figure 1. Correlation of calculated vertical transition energies ΔEvert and experimentally determined 0→0 transition energies ΔE0−0(exp) for the pure ROCIS method (a), UCIS (b), the DFT/ROCIS method (c), and TD-DFT (d).

Table 2. Scaled Excitation Energies for ROCIS, DFT/ ROCIS, and TD-DFT (B3LYP) compound

transition

1 2 3 4 5 6 7 8 9 9 10 10 11

12B2→12A2 12B2→22A2 12B2→22A2 12A2u→12A3g 12B3g→12A2u 12B2→12A2 12A2→12B2 12B2g→12B1u 12Au→12B3g 12Au→12B2g 12A2→12B2 12A2→32A2 12B2→12A2

experiment ROCISa 2.12 4.05 4.04 1.97 1.87 1.53 2.11 1.62 2.01 3.34 1.46 3.16 1.72

root mean square deviation mean abs error mean error max (+) error max (−) error

Chart 2. Lewis Structures of Iron and Chromium Carbonyl

B3LYP DFT/ B3LYP ROCISb TD-DFTc

2.43 3.94 4.02 1.86 1.94 1.60 2.01 1.44 2.23 3.00 1.57 3.38 1.56

2.02 3.53 4.37 2.01 1.81 1.60 2.01 1.56 1.99 3.37 1.76 3.19 1.78

2.21 4.04 3.83 2.07 1.79 1.52 1.95 1.69 1.94 3.46 1.48 3.31 1.70

0.18 0.16 0.00 0.31 −0.34

0.20 0.13 0.00 0.33 −0.52

0.10 0.09 0.00 0.15 −0.21

thus yields information about the applicability of the two methods to such transitions in complicated cases. The electronic ground state of [Fe(CO)5] features a (3dσ)4 (3dπ)4 (3dz2)0 valence orbital configuration where 3dσ and 3dπ refer to the E(x2 − y2, xy) and the E(xz, yz) orbital manifolds, respectively. In addition to the metal d-orbitals, the valence electronic shell comprises a number of unoccupied CO π*-orbitals with a considerable amount of metal 3d character. Accordingly, the low energy part of the electronic spectrum of [Fe(CO)5] shows two d−d transitions and a number of MLCT transitions. Due to the low resolution of the UV absorption experiment together with the high density of states, it is impossible to extract accurate excitation energies for all states from the experimental spectrum.70 Thus the results of the multireference contracted configuration interaction (MR-CCI)71 study of the electronic spectrum of [Fe(CO)5] that was published by Daniel and co-workers is used as a reference for the present study.72 It is assumed that these calculations cover all necessary static and the most significant part of the dynamic electron correlation present in [Fe(CO)5]. Table 3 presents the calculated excitation energies obtained from the MR-CCI, ROCIS, DFT/ ROCIS, and TD-DFT calculations. In the second column, the dominant electronic excitation of the observed excited states is listed under “dominant electronic transition”. Pure ROCIS,

Obtained from E = 0.769 × E(ROCIS) + 0.150. bObtained from E = 0.924 × E(DFT/ROCIS) + 0.489. cObtained from E = 1.076 × E(TDDFT) − 0.413.

a

the electronic spectra of the two transition metal carbonyl complexes feature d−d and metal-to-ligand charge transfer (MLCT) transitions. Calculation of the electronic excitation spectrum for the two complexes with ROCIS and DFT/ROCIS 3077

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

Table 3. Calculated Excitation Energies to the Low-Lying States of [Fe(CO)5] in Wavenumbers (cm−1) transition A1′ → a1E′ 1 A1′ → a1E″ 1 A′1 → 1A″1 1 A1′ → b1E″ 1 A1′ → 1A2″ 1 A′1 → 1A′2 1 A′1 → b1E′ 1 A1′ → 1A1′ 1

dominant electronic excitation MRCCI ROCIS 3dσ 3dπ 3dσ 3dσ 3dσ 3dσ 3dσ 3dσ

→ → → → → → → →

3dz2 3dz2 π*CO π*CO π*CO π*CO π*CO π*CO

27680 36760 33280 35540 37030 38050 39330 45030

31750 30900 30350 33240 37480 30000 33590 45383

DFT/ ROCIS

TD-DFT

30460 39470 30170 31160 32070 32910 37690 39050

31920 39130 30070 32180 32830 34490 35050 39910

Table 4. Calculated Excitation Energies to the Low Lying States of [Cr(CO)6] in Wavenumbers (cm1‑) transition A1g → a1Eu A1g → a1A2u 1 A1g → a1T2u 1 A1g → b1Eu 1 A1g → a1A1u 1 A1g → b1T2u 1 A1g → a1T1u 1

1

A1g → a1Eg A1g → a1T1g 1 A1g → b1T1g 1 A1g → a1T2g 1

1

DFT/ROCIS, and TD-DFT overestimate the excitation energy of the 3dσ → 3dz2 transition by 4000−4300 cm−1. Both orbitals involved in this transition are σ-antibonding, which indicates that the description of the σ-interaction between the metal and the ligands is similar for HF and DFT. In the case of the 3dπ → 3dz2 transition, the results obtained with the Hartree−Fock based ROCIS and the two DFT-based methods differ significantly. While ROCIS gives an excitation energy that is by 6000 cm−1 too low, DFT/ROCIS and TD-DFT overestimate the excitation energy by 2700 cm−1 and 2400 cm−1, respectively. This observation may be attributed to the π-bonding properties of the respective underlying reference methods. HF is known to underestimate the metal to ligand backbonding, whereas DFT tends to produce slightly too strong backbonds. Four of the six calculated MLCT transitions of [Fe(CO)5] are predicted too low in energy by ROCIS with the energetic deviation ranging from 2300 cm−1 to 8000 cm−1. The remaining two MLCT transitions are reproduced quite accurately with a deviation of 500 cm−1 and 300 cm−1, respectively. Similarly, DFT/ROCIS systematically underestimates the transition energy of all six MLCT transitions by 1600 cm−1 to 6000 cm−1. The TD-DFT excitation energies for MLCT transitions follow the numbers obtained with DFT/ ROCIS closely. The electronic ground state of [Cr(CO)6)] corresponds to a (t2g)6(eg)0 configuration. Hence, the electronic excitation spectrum features two 3-fold degenerate d−d states (b1T1g and a1T2g) alongside a larger number of MLCT states.70 For the same reason as discussed above, the reference values for electronic excitation energies of [Cr(CO)6] are taken from a correlated ab initio multireference study. Here, the CASPT2 calculations conducted by Pierloot and co-workers serve as a reference.73 Table 4 presents the CASPT2 excitation energies compared to the results obtained with ROCIS, DFT/ROCIS, and TD-DFT. For some of the MLCT transitions, an energy range is given for the CASPT2 reference energy. In these case, two calculations with different active spaces were performed. A detailed description of the various active spaces applied in the calculations can be found in ref 73. Pure ROCIS slightly underestimates the excitation energies of the two d−d excited states by about 2000 cm−1. Most of the MLCT states on the other hand are predicted 15000 cm−1 to 20000 cm−1 too high in energy. The only exceptions are energies of the states a1Eg and a1T1g that arise from the transition of an electron into the ligand based t2g orbitals, which are strongly involved in the metal-to-ligand backbonding. The deviations for these two transitions amount to −1500 cm−1 (a1Ea) and 3600 cm−1 (1T1g). Generally, the DFT/ROCIS excitation energies for [Cr(CO)6] are closer to the reference values than the pure ROCIS energies. DFT/ROCIS systematically overestimates the

dominant electronic excitation t2g → t1u t2g → t1u t2g → t1u t2g → t2u t2g → t2u t2g → t2u t2g → t1u t2g → t2u t2g → t2g t2g → t2g t2g → eg t2g → eg

CASPT2

DFT/ ROCIS ROCIS TD-DFT

27500−28960 28870 28710−29840 32020−32670 33070−33470 34840−35730 33150−36620

44200 45220 45440 52560 52950 56420 51403

33810 34420 32780 37100 39060 38060 35890

33580 34490 34110 38000 37860 38860 36510

36940 38880 39120 40970

38440 35262 37180 38940

34580 39070 42260 44790

41400 38450 40930 46750

excitation energies by about 2000 cm−1 to 5600 cm−1. Again, the transitions to the a1Eg and a1T1g states constitute the only exceptions with deviations of 100 cm−1 and −2500 cm−1, respectively. As in the case of [Fe(CO)5], most TD-DFT transition energies are close to the values obtained with DFT/ ROCIS. Partially large deviations of the calculated ROCIS excitation energies from the reference values reflect the necessity to include large amounts of electron correlation to obtain accurate results for the electronic spectra of transition metal complexes with complicated bonding patterns. The effect appears to be particularly pronounced for MLCT transitions in [Cr(CO)6] for which DFT/ROCIS and TD-DFT consequently yield significantly more accurate transition energies. However, despite being superior to the energies obtained with ROCIS, the DFT/ ROCIS and TD-DFT transition energies still deviate from the reference values by a few thousand wavenumbers. The deviation of the DFT/ROCIS energies from the reference values for all MLCT transitions in both complexes appears to be highly systematic. A special handling of the exchange parameter c2 in eq 1.10 that explicitly takes into account the magnitude of the respective exchange integrals might offer a route to more accurate and reliable charge transfer energies. Work along these lines is in progress. [Cr(acac)3]. Tris(acetylacetonato)chromium(III) ([Cr(acac)3]) is a spectroscopically well-studied model complex. The three acac ligands are arranged in such a way that the CrO6 core exhibits almost perfect octahedral symmetry (Figure 2). Assuming Oh symmetry, the electronic ground state transforms as 4A2g. The electronic absorption spectrum of [Cr(acac)3] features two bands in the visible region at ∼18100 cm−1 and ∼22700 cm−1. They have been assigned to the 4T2g and 4T1g excited states that follow from the (t2g)2(eg)1 excited configuration.74 However, the symmetry imposed on the chromium center by the π-electron system is D3 rather than Oh. This symmetry reduction affects the electronic structure and hence the electronic absorption spectrum: the 4T2g excited state splits up into states with 4A1 and 4E symmetry at 17700 cm−1 and 18500 cm−1, respectively. Using polarized light, transitions to the two states could be separated by their polarization parallel (4A1) and perpendicular (4E) to the molecular axis.74 With calculated bands at 22200 cm−1 and 35200 cm−1, the ROCIS method significantly overestimates the d−d excitation energies. This finding underlines the importance of dynamic correlation for a realistic description of excitation energies in 3078

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

were introduced and optimized. One of those structures, denoted the “hydrophobic model” (Figure 3), was chosen to

Figure 2. Ball and stick visualization of [Cr(acac)3]. The CrO6 core is almost perfectly Oh symmetric.

Figure 3. The “hydrophobic model” of the active site of [NiFe] hydrogenases76 including the inner core with the two metal centers, a CO, a hydride, and two CN− ligands as well as four cysteine residues. In addition, the model contains the Ala477, Pro478, Leu482, Val500, and Pro501 amino acid residues.

transition metal complexes. However, the splitting of the 4T2g state is qualitatively reproduced correctly. In the ROCIS calculation, the 4A1 is located at 22100 cm−1 while the 4E state comes at 22400 cm−1. Only the size of the splitting is underestimated. An inspection of the polarization of the transitions reveals that the character of the excited states is also reproduced correctly. The 4A2→4A1 transition is polarized along the molecular axis, whereas the 4A2→4E transition is polarized perpendicular to the molecular axis. Expectedly, the implicit description of dynamic correlation within the DFT/ ROCIS approach leads to improved excitation energies. The 4 T2g and 4T1g states are calculated around 15000 cm−1 and 17000 cm−1. As above, the splitting of the 4T2g state into a 4A1 and 4E state and the character of the two resulting excited states is predicted qualitatively correctly. The transition to the 4A1 state is polarized along the molecular z-axis and comes at 14300 cm−1, while the transition to the 4E state is located at 15300 cm−1 and polarized perpendicular to the molecular z-axis. In contrast to the above made observation for the pure ROCIS calculations, the splitting is slightly overestimated. TD-DFT predicts the 4A1 and 4E states at 21200 cm−1 and 21700 cm−1, respectively. It furthermore correctly reproduces the polarization of the two states but underestimates their energetic separation. For the 4T1g state, an excitation energy of 24400 cm−1 is obtained with TD-DFT. Hence, as was observed for the carbonyl complexes above, the TD-DFT and DFT/ROCIS calculations yield results of comparable accuracy. Generally, it is a challenging task to converge ROHF and ROKS calculations for transition metal complexes. Therefore, in our implementation it is also possible to use quasi-restricted orbitals (qro’s)75 of unrestricted Hartree−Fock and Kohn− Sham calculations for the reference function. Usually, the usage of qro’s has only a minor effect on the calculated transition energies and transition moments. In fact, as discussed elsewhere, the QRO and ROHF/ROKS determinants are nearly identical.52 Ni−C State of [NiFe] Hydrogenase and [MoHIPTN3N]NH: Timings. Hydrogenase enzymes catalyze the reversible formation of H2 in many microorganisms.77,78 Recently, a computational study of the Ni−C state of [NiFe] hydrogenases was performed by Kampa et al.76 In the course of the investigation, several model structures of the active site including several amino acid residues of the protein backbone

test the efficiency of the presented DFT/ROCIS method. The two metal centers with their first coordination shell, including a CO, a hydride and two CN− ligands as well as four cysteine residues, form the core of the model structure. As depicted in Figure 3 the model structure furthermore contains the Ala477, Pro478, Leu482, Val500 and Pro501 amino acid residues. In total, the hydrophobic model comprises 98 atoms with a total spin of S = 1/2 and a 2-fold negative charge (Figure 3). The def2-TZVP(-f) basis set47 was used for the Ni- and Fe-centers, the hydride atom, the two CO ligands, as well as all N and S atoms. For all other atoms, the def2-SVP basis set47 was employed resulting in a total number of 1011 basis functions. All calculations were carried out on a machine with 4 octacore Intel Xeon CPU E7-8837 with 2.67 GHz and 256 GB memory of which 4 GB per CPU core were used during the presented calculations. The RI approximation was applied during the integral transformation and the Fock matrix formation.53−56 Calculation of the ground and three excited states with DFT/ ROCIS on a single CPU core is achieved in 107751 s. The most time-consuming step of the calculation is the solution of the eigenvalue problem with the Davidson algorithm (95515 s), which needs 19 iterative cycles to reach convergence. Formation of the Fock matrices (4004 s) and the four-index integrals (8127 s) contributes only sparsely to the total CPU time. With increasing basis set size, the formation of the Fock matrices becomes more expensive relative to the other steps. Application of the, e.g., chainof-spheres approximation to the exchange (COSX)79 terms in the Fock and Kohn−Sham matrices should effectively reduce the computational cost for this step at almost no loss of accuracy. This logical extension will be pursued in the future. An equivalent TD-DFT calculation that uses a AO-based direct algorithm and invokes the RI approximation for Coulomb and exchange integrals80 requires only little more than half the time of the DFT/ROCIS calculation (59075 s). This speedup is closely related to the reduced number of iterative cycles needed to reach convergence (12), which we attribute to the generally smaller off-diagonal entries in the TD-DFT matrix as opposed to the off-diagonal entries in the DFT/ROCIS matrix. With increasing basis set size, the TD-DFT calculation will 3079

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

employed throughout utilizing the QZV/J auxiliary basis set.87 Figure 5 presents the computational speedup of ROCIS and

become more favorable in terms of computational costs due to its more advantageous scaling behavior. The effect will be further amplified if the COSX79 approximation to the exchange terms is invoked. These results demonstrate that the generation of a rigorously spin-adapted excitation space comes at the cost of more complicated equations than in the case of TD-DFT and hence a less efficient algorithm for their solution. However, albeit being less efficient than TD-DFT, DFT/ROCIS is applicable to medium-sized molecular systems with about 100 atoms and an open-shell ground state. The efficient parallelization of the program (vide inf ra) further increases its scope with respect to the system size and the number of calculated excited states. The effectiveness of the parallelization was tested for a model structure of the molybdenum trisamidoamine complex [Mo] (=(3,5-(2,4,6-i-Pr3C6H2)2C6H3NCH2CH2N)Mo) shown in Chart 3. [Mo] is the active species in the catalytic reduction Chart 3. Lewis Structure of [Mo]NH

Figure 5. Speed up of ROCIS calculations of the reduced model of [Mo]NH using the QZVP basis set with increasing number of CPU cores.

DFT/ROCIS versus the number of used CPU cores. Until a number of 16 CPU cores, the speedup increases strictly monotonically, reaching its maximum values of 12.7 and 11.8 for DFT/ROCIS and ROCIS, respectively. The increase is almost linearly until a number of 4 CPU cores. Above this number, the synchronization of and communication between the parallel cores compromise the speedup per CPU core.88 In general, DFT/ ROCIS scales slightly better than ROCIS due to the efficiently parallelized formation of the Kohn−Sham matrix. By construction, a DFT/ROCIS calculation always takes slightly longer than the equivalent ROCIS calculation.

of N2 to ammonia by protons and electrons at room temperature.81,82 Recently a combined experimental and theoretical study of the EPR spectroscopic properties of three derivatives of the [Mo] complex has been published.83 In the course of this study a model structure of this complex has been introduced. Instead of the large HIPT groups, the model structure features smaller ethyl residues (Figure 4). The ethyl groups serve the

■

DISCUSSION In this work we have presented an efficient implementation of a spin-adapted ROCIS method for the calculation of electronically excited states of open-shell molecules. The design of the multielectron spin-eigenfunction expansion basis for the excited states has received special attention in the treatment. Inspired by earlier works of Grimme and co-workers, a combination of the ROCIS expansion space with a DFT reference function in the framework of a parametrized CI scheme has been developed. The application of the ROCIS method to the calculation of electronically excited states of a series of organic radicals demonstrated some key properties of the method. All transition energies are calculated too high in energy for two reasons. First, the method calculates vertical transition energies, whereas experimentally only 0→0 transition energies are accessible. Moreover, the ROCIS method neglects electron correlation effects in the ground and all excited states. A linear correlation has shown that the introduced errors are highly systematic and can be accounted for by downscaling the calculated excitation energies by ∼20%. When applied to transition metal complexes, the ROCIS method shows, as expected, more significant deviations from experiment. For example, although the number, degeneracy, and character of the expected excited states for [Cr(acac)3] are reproduced correctly, the calculated transition energies are shifted by up to 12500 cm−1. In the case of iron and chromium carbonyl complexes, the ROCIS excitation

Figure 4. Model structure for [Mo]NH featuring ethyl residues instead of HIPT groups.

dual purpose of mimicking the electron donating effect of the HIPT groups and saving computational time. Here, serial and parallel ROCIS and DFT/ROCIS calculations of the ground and four excited states of [Mo]NH were performed. Scalar relativistic effects were taken into account within the zeroth order regular approximation (ZORA).84−86 Consequently, the scalarrelativistically recontracted QZVP basis set87 was used resulting in a total number of 1890 basis functions. The RI approximation was 3080

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

diagonal elements of the DFT/ROCIS matrix, equations for ROCIS and DFT/ROCIS density matrices, reference contributions to the electronic ground state of open-shell molecules, and DFT/ROCIS results with BhLYP. This material is available free of charge via the Internet at http://pubs.acs.org.

energies deviate from highly correlated multireference ab initio data by up to 20000 cm−1. The most important reason for this erroneous behavior is the lack of dynamic correlation in the ROCIS treatment, which has an even more pronounced and complex effect on the electronic structure of transition metal complexes than in the case of organic compounds. Dynamic correlation effects are partially taken into account by the usage of DFT orbitals resulting in the DFT/ROCIS method. In order to compensate the generally too large two-electron integrals of DFT orbitals occurring in the CI matrix, DFT/ROCIS utilizes a parametrizing scheme following the ideas of Grimme’s DFT/SCI and DFT/MRCI methods. It features the reduction of Coulomb and exchange integrals by scaling parameters that are purely empirical but also applicable to a wide range of molecules including organic compounds and transition metal complexes. These parameters have been determined by optimizing the calculated transition metal L-edges for a test set of first-row transition metal complexes. As expected, DFT/ROCIS yields improved excitation energies compared to pure ROCIS. The DFT/ROCIS transition energies of organic radicals are significantly closer to the experimental values. However, in contrast to the ROCIS transition energies, they are systematically too low, which can be accounted for by an upscaling of the calculated values. Furthermore the predicted d−d transition energies of [Cr(acac)3] agree well with experiment. As for the pure ROCIS calculations, the splitting of the lowest 3T2g state into polarized components is reproduced correctly. The most prominent improvement of calculated transition energies was observed for [Fe(CO)5] and [Cr(CO)6]. While pure ROCIS produces large errors for the transition energies of the observed charge transfer transitions the DFT/ROCIS energies deviate only on the order 0−6000 cm−1 from the reference values in a very systematic fashion. A special handling of the exchange parameter in the present parametrization scheme might help to further improve the calculated excitation energies for charge transfer transitions. Our results demonstrate that the usage of DFT orbitals in combination with the presented parametrization scheme significantly improves the obtained results. Due to its inherent single reference character and the lack of higher excitations than single excitations, it does not reach the same accuracy as DFT/MRCI or highly correlated ab initio methods such as CASPT2 or MRCI, which cover a much larger expansion space. DFT/ ROCIS on the other hand has the advantage of being applicable to large systems, which at the moment cannot be treated by these highly accurate methods. In terms of accuracy, the DFT/ ROCIS results compare well with the TD-DFT results for the presented cases. Test calculations on large models of the active site of [NiFe] hydrogenases demonstrate that DFT/ROCIS is suitable for application to large systems with about 100 atoms and more. Thus DFT/ROCIS constitutes a method of moderate accuracy for large transition metal complexes and organic compounds. It has to be noted that TD-DFT achieves a comparable accuracy at a lower cost and better scaling with increasing system size. However, DFT/ROCIS has the advantage that it keeps spin symmetry in the excited CSF manifold, making it a suitable basis for the treatment of spin-related interactions such as spin−orbit coupling on a multielectron level. Extensions of the basic method in this direction are readily developed and will be reported in due course.

■

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Dreuw, A.; Head-Gordon, M. Single-Reference Ab Initio Methods for the Calculation of Excited States of Large Molecules. Chem. Rev. 2005, 105, 4009−4037. (2) Roos, B. O.; Andersson, K.; Fülscher, M. P.; Malmqvist, P.; Serrano-Andrés, L.; Pierloot, K.; Merchán, M. Multiconfigurational Perturbation Theory: Applications in Electronic Spectroscopy. Adv. Chem. Phys. 1996, 219−331. (3) Neese, F.; Petrenko, T.; Ganyushin, D.; Olbrich, G. Advanced Aspects of Ab Initio Theoretical Optical Spectroscopy of Transition Metal Complexes: Multiplets, Spin-Orbit Coupling and Resonance Raman Intensities. Coord. Chem Rev 2007, 251, 288−327. (4) Buenker, R. J.; Peyerimhoff, S. D. Individualized Configuration Selection in CI Calculations with Subsequent Energy Extrapolation. Theor. Chim. Acta 1974, 35, 33−58. (5) Miralles, J.; Daudey, J.-P.; Caballol, R. Variatonal Calculation of Small Energy Differences. The Singlet−Triplet Gap in [Cu2Cl6]2−. Chem. Phys. Lett. 1992, 198, 555−562. (6) Miralles, J.; Castell, O.; Caballol, R.; Malrieu, J.-P. Specific CI Calculation of Energy Differences: Transition Energies and Bond Energies. Chem. Phys. 1993, 172, 33−43. (7) Neese, F. A Spectroscopy Oriented Configuration Interaction Procedure. J. Chem. Phys. 2003, 119, 9428−9443. (8) Andersson, K.; Malmqvist, P. A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. Second-Order Perturbation Theory with a CASSCF Reference Function. J. Phys. Chem. 1990, 94, 5483−5488. (9) Andersson, K.; Malmqvist, P.-A.; Roos, B. O. Second-Order Perturbation Theory with a Complete Active Space Self-Consistent Field Reference Function. J. Chem. Phys. 1992, 96, 1218−1226. (10) Andersson, K.; Roos, B. O. In Modern Electronic Structure Theory; Yarkony, D., Ed.; World Scientific: Singapore, 1995; p 55. (11) Head-Gordon, M.; Rico, R. J.; Oumi, M.; Lee, T. J. A Doubles Correction to Electronic Excited States from Configuration Interaction in the Space of Single Substitutions. Chem. Phys. Lett. 1994, 219, 21− 29. (12) Head-Gordon, M.; Maurice, D.; Oumi, M. A Perturbative Correction to Restricted Open Shell Configuration Interaction with Single Substitutions for Excited States of Radicals. Chem. Phys. Lett. 1995, 246, 114−121. (13) Sekino, H.; Bartlett, R. J. A Linear Response, Coupled-Cluster Theory for Excitation Energy. Int. J. Quantum Chem. 1984, 26, 255− 265. (14) Koch, H.; Jorgensen, P. Coupled Cluster Response Functions. J. Chem. Phys. 1990, 93, 3333−3344. (15) Christiansen, O.; Koch, H.; Jørgensen, P. The Second-Order Approximate Coupled Cluster Singles and Doubles Model CC2. Chem. Phys. Lett. 1995, 243, 409−418. (16) Hattig, C.; Weigend, F. CC2 Excitation Energy Calculations on Large Molecules Using the Resolution of the Identity Approximation. J. Chem. Phys. 2000, 113, 5154−5161. (17) Hattig, C.; Kohn, A.; Hald, K. First-Order Properties for Triplet Excited States in the Approximated Coupled Cluster Model CC2 Using an Explicitly Spin Coupled Basis. J. Chem. Phys. 2002, 116, 5401−5410.

ASSOCIATED CONTENT

S Supporting Information *

Calculated transition energies of organic closed shell compounds with DFT/ROCIS, DFT/MRCI, and TD-DFT, Parameterized 3081

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

(18) Geertsen, J.; Rittby, M.; Bartlett, R. J. The Equation-of-Motion Coupled-Cluster Method: Excitation Energies of Be and CO. Chem. Phys. Lett. 1989, 164, 57−62. (19) Stanton, J. F.; Bartlett, R. J. The Equation of Motion CoupledCluster Method. A Systematic Biorthogonal Approach to Molecular Excitation Energies, Transition Probabilities, and Excited State Properties. J. Chem. Phys. 1993, 98, 7029−7039. (20) Runge, E.; Gross, E. K. U. Density-Functional Theory for TimeDependent Systems. Phys. Rev. Lett. 1984, 52, 997−1000. (21) Neese, F. Prediction of Molecular Properties and Molecular Spectroscopy with Density Functional Theory: From Fundamental Theory to Exchange-Coupling. Coord. Chem. Rev. 2009, 253, 526−563. (22) Hirata, S.; Head-Gordon, M. Time-Dependent Density Functional Theory within the Tamm−Dancoff Approximation. Chem. Phys. Lett. 1999, 314, 291−299. (23) Schirmer, J.; Dreuw, A. Critique of the Foundations of TimeDependent Density-Functional Theory. Phys. Rev. A 2007, 75, 022513. (24) Grimme, S.; Parac, M. Substantial Errors from Time Dependent Density Functional Theory for the Calculation of Excited States of Large Systems. ChemPhysChem 2003, 4, 292−295. (25) Tozer, D. J.; Amos, R. D.; Handy, N. C.; Roos, B. O.; SerranoAndres, L. Does Density Functional Theory Contribute to the Understanding of Excited States of Unsaturated Organic Compounds? Mol. Phys. 1999, 97, 859−868. (26) Grimme, S. Density Functional Calculations with Configuration Interaction for the Excited States of Molecules. Chem. Phys. Lett. 1996, 259, 128−137. (27) Grimme, S.; Waletzke, M. A Combination of Kohn−Sham Density Functional Theory and Multi-Reference Configuration Interaction Methods. J. Chem. Phys. 1999, 111, 5645−5655. (28) Silva-Junior, M. R.; Schreiber, M.; Sauer, S. P. A.; Thiel, W. Benchmarks for Electronically Excited States: Time-Dependent Density Functional Theory and Density Functional Theory based Multireference Configuration Interaction. J. Chem. Phys. 2008, 129, 104103. (29) Edwards, W. D.; Weiner, B.; Zerner, M. C. On the Low-Lying States and Electronic Spectroscopy of Iron(II) Porphine. J. Am. Chem. Soc. 1986, 108, 2196−2204. (30) Anderson, W. P.; Edwards, W. D.; Zerner, M. C. Calculated Spectra of Hydrated Ions of the First Transition-Metal Series. Inorg. Chem. 1986, 25, 2728−2732. (31) Maurice, D.; Head-Gordon, M. Configuration Interaction with Single Substitutions for Excited States of Open-Shell Molecules. Int. J. Quantum Chem. 1995, 56, 361−370. (32) Edwards, W. D.; Zerner, M. C. A Generalized Restricted OpenShell Fock Operator. Theor. Chim. Acta 1987, 72, 347−361. (33) Ridley, J.; Zerner, M. An Intermediate Neglect of Differential Overlap Technique for Spectroscopy: Pyrrole and the Azines. Theor. Chim. Acta 1973, 32, 111−134. (34) Bacon, A. D.; Zerner, M. C. An Intermediate Neglect of Differential Overlap Theory for Transition Metal Complexes: Fe, Co and Cu Chlorides. Theor. Chim. Acta 1979, 53, 21−54. (35) Zerner, M. C.; Loew, G. H.; Kirchner, R. F.; MuellerWesterhoff, U. T. An Intermediate Neglect of Differential Overlap Technique for Spectroscopy of Transition-Metal Complexes. Ferrocene. J. Am. Chem. Soc. 1980, 102, 589−599. (36) Ake, R. L.; Gouterman, M. Porphyrins XIV. Theory for the Luminescent State in VO, Co, Cu Complexes. Theor. Chim. Acta 1969, 15, 20−42. (37) Petrenko, T.; Ray, K.; Wieghardt, K. E.; Neese, F. Vibrational Markers for the Open-Shell Character of Transition Metal BisDithiolenes: An Infrared, Resonance Raman, and Quantum Chemical Study. J. Am. Chem. Soc. 2006, 128, 4422−4436. (38) Cory, M. G.; Zerner, M. C. Metal−Ligand Exchange Coupling in Transition-Metal Complexes. Chem. Rev. 1991, 91, 813−822. (39) Maurice, D.; Head-Gordon, M. On the Nature of Electronic Transitions in Radicals: An Extended Single Excitation Configuration Interaction Method. J. Phys. Chem. 1996, 100, 6131−6137.

(40) McLean, A. D.; Liu, B. Classification of Configurations and the Determination of Interacting and Noninteracting Spaces in Configuration Interaction. J. Chem. Phys. 1973, 58, 1066−1078. (41) Neese, F. The ORCA Program System. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 73−78. (42) Hsiao, Y. W.; Zerner, M. C. Calculating ESR G Tensors of Doublet Radicals by the Semiempirical INDO/S Method. Int. J. Quantum Chem. 1999, 75, 577−584. (43) McWeeny, R. Methods of Molecular Quantum Mechanics, 2 ed.; Academic Press: London, 1992. (44) Roemelt, M.; Maganas, D.; DeBeer, S.; Neese, F. submitted for publication, 2013. (45) Becke, A. D. Density-Functional Thermochemistry 0.3. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (46) Lee, C. T.; Yang, W. T.; Parr, R. G. Development of the Colle− Salvetti Correlation-Energy Formula into a Functional of the ElectronDensity. Phys. Rev. B 1988, 37, 785−789. (47) 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. (48) Becke, A. D. A New Mixing of Hartree−Fock and Local Density-Functional Theories. J. Chem. Phys. 1993, 98, 1372. (49) Davidson, E. R. The Iterative Calculation of a Few of the Lowest Eigenvalues and Corresponding Eigenvectors of Large Real-Symmetric Matrices. J. Comput. Phys. 1975, 17, 87−94. (50) Crouzeix, M.; Philippe, B.; Sadkane, M. The Davidson Method. SIAM J. Sci. Comput. 1994, 15, 62−76. (51) Lowdin, P. O. On The Non-orthogonality Problem Connected with the Use of Atomic Wave Functions in the Theory of Molecules and Crystals. J. Chem. Phys. 1950, 18, 365−375. (52) Sinnecker, S.; Neese, F. Spin−Spin Contributions to the ZeroField Splitting Tensor in Organic Triplets, Carbenes and Biradicals − A Density Functional and Ab Initio Study. J. Phys. Chem. A 2006, 110, 12267−12275. (53) Baerends, E. J.; Ellis, D. E.; Ros, P. Self Consistent Molecular Hartree−Fock−Slater Calculations − I. The Computational Procedure. Chem. Phys. 1973, 2, 41−51. (54) Whitten, J. L. Coulombic Potential-Energy Integrals and Approximations. J. Chem. Phys. 1973, 58, 4496−4501. (55) Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. Some Approximations in Applications of X-Alpha Theory. J. Chem. Phys. 1979, 71, 3396−3402. (56) Vahtras, O.; Almlö f , J.; Feyereisen, M. W. Integral Approximations for LCAO-SCF Calculations. Chem. Phys. Lett. 1993, 213, 514−518. (57) Eichkorn, K.; Treutler, O.; Ohm, H.; Haser, M.; Ahlrichs, R. Auxiliary Basis-Sets to Approximate Coulomb Potentials. Chem. Phys. Lett. 1995, 240, 283−289. (58) Craig, D. P.; Thirunamachandran, T. Molecular Quantum Electrodynamics; Academic Press: London, 1984. (59) Griffith, J. S. The Theory of Transition Metal Ions; Cambridge University Press: Cambridge, U.K., 1964. (60) DeBeer, S.; Petrenko, T.; Neese, F. Time-Dependent Density Functional Calculations of Ligand K-Edge X-ray Absorption Spectra. Inorg. Chim. Acta 2008, 361, 965−972. (61) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098. (62) Perdew, J. P. Density-Functional Approximation for the Correlation Energy of the Inhomogeneous Electron Gas. Phys. Rev. B 1986, 33, 8822. (63) Weigend, F. Accurate Coulomb-Fitting Basis Sets for H to Rn. Phys. Chem. Chem. Phys. 2006, 8, 1057−1065. (64) Dierksen, M.; Grimme, S. The Vibronic Structure of Electronic Absorption Spectra of Large Molecules: A Time-Dependent Density Functional Study on the Influence of “Exact” Hartree−Fock Exchange. J. Phys. Chem. A 2004, 108, 10225−10237. 3082

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083

The Journal of Physical Chemistry A

Article

(65) Dewar, J. S. A Review of the π-Complex Theroy. Bull. Soc. Chim. Fr. 1951, 18, C71−C79. (66) Chatt, J.; Duncanson, L. A. 586. Olefin Co-Ordination Compounds. Part III. Infra-Red Spectra and Structure: Attempted Preparation of Acetylene Complexes. J. Chem. Soc. (Resumed) 1953, 2939−2947. (67) Elian, M.; Hoffmann, R. Bonding Capabilities of Transition Metal Carbonyl Fragments. Inorg. Chem. 1975, 14, 1058−1076. (68) Persson, B. J.; Roos, B. O.; Pierloot, K. A Theoretical Study of the Chemical Bonding in M(CO) (M = Cr, Fe, and Ni). J. Chem. Phys. 1994, 101, 6810. (69) Ehlers, A. W.; Dapprich, S.; Vyboishchikov, S. F.; Frenking, G. Structure and Bonding of the Transition-Metal Carbonyl Complexes M(CO)5L (M = Cr, Mo, W) and M(CO)3L (M = Ni, Pd, Pt; L = CO, SiO, CS, N2, NO+, CN−, NC−, HCCH, CCH2, CH2, CF2, H2)1. Organometallics 1996, 15, 105−117. (70) Kotzian, M.; Roesch, N.; Schroeder, H.; Zerner, M. C. Optical Spectra of Transition-Metal Carbonyls: Chromium Hexacarbonyl, Iron Pentacarbonyl, and Nickel Tetracarbonyl. J. Am. Chem. Soc. 1989, 111, 7687−7696. (71) Werner, H. J.; Knowles, P. J. An Efficient Internally Contracted Multiconfiguration−Reference Configuration Interaction Method. J. Chem. Phys. 1988, 89, 5803. (72) Rubner, O.; Engel, V.; Hachey, M.; Daniel, C. A CASSCF/MRCCI Study of the Excited States of Fe (CO) 5. Chem. Phys. Lett. 1999, 302, 489−494. (73) Pierloot, K.; Tsokos, E.; Vanquickenborne, L. Optical Spectra of Ni(CO)4 and Cr(CO)6 Revisited. J. Phys. Chem. 1996, 100, 16545− 16550. (74) Atanasov, M. A.; Schönherr, T.; Schmidtke, H. H. The Role of the π-Bonding Network for Trigonal Level Splittings of Tris-Bidentate [Cr(acac)3] and [Cr(ox)3]3−. Theor. Chim. Acta 1987, 71, 59−73. (75) Ray, K.; Begum, A.; Weyhermüller, T.; Piligkos, S.; van Slageren, J.; Neese, F.; Wieghardt, K. The Electronic Structure of the Isoelectronic, Square Planar Complexes [Fe II (L) 2 ] 2− and [CoIII(LBu)2]− (L2− and (LBu)2− = Benzene-1,2-dithiolates): An Experimental and Density Functional Theoretical Study. J. Am. Chem. Soc. 2005, 127, 4403−4415. (76) Kampa, M.; Lubitz, W.; Gastel, M.; Neese, F. Computational Study of the Electronic Structure and Magnetic Properties of the Ni− C State in [NiFe] Hydrogenases Including the Second Coordination Sphere. J. Biol. Inorg. Chem. 2012, 17, 1269−1281. (77) Adams, M. W. W. The Structure and Mechanism of IronHydrogenases. Biochim. Biophys. Acta, Bioenerg. 1990, 1020, 115−145. (78) Vignais, P. M.; Billoud, B.; Meyer, J. Classification and Phylogeny of Hydrogenases1. FEMS Microbiol. Rev. 2001, 25, 455− 501. (79) Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Efficient, Approximate and Parallel Hartree−Fock and Hybrid DFT Calculations. A “Chain-of-Spheres” Algorithm for the Hartree−Fock Exchange. Chem. Phys. 2009, 356, 98−109. (80) Neese, F.; Olbrich, G. Efficient Use of the Resolution of the Identity Approximation in Time-Dependent Density Functional Calculations with Hybrid Density Functionals. Chem. Phys. Lett. 2002, 362, 170−178. (81) Yandulov, D. V.; Schrock, R. R. Reduction of Dinitrogen to Ammonia at a Well-Protected Reaction Site in a Molybdenum Triamidoamine Complex. J. Am. Chem. Soc. 2002, 124, 6252−6253. (82) Yandulov, D. V.; Schrock, R. R. Catalytic Reduction of Dinitrogen to Ammonia at a Single Molybdenum Center. Science 2003, 301, 76−78. (83) McNaughton, R. L.; Roemelt, M.; Chin, J. M.; Schrock, R. R.; Neese, F.; Hoffman, B. M. Experimental and Theoretical EPR Study of Jahn−Teller-Active [HIPTN3N]MoL Complexes (L = N2, CO, NH3). J. Am. Chem. Soc. 2010, 132, 8645−8656. (84) Lenthe, E. v.; Baerends, E. J.; Snijders, J. G. Relativistic Regular Two-Component Hamiltonians. J. Chem. Phys. 1993, 99, 4597−4610. (85) Heully, J.; Lindgren, I.; Lindroth, E.; Lundqvist, S.; MartenssonPendrill, A. Diagonalisation of the Dirac Hamiltonian as a Basis for a

Relativistic Many-Body Procedure. J. Phys. B: At. Mol. Phy. 1986, 19, 2799−2815. (86) van Wüllen, C. Molecular Density Functional Calculations in the Regular Relativistic Approximation: Method, Application to Coinage Metal Diatomics, Hydrides, Fluorides and Chlorides, and Comparison with First-Order Relativistic Calculations. J. Chem. Phys. 1998, 109, 392. (87) Pantazis, D. A.; Chen, X. Y.; Landis, C. R.; Neese, F. AllElectron Scalar Relativistic Basis Sets For Third-Row Transition Metal Atoms. J. Chem. Theory Comput. 2008, 4, 908−919. (88) Amdahl, G. M. In Proceedings of the April 18−20, 1967, Spring Joint Computer Conference; ACM: New York, 1967; p 483−485.

3083

dx.doi.org/10.1021/jp3126126 | J. Phys. Chem. A 2013, 117, 3069−3083