Explicitly correlated internally contracted multireference coupled-cluster singles and doubles theory: ic-MRCCSD(F12∗)

Chemical Physics Letters 565 (2013) 122–127

Contents lists available at SciVerse ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

Explicitly correlated internally contracted multireference coupled-cluster singles and doubles theory: ic-MRCCSD(F12*) Wenlan Liu ⇑, Matthias Hanauer, Andreas Köhn Institut für Physikalische Chemie, Johannes Gutenberg-Universität Mainz, Duesbergweg 10-14, D-55128 Mainz, Germany

a r t i c l e

i n f o

Article history: Received 27 November 2012 In final form 28 December 2012 Available online 19 January 2013

a b s t r a c t An explicitly correlated ansatz employing Slater-type geminals and cusp conditions is developed for the internally contracted multireference coupled-cluster singles and doubles method. Only the most important geminal terms are retained in the spirit of earlier work for single-reference theory. Throughout all our test calculations, the new ic-MRCCSD(F12*) method improves the basis set convergence of many properties, e.g., spectroscopic constants or singlet–triplet splittings, with only little extra computational cost. If a perturbative correction for connected triples is included (the ic-MRCCSD(F12*)+(T) method), very accurate results can be obtained even with minimal active spaces. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction For the quantitative treatment of systems with a complex electronic structure, such as diradicals, transition metal compounds, or intermediates during bond breaking processes, the accurate treatment of both dynamic and static correlation effects is equally important [1,2]. The appropriate choices for such systems are multireference (MR) methods, e.g., complete active space based second order perturbation (CASPT2) theory [1,3], the multireference configuration interaction (MRCI) method [1,4], the multireference averaged coupled pair functional (MRACPF) [1,5], or multireference coupled-cluster (MRCC) theory [2,6–8]. For all these methods, however, the description of dynamic correlation is hampered by its slow basis set convergence, due to the use of orbital product expansions, which are not well-suited for the short-range part of the Coulomb correlation. As an alternative to using huge basis sets, geminal functions that explicitly depend on the interelectronic distance may be employed. For molecular systems, the most promising strategy is the F12 method [9–12], which applies additional Slater-type geminal functions for describing the cusp region [13]. This approach has been successfully combined with single-reference methods, e.g., second-order Møller-Plesset perturbation (MP2) theory, or coupled-cluster (CC) theory [10,13–17]. In particular for CC theory, the newest generation of methods, CCSD-F12b [18] and CCSD(F12*) [19], already gives near basis-set limit results for triple-f type basis sets while the computational cost increases by typically 10–50 percent only, compared to a run without geminal correlation factors [12].

⇑ Corresponding author. E-mail addresses: [email protected] (W. Liu), (M. Hanauer), [email protected] (A. Köhn).

[email protected]

0009-2614/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2012.12.052

Explicit electron correlation for multireference methods has been pioneered by Gdanitz for the MRACPF method [20,21] (using linear correlation factors, then known as R12 theory). F12 theory has more recently been applied to MR-MP2 [22], CASPT2 [23], MRCI [24], Mukherjee’s state-specific multireference coupled-cluster singles and doubles (Mk-MRCCSD) method [25], and to a perturbative variant thereof [26]. In this letter, we present an F12 ansatz for the internally contracted MRCCSD (ic-MRCCSD) method [7,27–30]. This method appears particularly promising due to its size-extensivity and its robustness with respect to the choice of active spaces. 2. Theoretical approach 2.1. ic-MRCC theory The ic-MRCC theory employs a complete active space self-consistent field (CASSCF) wave function jW0 i as its reference function, which is a linear combination of determinants jUl i,

jW0 i ¼

X cl jUl i:

ð1Þ

l

In analogy to single-reference coupled-cluster theory, the ic-MRCC wavefunction jWi is given by ^

jWi ¼ eT jW0 i:

ð2Þ

In the present formulation of the theory [28], we use a Hamiltonian normal ordered with respect to the core determinant j0i (the determinant of the CASSCF inactive space),

^ ¼ h0jHj0i ^ H þ

X 1 X pq rs ^qp þ ^ ; fqp a g a 4 pqrs rs pq pq

ð3Þ

123

W. Liu et al. / Chemical Physics Letters 565 (2013) 122–127

^1 þ K ^2 ^ ¼K K X 0 A X 1 0 AB ^0y ÞIA þ ^0y ÞIJAB ðk ÞIJ ðs ¼ ðk ÞI ðs 4 IA IJAB

Table 1 Index convention used in this work. p; q; r; . . . i; j; k; . . . u; v ; w; . . . a; b; c; . . . 0 a0 ; b ; c0 ; . . . 00 00 a ; b ; c00 ; . . . I; J; K; . . . A; B; C; . . .

