The local CC2 model

The CC2 model has been proposed by Christiansen and coworkers [21]

as an approximation to the CCSD model, where the singles are treated as zeroth-order parameters in terms of the fluctuation potentialW(Wis a dif-ference between the HamiltonianHand the Fock operatorF). The doubles excitations are treated at first order inW(as in MP2). This parametrization was rationalized in Ref. 21 by the fact that singles respond to an exter-nal perturbation in zeroth W order. To keep the singles to all orders is crucial for response theory where un-relaxed orbitals are usually used.

The orbital-relaxed response approach leads to a pole structure inconsis-tent with exact theory, as discussed e.g. in Ref. 20. The equations for the ground-state CC2 amplitudes are with the exponent of the singles cluster operator, and ˜Ψ^a

i, ˜Ψ^ab

ij are con-travariantconfiguration state functions [51] projecting onto the singles and doubles manifold. The singles and doubles operators are described in terms of usual single and double orbital substitution operators [52],

T1 = X We denote the occupied orbitals by indicesi, j,k,land the virtual (or PAO) orbitals by a,b,c,d. We have also tested a slightly modified model, where the similarity transformed Fock operator ˆF(for a definition of ˆFsee Ref. 53)

was used also in the commutator of eq. (2.2). We will use in the following the acronymCC2-bfor this modification of the CC2 model. For CC2-b the doubles-doubles block of the Jacobian is no longer diagonal in a canonical MO basis, what is an essential feature for previous CC2 implementations.

On the other hand the doubles-doubles block of the local CC2 Jacobian is non-diagonal anyhow, hence local CC2-b does not involve any additional cost relative to local CC2.

To devise a local CC2 method, orthogonal localized molecular orbitals and mutually non-orthogonal projected atomic orbitals are introduced to span occupied and virtual orbital space, respectively. The former are con-structed from the parental occupied canonical Hartree-Fock orbitals by an unitary transformation according to a localization criterion [54,55], the lat-ter by projection of the AOs onto the virtual space. By construction the PAOs can still be assigned to individual atoms. In the spirit of the original CC2 idea no local restrictions are applied for the singles, whereas double excitations from LMOsiand jare restricted to pair-specific subspaces (do-mains) [i j] of the virtual space. These pair domains contain only PAOs belonging to a certain subset of atoms in the spatial vicinity of the cor-responding LMOs and can be determined, e.g. by following a procedure devised by Boughton and Pulay [56]. In addition to the truncation of the virtual space, restrictions on the pair list (based e.g. on a distance criterion between the related LMOs) can be imposed, such that the number of re-maining doubles amplitudes scales asymptotically linearly with molecular size.

By employing diagrammatic techniques, working equations for local CC2 in terms of amplitudes and integrals are conveniently obtained from eqs. (2.1,2.2) as

whereSandFare the PAO overlap and Fock matrices and

T˜^ij =2T^ij−T^ji. (2.6)

Here and in the following we employ the convention to denote matrices and vectors in the virtual space by capital and small bold letters, respec-tively, and scalars by normal (no-bold) letters [57]. Although not shown explicitly, the summations in eqs. (2.4) and (2.5) involving double ampli-tudes are restricted to their respective domains and pair lists.

In eqs. (2.4,2.5) modified (dressed) integrals (all objects with a hat) appear due to the similarity transformation ofH with the exponent ofT1. Such dressed integrals are defined as (see, e.g. [53])

(pqˆ|rs) X

µνρσ

(µν|ρσ)Λ^p_µpΛ^h_νqΛ^p_ρrΛ^h_σs, (2.7) where p,q, . . . and µ, ν, . . . denote general MO (LMO or PAO) and AO indices, respectively. The dressed Fock integrals are defined as

fˆpq hˆpq+X

k

(2(kkˆ|pq)−(kqˆ|pk)). (2.8) In eq. (2.7) modified MO coefficient matrices

Λ^p h

L|| (P−L^†t^†₁S)i

, Λ^h h

(L+Pt1)|| Pi

(2.9) are used in the integral transformation instead of the usual LMO/PAO coefficient matrices Land P(see e.g. [48]). Here the ”k” sign means that matrices on the left (n_AO ×n_LMO) and on the right (n_AO×n_PAO) are glued together,t1denotes thenPAO×nLMOmatrix of the singles amplitudes, where n_LMOis the number of occupied andn_PAO =n_AO– the number of projected atomic orbitals. The Lmatrix is the rectangular matrix of LMO transfor-mation coefficients, while the P matrix is the projection matrix from the AO to PAO basis. Therefore, theΛ^pand Λ^hmatrices are rectangular with the dimensions nAO ×(nLMO +nPAO). Due to the structure of the integral transformation in eq. (2.7) the permutational symmetry between orbitals related to the same electron is lost. In a diagrammatical context, the first and the second orbital index of an electron correspond to anoutgoingand an incoming line, respectively, relative to the interaction vertex. Thus, in eq. (2.7) the indicesp,rcorrespond to outgoing, andq,sto incoming lines.

Taking into account the structure of the modified MO coefficient matrices, eq. (2.9), it is immediately clear, that only quasi-particle creation lines, i.e., lines connecting the interaction vertices with the bra side, are effectively dressed. Therefore, the last term of eq. (2.4) involving dressed 3-external integrals can be conveniently rewritten in terms of bare (undressed)

inte-grals,

The only remaining dressed integrals occurring in eqs. (2.4,2.5) are (ˆfⁱ)a fˆai = X

k

2(aiˆ|kk)−(akˆ|ki)

, (2.11)

( ˆK^ij)ab Kˆ^ij_ab = (aiˆ|bj), and ( ˆk^ijk)a kˆ^ijk_a =(iaˆ|jk).

The computational cost for obtaining the singles and doubles residuals in eqs. (2.4,2.5) can be significantly reduced by invoking the DF approxima-tion for the ERI integrals [32],

(pqˆ|rs)=X are 2-index and dressed 3-index ERIs, and ˆc^P_rs the corresponding (dressed) fitting coefficients. With the definition

V_ia^P

eq. (2.4) can be rewritten in the form Ωⁱ

Note that the summation over the PAO indexcin the last term of eq. (2.14) is restricted to the united pair domain [i]U = ∪[ik],∀(ik), i.e., the union of all pair domains with a fixed LMO i. The quantity V_ia^P is identical to the important intermediate in DF–LMP2 gradient calculations, defined in

eq. (51) of Ref. 43 and involves only bare fitting coefficients. The virtual-occupied part of the dressed Fock matrix defined in eq. (2.11) is calculated as

fˆµν = 2X

P

ˆ

cP(P|µν)−X

kP

(µkˆ|P)c^P_kν, with (2.15) ˆ

cP X

kν

c^P_kνΛ^h_νk.

The doubles amplitude equations (2.5) are virtually identical to the local MP2 equations with the sole difference that dressed rather than bare ex-change integrals ˆK^ij_abappear. The latter are obtained according to eq. (2.12) in the usual way [58]. For CC2-b, in addition, also the Fock matrix elements are replaced by their dressed counterparts.

The DF–LCC2 amplitude equations (2.14,2.5) are solved iteratively us-ing DIIS convergence acceleration [59]. Redundancies in the PAO basis are eliminated at the stage of forming the amplitude updates from the residual vectors Ωⁱ,Ω^ij, as described elsewhere [37, 40]. As it is evident from eqs. (2.14,2.15), only dressed 3-index integrals with at least one index transformed to the LMO basis do occur. In particular the bare 2-external integrals (P|ab) can be generatedoutsidethe CC2 iterations. The use of this integral set could be entirely avoided at the price of back-transforming the PAO index ofV^P_iato the AO basis, carrying out the contraction with the inte-grals (the last term of eq. (2.14)) entirely in the AO basis, and transforming the result again to the PAO basis in each iteration. Assuming 15 iterations until convergence the break even point for this alternative algorithm with respect to the CPU time would be arrived at a basis set size with a ratio nAO/nLMO ≈ 30, which corresponds roughly to the aug-cc-pVTZ basis. As a further disadvantage, the domain restriction in the contraction step is lost in the AO basis.

The formally most expensive step is the formation of theV_ia^P intermedi-ate according to eq. (2.13) with a nominal scaling of ∝n_Pn²_AOn_FIT =O(N⁵) (n_P: number of pairs, n_FIT: number of fitting functions), yet by virtue of the local approximation this scaling reduces to ∝ nPL²nFIT = O(N³) (L is the average size of a pair domain). Restrictions in the pair list reduces the scaling further toO(N²), and linear scaling can be achieved by introducing local fit domains [42].