Orbitals in finite basis Occupied orbitals Active orbitals Virtual orbitals in finite basis Virtual orbitals in formally complete basis Dto. restricted to complementary space Union of occupied and active orbitals (internal orbitals) Union of active and virtual orbitals in finite basis

ð8Þ

and

hW0 j ¼

X cl hUl j

ð9Þ

l

 contain the Lagrange multipliers ðk0 ÞAI ; ðk0 ÞAB IJ , and c l . 2.2. Merge with F12 theory

where the indices p; q; r; s; . . . label the orbitals in a finite basis. Our index convention is given in Table 1. fqp and g pq rs are matrix elements of the core Fock operator and the two-electron operator, respec^qp ¼ fa ^p g and ^q a tively, both are defined in Table 2; a qs q s ^pr ¼ fa ^a ^r a ^p g are normal-ordered strings of creation (a ^ a ^p ) and a ^ annihilation (aq ) operators. With these definitions, the cluster operator T^ is defined for ic-MRCCSD as

X 1 X IJ AB ^AI þ ^ ; T^ ¼ t IA a t a 4 IJAB AB IJ IA

In the following, we will generalize the F12 ansatz developed for single-reference theory, specifically CCSD(F12*), [19,31], to icMRCCSD theory. We introduce a geminal function

^ 12 f ðr 12 ÞjKLi juKL i ¼ Q

where f ðr12 Þ ¼  1c ecr12 is a Slater-type correlation factor with the interelectronic distance r12 and a length scale parameter c. The pro^ 12 ensures strong orthogonality to the space spanned by jector Q orbital product wavefunctions

ð4Þ

^uv ). Inwhere we implicitly exclude purely internal excitations (e.g., a stead, the expansion coefficients of jW0 i, cl , are obtained from

X ^ ^ T^ hUl jeT He jUm icm ¼ Ecl :

^ 1 Þð1  O ^ 2 Þ  v^ 1 v^ 2 ; ^ 12 ¼ ð1  O Q P

m

The amplitudes of T^ are determined from ^

^

fð^0 ÞAI ; . . . ;

ð6Þ

^ 12 ¼ 1 R 8

ð^0 ÞAB IJ ; . . .g

where sq 2 s s is an element of the linearly independent excitation operator manifold (see Ref. [28]). For the subsequent extension to F12 theory, it is useful to merge Eqs. (5) and (6) into a stationary energy functional. The corresponding Lagrange function for the ic-MRCCSD energy is then given by ^0

X

0 0 e KL0 0 a ^a b ; cIJKL R a b IJ

ð12Þ

IJKLa0 b0

e KL0 0 ¼ S KL ha0 b0 juKL i contains the special where the matrix element R ab symmetrizer

S KL ¼

^ ^ ^ T^ þT^ ^ eT^ 1 T^ 2 He ^ T^1 þT^2 jW0 i 1 2 L ¼ hW0 jeT 1 T 2 He jW0 i þ hW0 jK

 EðhW0 jW0 i  1Þ

ð11Þ P

^n ¼ ^ where O I jIðnÞihIðnÞj and v n ¼ a jaðnÞihaðnÞj are projectors on the internal and finite virtual orbital subspaces, respectively. We can then construct geminal operators as double excitation operators from the internal space (I; J; . . .) into the formally com0 plete virtual space (a0 ; b ; . . .),

ð5Þ

^ T jW0 i ¼ 0; ^0yq eT He hW0 js

ð10Þ

3 þ P~KL : 8

ð13Þ

~ KL permutes the spatial part of spin orbitals u and u . The coeffiP K L cients cIJKL ¼ dIK dJL  dIL dJK can be determined from the cusp conditions for both singlet and triplet pairs [13]. Eq. (12) hence simplifies to

ð7Þ

where

Table 2 ^ 12 is the exchange operator that permutes electrons (P ^ 12 jpð1Þqð2Þi ¼ jpð2Þqð1Þi) and P ^ ðpjqÞ is the Definition of the integrals and intermediates used in this work. P ^ ðpjqÞ Apq ¼ Apq  Aqp ). The special symmetrizer S pq is defined in the text. ðc0 Þu and ðc0 Þuv are the zeroth-order density antisymmetrization operator on a pair of spin–orbital indices (P v wx ij ij ij ð0Þ matrix elements of the initial CASSCF wavefunction, jW0 i. ^ qp ð1Þ ¼ hpð2Þj^r 1 jqð2Þi J 12 i P h^i ^ ^ ^i Fð1Þ ¼ hð1Þ þ i J i ð1Þ  Ki ð1Þ h i P ^ ^ e ^ J ð1Þ ^ J ð1Þ  K F ð1Þ ¼ hð1Þ þ IJ ðc0 ÞIJ J I I

^ W0 i E0 ¼ hW0 jHj ^ fpq ¼ hpjFð1Þjqi g rs pq

¼

hpqj^r1 12 jrsi

BIJKL

V Ip ¼ CIJab

P

^ ðajbÞ ¼P

P

c ~ IJ c00 F a Rc00 b 00

ð0Þ

P

00

~ IJ00 CIJau ¼ b00 F bu R b a   ^ ðajbÞ P ^ ðIjJÞ P 00 g Ic00 R e Jk 00 þ P 00 g Ic00 R e Jv 00 ðc0 Þu U IJab ¼ P kc ak c uv au v bc bc P  P 00 00 ib e jk ib e Jw 0 v ^ ðijjÞ 00 00 U ijua ¼ P kb g uk R ab00 þ b v w g uv R ab00 ðc Þw P P P 00 00 00 jk j v u i 1 ib ib 0 1 ib e wx 0 uv e 00 þ e 00 00 U a ¼ 2 jkb00 g jk R jb uv g ju R ab00 ðc Þv þ 4 b uv wx g uv R ab00 ðc Þwx ab

ð0Þ

^uv jW0 i ðc0 Þuv ¼ hW0 ja P q q F p ¼ fp þ uv g qpuv ðc0 Þuv

^ 12 ^f ðr12 ÞjIJi e IJ0 0 ¼ S IJ ha0 b0 jQ R ab   ^ ^ ^ 12 ^f ðr12 Þ e ^ 12 jIJi F ð2Þ ^f ðr 12 ÞQ F ð1Þ þ e ¼ S IJ S KL hKLjQ

^ 12 ^f ðr12 ÞjIJi X IJKL ¼ S IJ S KL hKLj^f ðr 12 ÞQ ^ ^ V IJpq ¼ S IJ hpqj^r1 12 Q12 f ðr 12 ÞjIJi w ~ uw ¼ Buw þ P ^ ðujwÞ ðP X uk F w  P X uy B vx vx k vx k y vx F y Þ P P P P lk u P uj uk l uk w uw k ~u ¼ B B þ X F  X F þ X F  v j vj kl v l kw wk v kl v k l kw X v k F w k ~ 0 ¼ 1 P Bij þ P X ik F l  P X iw F l B w ij ij ikl il ilw il k 2 Ij j V pj

^ ^ qp ð1Þ ¼ hpð2Þj^r 1 P K 12 12 jqð2Þi

P V 0 ¼ 12 ij V ijij 00 P ~ Ik 00 CIa ¼ b00 k F bk R ab

124

W. Liu et al. / Chemical Physics Letters 565 (2013) 122–127

SP ^ 12 ¼ R

X1 0 0

IJa b

4

0 0 e IJ0 0 a ^a b : R a b IJ

ð14Þ

^ y ½F ^;R ^ 12 jW0 i: þ hW0 jR 12

The choice of the projector Eq. (11), implies that we neglect the 0 ^ 12 operators like R e IJ 0 a ^ub . This can be justified, besemi-internal R ub IJ cause in the atomic case the partial wave expansion for the semiinternal part of the correlation energy terminates for finite L-quantum numbers [22]. In order to arrive at an energy expression for ic-MRCC with F12 terms and the SP ansatz, we need a stationary energy functional [18,32]. Analogous to derivation of the single reference energy expression [32], we introduce an operator with additional La0 ^ F12 ¼ 1 P cKL ðR ~ y Þa0 b a ^IJ0 0 . Employing the SP grange multipliers, K 8

ansatz these take the form

1 4

IJKLa0 b0

P

IJa0 b0

IJ

~y

KL

a0 b 0

ðR ÞIJ

^IJ0 0 a ab

ab

and we obtain

LF12 ¼ L þ DLF12 ; ^ R^ 12 jW i þ hW jK ^ R^ 12 jW i ^ eR^ 12 He e y ÞeR^ 12 He ¼ hW0 jð1 þ R 0 0 0 12    E hW0 jW0 i  1 ;

ð15Þ

^ is the short form of eT^1 T^ 2 He ^ T^ 1 þT^ 2 . Note that we put R ey where H 12 into the first term of Eq. (15) instead of the second term, although the latter possibility may seem more natural at first sight. Our choice, however, ensures that the stationarity condition

@LF12 @ cl

X ^ eff cm  Ecl ¼ 0 H ¼ lm

ð16Þ

m

  ^ R ^ þ ½H; ^ R ^ eff ¼ hUl jð1 þ R ^y Þ H ^ 12  þ 1 ½½H; ^ 12 ; R ^ 12  jUm i; H lm 12 2

ð17Þ

which directly includes the F12 correction terms and yields the corrected energy as its lowest eigenvalue E. The T^ 1 and T^ 2 amplitude equations are obtained by differentiating Eq. (15) with respect to the Lagrange multipliers

@LF12 ¼ 0; @ðk0 ÞAI @LF12 ¼ ¼ 0: @ðk0 ÞAB IJ

XIA ¼

ð18Þ

IJ AB

ð19Þ

In order to arrive at a computationally tractable theory, we employ the same set of approximations that for single-reference theory lead ^ 12 terms are kept, and to the CCSD(F12*) model [31]. Only linear R ^ 12 only couple through the zeroth-order Hamiltonian ^ y and R R 12 which is chosen as the effective Fock operator,

^ ¼ F

X

^pq ; F qp a

F qp ¼ fpq þ

X g qpuv ðc0 Þuv ;

ð20Þ

uv

pq

where ðc0 Þvu is the one-particle density matrix of the initial CASSCF wavefunction (note that the cl can be different in the final icMRCCSD wavefunction). Thus, the Lagrange function simplifies to

^ R e F12SP ¼ hW0 j½H; ^ 12 jW0 i DL ^ W i þ hW jR ^ y Hj ^ y ½F ^;R ^ 12 jW0 i þ hW0 jR 0 0 12 12 ^ R ^ þK ^ Þ½H; ^ jW i: þ hW jðK 0

1

2

12

0

ð21Þ

We note that when Eqs. (16), (18), and (19) are fulfilled, the corresponding Lagrange multipliers need not to be determined and the energy expression can instead be written as

ð22Þ

Analogous to the single-reference CCSD(F12*) method, we only keep terms up to third order in perturbation theory and some selected higher order terms (those that do not vanish under the assumptions of the standard approximation [31]). The intermediates are defined and evaluated as in single-reference theory, where the occupied orbitals are replaced by internal orbitals, see Table 2. We also note that for the U intermediate some additional considerations are necessary. In order to arrive at an intermediate that can be evaluated prior to solving the cluster equations, we adopt some elements from generalized normal ordering [34]. This allows us to approximately factorize contributions from higher-order density p r p r matrices, e.g., cpr qs  cq cs  cs cq and to absorb contributions from low-order density matrices into the U intermediate. Note, that this is also the background for defining the effective Fock operator, Eq. (20) [34]. A huge number of terms occur in the energy and residual equations, as expected for ic-MRCC theory [28]. However, all of them are analogous to terms appearing in the conventional theory and can thus be absorbed by defining appropriate intermediates. Here, we list only the leading-order terms of the energy expression (those from second-order perturbation theory),

^R ^ y ½F ^ ; ðR ^ 12 þ T^ 2 Þ þ 2H ^ 12 ÞjW0 i hW0 jðR 12 X y ab ij X y aw 1 1 ~0 þ ¼B ðC Þij t ab þ ðC Þij tijaw þ 2V 0 4 ijab 2 ijaw 1 X y ab iv X y aw iv þ cv B~ vu þ ðC Þiu tab þ ðC Þiu t aw þ 2V vu 2 ija uv iaw X

leads to an eigenvalue equation for an effective Hamiltonian

X

^ þ ½H; ^ R ^W i ^ 12 jW0 i þ hW0 jR ^ y Hj E ¼ hW0 jH 0 12

!

u

! 1 X uw ~ v x 1 X y ab v x X y ay v x v x c Buw þ ðC Þuw t ab þ ðC Þuw t ay þ 2V uw : þ 4 uv wx v x 2 ab ay ð23Þ

Similar expressions have been reported in the literature, e.g., Ref. [22,23,26]. 2.3. Implementation and calculation details The ic-MRCCSD(F12*) method has been implemented in the General Contraction Code (GeCCo), making use of automated term generation and string-based generalized tensor contractions. This approach has been employed previously in the implementation of single-reference based F12 methods [35–37] and ic-MRCC methods [28–30]. The CASSCF orbitals and the required integrals are obtained from a local modification of the DALTON program [38]. For all our test examples, we use a commutator approximation to keep the number of terms of ic-MRCCSD at a tractable level (cf. Refs. [27,28]). In this work, we truncated both the energy expression and the amplitude equations after the double commutator terms. Removal of redundant amplitudes was done as in Ref. [29]. Peterson’s correlation consistent F12-optimized basis sets [39,40] cc-pVXZ-F12 (X = D, T, Q) were employed in conjunction with the appropriate cc-pVXZ-F12-RI complementary auxiliary basis sets. In the following, we use the abbreviations DZ, TZ and QZ for these pairs of basis sets. For the Slater-type correlation factor, we take a constant value of c ¼ 1:0 throughout all our calculations. For the test systems we carried out both ‘‘conventional’’ icMRCCSD and ic-MRCCSD(T) calculations, as well as calculations with explicitly correlated geminals, denoted ic-MRCCSD(F12*) and ic-MRCCSD(F12*)+(T). The perturbative triples correction, EðTÞ , for ic-MRCCSD(F12*) is obtained by using the conventional energy expression, while the T^ 1 and T^ 2 amplitudes are calculated by

125

W. Liu et al. / Chemical Physics Letters 565 (2013) 122–127

-0.28

-0.32

Correlation energy (Eh)

each molecule. For the triatomic CH2 molecule, a numerical determination of the equilibrium structure for both the singlet and the triplet state was carried out for each method and basis set combination.

ic-MRCCSD/CAS(2,2) -’’/CAS(6,4) -’’/CAS(8,6) ic-MRCCSD(F12*)/CAS(2,2) -’’/CAS(6,4) -’’/CAS(8,6)

3. Results and discussions

-0.36

3.1. Performance on an open shell system: O2 For the O2 molecule, we have chosen three different types of complete active spaces (CAS), denoted CAS(2,2), CAS(6,4), and CAS(8,6). The CAS(2,2) consists of two electrons distributed among the two 1p orbitals, while the CAS(6,4) includes the two additional 1p orbitals and the CAS(8,6) further adds the 2r and 2r orbitals. First, we performed for this system a test of the (F12*) truncation scheme and implemented to this end the full (F12) ansatz. We calculated with both methods the total energy of the ground state of O2 employing the DZ basis set (for which the F12 corrections are largest). The energy differences between the two methods are about 0.038 and 0.089 mEh for CAS(2,2) and CAS(8,6), respectively. This energy difference is comparable to its single-reference analogues CCSD(F12*) and CCSD(F12) [19,33] and already much less than the typical accuracy of the ic-MRCCSD method [28,41]. In terms of computational effort, the (F12*) truncation greatly reduces the overhead (with respect to a conventional calculation) during the iterative solution of the coupled-cluster equations. In the present pilot implementation, however, due to short-comings of our automated intermediate recognition procedure, we still observe a small overhead of approximately 30 percent in each iteration, although it should be possible to reduce the overhead to zero, as in the single reference case [12]. The overhead should then be solely determined by the non-iterative steps, in particular the calculation of the special F12 intermediates. Next, we turn to the basis set convergence of the methods, using the (F12*) scheme throughout. For the sake of the following discussion, we define the dynamic correlation energy as the energy difference between the ic-MRCCSD energy (or ic-MRCCSD(F12*) energy, respectively) and the CASSCF energy. As shown in Figure 1,

-0.40

-0.44

-0.48

-0.52

cc-p

VDZ

cc-p

-F1

2

VTZ

cc-p

-F1

VQZ

2

-F1

2

Figure 1. Correlation energies (defined as difference to CASSCF energy, see text) of O2 at the equilibrium structure, calculated by the ic-MRCCSD and ic-MRCCSD(F12*) methods with different active spaces.

the ic-MRCCSD(F12*) method. We adopted the {3} triples model [29], which neglects those triple excitation operators that contain more than three active indices. The frozen core approximation was used in all cases. For diatomic systems, we have calculated a part of the potential energy curves (PECs) of the ground states for each above mentioned method and basis set combination. By fitting these PECs to fifth-order polynomials, we obtained the equilibrium bond distances, harmonic frequencies and the anharmonicity constant for

-150.06

ic-MRCCSD/CAS(2,2) -’’/CAS(6,4) -’’/CAS(8,6) ic-MRCCSD(F12*)/CAS(2,2) -’’/CAS(6,4) -’’/CAS(8,6)

-150.08

ic-MRCCSD(T)/CAS(2,2) -’’/CAS(6,4) -’’/CAS(8,6) ic-MRCCSD(F12*)+(T)/CAS(2,2) -’’/CAS(6,4) -’’/CAS(8,6)

Total energy (Eh)

-150.10 -150.12 -150.14 -150.16 -150.18 -150.20

cc-p

VDZ

cc-p

-F1

2

VTZ

cc-p

-F1

2

VQZ

cc-p

-F1

2

VDZ

cc-p

-F1

2

VTZ

cc-p

-F1

2

VQZ

-F1

2

Figure 2. Comparison between the calculated total energies of O2 at the equilibrium structure, employing the ic-MRCCSD and ic-MRCCSD(F12*) methods (left panel), and the ic-MRCCSD(T) and ic-MRCCSD(F12*)+(T) methods (right panel).

126

W. Liu et al. / Chemical Physics Letters 565 (2013) 122–127

Table 3 Calculated equilibrium bond lengths (r e ), harmonic frequencies (xe ), and vibrational anharmonicity (xe xe ) of the 3 R g ground state of O2, in comparison to experimental values. We denote the ic-MRCCSD(T) and ic-MRCCSD(F12*)+(T) methods as ‘‘(T)’’ and ‘‘F12*+(T)’’, respectively. O2 ; X 3

P

Basis set

g

CAS(2,2)

r e (pm)

DZ TZ QZ DZ TZ QZ DZ TZ QZ

CAS(6,4)

CAS(8,6)

Experiment [43]

xe (cm1 )

xe xe (cm1 )

(T)

F12*+(T)

(T)

F12*+(T)

(T)

F12*+(T)

122.07 121.06 120.88 122.14 121.12 120.93 122.28 121.29 121.11

121.30 120.69 120.66 121.38 120.76 120.72 121.62 120.97 120.91

1545.1 1586.0 1593.4 1536.5 1578.3 1585.8 1522.6 1563.1 1570.9

1586.3 1602.0 1602.2 1576.8 1593.7 1594.4 1555.9 1575.6 1578.1

10.2 11.3 11.1 10.9 11.6 11.0 11.5 11.9 11.6

10.7 11.1 11.0 10.9 11.4 11.3 11.5 11.8 11.7

120.75

1580.2

12.0

Table 4 Calculated equilibrium bond lengths (re ), harmonic frequencies (xe ), and vibrational anharmonicity (xe xe ) of the 1 Rþ g ground state of C2, in comparison with experimental values. We denote the ic-MRCCSD, ic-MRCCSD(F12*), ic-MRCCSD(T), and ic-MRCCSD(F12*)+(T) methods as ‘‘conv.’’, ‘‘F12*’’, ‘‘(T)’’, and ‘‘F12*+(T)’’, respectively. C2 ; X 1

Pþ g

CAS(2,2)

CAS(8,8)

Basis set

DZ TZ QZ DZ TZ QZ

Experiment [43]

r e (pm)

xe (cm1 )

xe xe (cm1 )

Conv.

F12*

(T)

F12*+(T)

Conv.

F12*

(T)

F12*+(T)

Conv.

F12*

(T)

F12*+(T)

124.13 123.76 123.56 125.26 124.93 124.73

123.60 123.46 123.42 124.81 124.68 124.62

125.14 124.80 124.60 125.25 124.91 124.71

124.52 124.46 124.45 124.80 124.66 124.59

1898.5 1914.4 1920.7 1821.3 1835.2 1841.0

1926.1 1927.8 1926.3 1844.3 1845.9 1845.6

1827.7 1842.7 1848.8 1822.2 1836.7 1842.8

1862.8 1859.4 1855.9 1845.0 1847.3 1847.4

12.6 12.1 12.1 13.9 13.5 13.5

12.3 12.0 12.0 13.6 13.5 13.5

13.7 13.3 13.3 13.9 13.5 13.9

13.1 13.1 13.2 13.5 13.6 13.7

124.25

1854.7

13.3

Table 5 Adiabatic singlet–triplet separation energies of methylene (in kcal/mol). MRCC and MRCI calculations in comparison to experimental results. The equilibrium structures are optimized individually for each method and basis set combination. We denote the ic-MRCCSD(T) and icMRCCSD(F12*)+(T) methods as ‘‘(T)’’ and ‘‘F12*+(T)’’, respectively. The MRCI and MRCI-F12(FIX) results are taken from Ref. [24]. Basis set

CAS(2,2)

CAS(6,15)

CASSCF

Conv.

F12*

(T)

F12*+(T)

CASSCF

MRCI

MRCI-F12(FIX)

DZ TZ QZ

10.85 10.41 10.40

10.31 9.26 8.99

9.44 8.88 8.80

10.36 9.36 9.11

9.50 8.99 8.92

11.21 10.61 10.54

10.35 9.33 9.07

9.36 8.94 8.86

Experiment [44]

9.12  0.20

a clearly much better basis set convergence of the correlation energy was found for the ic-MRCCSD(F12*) method compared to the conventional ic-MRCCSD method. The correlation energies computed at the ic-MRCCSD(F12*)/DZ level are already comparable to the ones at the ic-MRCCSD/QZ level of theory. Of course, this definition of the correlation energy has its limits, since the CASSCF wavefunction partially includes dynamic correlation if larger active spaces are used. In Figure 1 the absolute value of the above defined correlation energy thus decreases with the active space size. In Figure 2, we show the basis set convergence of the total energies, with and without the (T) correction. Comparing to the previous figure, it becomes apparent that the remaining basis set error is dominated by the CASSCF energy. This has been noted earlier [23,24,44] and schemes for a perturbative correction, making use of the auxiliary basis set, exist [42]. Another interesting observation is the following: including the (T) correction, the calculated total energies are nearly independent of the choice of the CAS (Figure 2, right panel). For instance, using the QZ basis the maximum energy difference due to the active spaces is around 1 mEh for both ic-MRCCSD(T) and ic-MRCCSD(F12*)+(T), while it is more than 10 mEh when the (T) correction is not included. This is in line with previous experience with the ic-MRCCSD(T) method [29] and

our example shows that this behavior is inherited by its F12 variant. Other properties (listed in Table 3), e.g., equilibrium bond lengths (r e ), harmonic frequencies (xe ), and vibrational anharmonicity (xe xe ), give further support to the excellent performance of the ic-MRCCSD(F12*)+(T) method.

3.2. Performance on a system with strong multireference character: C2 Since the ic-MRCC method is designed for systems with strong multireference character, we have chosen the C2 molecule as our next sample system. Analyzing the full-valence CASSCF wavefunction for the ground state of C2, we find that a CAS(2,2) should be sufficient. It consists of 2 electrons distributed among the 2rs and the 2rs orbitals. Results are listed in Table 4. When the (T) correction is included, the CAS(2,2) based calculations give indeed very accurate results. In particular when the (F12*) correction is applied, the results are in excellent agreement with experimental values. We also performed calculations with the full valence CAS(8,8). Here, the (T) correction has only minor impact on the calculated values, which otherwise nicely corroborate the CAS(2,2) results.

W. Liu et al. / Chemical Physics Letters 565 (2013) 122–127

3.3. CH2 singlet–triplet separation The singlet–triplet splitting of methylene (CH2) is another challenging example. We make a comparison in Table 5 between the experimental value (for the purely electronic term) and the calculated ones from different methods. For the ic-MRCC calculations, we have chosen the minimal active space, CAS(2,2), which consists of two electrons distributed among the 3a1 and the 1b1 orbitals. For MRCI based methods, we take the values from Ref. [24] where a large CAS(6,15) including Rydberg orbitals was used. The singlet (1 A1 ) and triplet (3 B1 ) equilibrium structures are optimized numerically for each method and basis set (see Supplementary information). The two sets of methods, ic-MRCC and MRCI, provide very similar results, and all agree well with the experimental value. Compared with the MRCI-F12(FIX) results, the ones of the icMRCCSD(F12*) method show a very similar basis set convergence (the energy difference between DZ and TZ is 0.52 kcal/mol for icMRCCSD(F12*) and 0.42 kcal/mol for MRCI-F12(FIX), and the one between TZ and QZ is 0.08 kcal/mol for both methods). On the other hand, ic-MRCC based methods benefit from the much smaller CAS, which gives a more compact and computationally more tractable wavefunction and also simplifies the construction and optimization of the correct CASSCF wavefunction.

4. Conclusions We have proposed and implemented an explicitly correlated icMRCCSD(F12*) method. The method employs geminal functions with a Slater-type correlation factor for describing short-range electron correlation. The amplitudes of the geminal functions are determined from the cusp conditions (SP ansatz) [13]. Based on perturbative arguments, analogous to the single-reference method CCSD(F12*) [31], we introduced a truncation scheme for the F12 correction terms, resulting in a computationally efficient and numerically accurate method. The accuracy of the ansatz was tested by comparing to the full (F12) ansatz and by studying the basis set convergence of the energies and spectroscopic parameters of O2 and C2 and the singlet–triplet splitting of CH2. Throughout our test examples, the new method icMRCCSD(F12*) greatly improves the basis set convergence for the calculated properties, compared with the conventional approach without geminal terms. The remaining basis set incompleteness errors are mainly due to the CASSCF wavefunction. Methods to handle this problem are known in the literature [42] and should be compatible with the ic-MRCCSD(F12*) method. Our results also underline that the ic-MRCCSD(T) method (which perturbatively accounts for connected three-electron clusters) [29] requires much smaller active spaces to achieve the same level of accuracy as traditional approaches, e.g., MRCI. This property is also inherited by the ic-MRCCSD(F12*)+(T) method, which is proposed in this letter. It augments the ic-MRCCSD(F12*) wavefunction with a conventional correction for the effect of connected three-electron clusters. A possibility to account for short-range correlation in connected triples has been previously discussed for single-reference methods [35] and might be of interest in the multireference context, too.

127

In summary, we believe that the presented method constitutes an important step towards the development of compact but highly accurate wavefunctions for molecular systems. Acknowledgments This work was supported by the Deutsche Forschungsgemein-schaft (DFG, Grant KO 2337/2-1). A.K. acknowledges a Heisenberg fellowship of the DFG (Grant KO 2337/3-1). Appendix A. Supplementary data Supplementary data associated with this article can be found in the online version, at http://dx.doi.org/10.1016/j.cplett.2012.12. 052. References [1] P.G. Szalay, T. Müller, G. Gidofalvi, H. Lischka, R. Shepard, Chem. Rev. 112 (2012) 108. [2] D.I. Lyakh, M. Musiał, V.F. Lotrich, R.J. Bartlett, Chem. Rev. 112 (2012) 182. [3] K. Andersson, P.Å. Malmqvist, B.O. Roos, J. Chem. Phys. 96 (1992) 1218. [4] H.-J. Werner, E.-A. Reinsch, J. Chem. Phys. 76 (1982) 3144. [5] R.J. Gdanitz, R. Ahlrichs, Chem. Phys. Lett. 143 (1988) 413. [6] D. Mukherjee, R.K. Moitra, A. Mukhopadhyay, Mol. Phys. 30 (1975) 1861. [7] A. Banerjee, J. Simons, Int. J. Quantum Chem. 19 (1981) 207. [8] B. Jeziorski, H.J. Monkhorst, Phys. Rev. A 24 (1981) 1668. [9] W. Kutzelnigg, Theor. Chim. Acta 68 (1985) 445. [10] W. Klopper, W. Kutzelnigg, Chem. Phys. Lett. 134 (1987) 17. [11] W. Kutzelnigg, W. Klopper, J. Chem. Phys. 94 (1991) 1985. [12] C. Hättig, W. Klopper, A. Köhn, D. Tew, Chem. Rev. 112 (2012) 4. [13] S. Ten-no, Chem. Phys. Lett. 398 (2004) 56. [14] W. Klopper, R. Röhse, W. Kutzelnigg, Chem. Phys. Lett. 178 (1991) 455. [15] J. Noga, W. Kutzelnigg, W. Klopper, Chem. Phys. Lett. 199 (1992) 497. [16] A.J. May, F.R. Manby, J. Chem. Phys. 121 (2004) 4479. [17] D.P. Tew, W. Klopper, J. Chem. Phys. 123 (2005) 074101. [18] T.B. Adler, G. Knizia, H.-J. Werner, J. Chem. Phys. 127 (2007) 221106. [19] A. Köhn, J. Chem. Phys. 130 (2009) 104104. [20] R.J. Gdanitz, Chem. Phys. Lett. 210 (1993) 253. [21] R.J. Gdanitz, Chem. Phys. Lett. 283 (1998) 253. [22] S. Ten-no, Chem. Phys. Lett. 447 (2007) 175. [23] T. Shiozaki, H.-J. Werner, J. Chem. Phys. 133 (2010) 141103. [24] T. Shiozaki, G. Knizia, H.-J. Werner, J. Chem. Phys. 134 (2011) 034113. ˇ a, S. Ten-no, J. Pittner, J. Noga, Phys. Chem. Chem. [25] O. Demel, S. Kedzˇuch, M. Švan Phys. 14 (2012) 4753. [26] R. Haunschild, S. Mao, D. Mukherjee, W. Klopper, Chem. Phys. Lett. 531 (2012) 247. [27] F.A. Evangelista, J. Gauss, J. Chem. Phys. 134 (2011) 114102. [28] M. Hanauer, A. Köhn, J. Chem. Phys. 134 (2011) 204111. [29] M. Hanauer, A. Köhn, J. Chem. Phys. 136 (2012) 204107. [30] M. Hanauer, A. Köhn, J. Chem. Phys. 137 (2012) 131103. [31] C. Hättig, D.P. Tew, A. Köhn, J. Chem. Phys. 132 (2010) 231102. [32] D.P. Tew, W. Klopper, C. Hättig, Chem. Phys. Lett. 452 (2008) 7. [33] A. Köhn, D.P. Tew, J. Chem. Phys. 133 (2010) 174117. [34] W. Kutzelnigg, D. Mukherjee, J. Chem. Phys. 107 (1997) 432. [35] A. Köhn, J. Chem. Phys. 130 (2009) 131101. [36] M. Hanauer, A. Köhn, J. Chem. Phys. 131 (2009) 124118. [37] A. Köhn, J. Chem. Phys. 133 (2010) 174118. [38] DALTON, A Molecular Electronic Structure Program, Release 2.0, 2005. Available from: . [39] K.A. Peterson, T.B. Adler, H.-J. Werner, J. Chem. Phys. 128 (2008) 084102. [40] K.E. Yousaf, K.A. Peterson, J. Chem. Phys. 129 (2008) 184108. [41] F.A. Evangelista, M. Hanauer, A. Köhn, J. Gauss, J. Chem. Phys. 136 (2012) 204108. [42] L. Kong, E.F. Valeev, J. Chem. Phys. 133 (2010) 174126. [43] K.P. Huber, G. Herzberg, Molecular Spectra and Molecular Structure IV. Constants of Diatomic Molecules, Van Nostrand, Princeton, 1979. [44] I. Shavitt, Tetrahedron 41 (1985) 1531.