Hybrid and Constrained Resolution-of-Identity ... - ACS Publications

Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Article

Hybrid and constrained resolution-ofidentity techniques for Coulomb integrals Ivan Duchemin, Jing Li, and Xavier Blase J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b01215 • Publication Date (Web): 17 Jan 2017 Downloaded from http://pubs.acs.org on January 23, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33

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

Journal of Chemical Theory and Computation

Hybrid and constrained resolution-of-identity techniques for Coulomb integrals Ivan Duchemin,∗,† Jing Li,‡ and Xavier Blase‡ INAC, SP2M/L_Sim, CEA/UJF Cedex 09, Université Grenoble Alpes, 38054 Grenoble, France, and Université Grenoble Alpes, CNRS, Inst NEEL, F-38042 Grenoble, France E-mail: [email protected]

∗ To

whom correspondence should be addressed SP2M/L_Sim, CEA/UJF Cedex 09, Université Grenoble Alpes, 38054 Grenoble, France ‡ Université Grenoble Alpes, CNRS, Inst NEEL, F-38042 Grenoble, France † INAC,

1 ACS Paragon Plus Environment


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

Graphical TOC Entry


Page 2 of 33

Page 3 of 33

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


Abstract The introduction of auxiliary bases to approximate molecular orbital products has paved the way to significant savings in the evaluation of four-center two-electron Coulomb integrals. We present a generalized dual space strategy that shed a new light on variants over the standard density and Coulomb-fitting schemes, including the possibility of introducing minimization constraints. We improve in particular the charge- or multipole-preserving strategies introduced respectively by Baerends and Van Alsenoy that we compare to a simple scheme where the Coulomb metric is used for lowest angular momentum auxiliary orbitals only. We explore the merits of these approaches on the basis of extensive Hartree-Fock and MP2 calculations over a standard set of medium size molecules.

1

Introduction

Since the electronic Hamiltonian only contains one- and two-body operators, its expectation value over Slater determinants leads to the well-known one- and two-electron integrals that stand at the heart of all quantum chemistry code. Taking as a paradigmatic example the Hartree and exchange terms involving the product of four one-body spin-orbitals, one obtains a set of four-index electronrepulsion integrals (ERI) that needs to be either stored or re-calculated on the fly, in order to be repeatedly used in various forms depending 1,2 on the chosen formalism (HF, DFT, MP2, CC, etc.) Alternatively, with the use of atom-centered bases such as Gaussian atomic orbitals (AO), these four-indices ERIs can be written in terms of the standard four-center AO integrals of which the number increases rapidly with the basis size, in particular in the case of correlated calculations requiring AO basis with large ζ -multiplicity. While the localized nature of the Gaussian AO basis allows to somehow reduce the number of integrals to be stored and calculated, with the help of elegant pre-screening detection criteria, 3 it remains however that the number of relevant fourcenter integrals is extremely large, scaling as O(N 2 ) with N the number of AO in the system, leading to substantial memory requirements. In order to reduce such memory constraints, and further considerably reduce the computational 3 ACS Paragon Plus Environment


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

cost associated with the use of such integrals in the evaluation of Hartree, exchange or higherorder correlation energy terms, efficient techniques have been introduced, including in particular the resolution-of-the-identity (RI) family of approaches. 4–12 The RI approaches proceed by approximating the four-centers ERI by a combination of three- and two-centers integrals, relying on the approximate expansion of AO direct products over auxiliary (or fitting) basis functions (ABFs). Unless an efficient error correction scheme is used, such as in the original Van Alsenoy paper, 9 the chosen ABF set generally needs to be larger than the AO basis used to perform the desired calculations, as they must span the resulting molecular orbitals product space. Many studies have been dedicated to the design of optimized and accurate auxiliary Gaussian bases, allowing to preserve a remarkable accuracy with a reasonable ABF set size. 13–17 While alternative approaches such as Cholesky decomposition 18,19 have been developed to reduce the number of ERIs, RI techniques remain central to a large number of modern quantum chemistry calculations. In particular, it has lead to the computationally efficient RI-MP2 approach. 10,14,20,21 The fitting procedure that can be used to approximate the direct product of two molecular orbitals by a sum over ABFs, such as Gaussian bases, is not unique, 14–16,22 and complementary strategies such as robust fitting allow to further reduce the error induced by the fit. 23–28 As already described in literature 8,29 and re-detailed below, the metric used to minimize the error introduced by the fitting procedure leads to different approximations, such as the density-fitting or the Coulomb-fitting procedure, allowing to balance the accuracy to computer cost ratio. It is now recognized that the Coulomb-fitting scheme that has been developed to minimize the residual density coulomb energy rather than the residual density itself, leads consistently to better results than the density-fitting scheme for a given ABF set. 11,23 However, Coulomb fitting comes with the disadvantage that three-center ERI need to be either re-calculated or stored, instead of the less numerous three-center overlap integrals associated with the density-fitting scheme (see Figure 1 for an illustration). Even though providing a considerable improvement over four-center ERI, the three-center Coulomb integrals still represent a significant memory requirement. In particular, the Coulomb fitting scheme retains a O(N 2 ) scaling in the number of non-zero three center integrals,


Page 4 of 33

Page 5 of 33

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


reducing only their prefactor significantly. The quality of the auxiliary basis is also of great importance when performing different RI approximations. The differences between density and Coulomb fitting procedures come from the fact that both approach behave differently when the exact density deviates from the ones spanned by the auxiliary basis: in principle, all the fitting techniques are equivalent when the density being fitted is exactly contained within the auxiliary basis span. Exploring the merits of a fitting procedure involves thus consideration of its accuracy, balanced with its computational cost, both depending on the size and quality of the auxiliary basis used. In this study, we review two variations, and their combination, over the density-fitting and Coulomb-fitting standard schemes. We wish in particular to explore the quality of a constrained density-fitting scheme where the error minimization scheme includes a generalization of Baerends charge-preserving 5 scheme, that has been further extended to multipole-preserving constraints by Van Alsenoy. 9 In a second direction, we study a hybrid scheme that lays in between densityand Coulomb-fitting approaches by differentiating the treatment of ABF (l=0) angular momentum components from higher order terms. These different approaches are tested within standard variations over the auxiliary basis sets by exploring the impact of using spherical or Cartesian auxiliary basis functions, applied with the standard cc-pVTZ-RI and aug-cc-pVTZ-RI basis sets 30–33 in conjunction with the standard cc-pVTZ AO basis. 30 We demonstrate in particular how, even when associated with larger augmented auxiliary basis sets, both of these approach can lead to a significant decrease in the number of required three-center integrals while maintaining an excellent accuracy. Our work is organized as follows: the theory part starts with a brief review of the well known density and Coulomb fitting formalisms (Section 2.1). We then introduce in Section 2.2 an alternative and more general point of view on fitting techniques, namely a dual space analysis, that encloses both cases and allows to treat on the same footing the choice of the fitting metric and the imposed constraints. This alternative approach leads in particular to define two computationally advantageous variations on the standard fitting techniques, namely an hybrid density and Coulomb



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

Page 6 of 33

fitting approach, and an improved multipole-preserving constrained technique that generalizes Van Alsenoy’s original formalism. 9 The introduced techniques are tested (Section 3) on the basis of Hartree-Fock and correlated MP2 calculations performed on the standard set of molecules proposed by Thiel and coworkers. 34

2 2.1

Theory Standard RI approaches

The standard resolution-of-the-identity (RI) techniques can be introduced in an attempt to express the product of two molecular orbitals ρi j (r) = φi (r)φ j (r) over an auxiliary basis set labeled {β }, namely by finding the optimal cβ coefficients such that:

ρ(r) ' ∑ cβ β (r)

(1)

β

where we droped the (i j) indices for brievety. The least-square fit of the (co)density ρ by the auxiliary basis starts with the minimization with respect to the cβ coefficients of the integrated residual error: 2

|| δ ρ || =

Z

2 dr ρ(r) − ∑ cβ β (r)

(2)

β

leading to the so-called density-fitting RI (RI-SVS) approach : cβ = ∑[S−1 ]β β 0 h β 0 | ρ i β0

hβ0 | ρ i =

Z

(3) dr β 0 (r)φi (r)φ j (r)

where S−1 denotes the inverse of the overlap matrix S of the auxiliary basis, i.e., with coefficients [S]β β 0 = h β | β 0 i. The second most common RI approach proceeds by minimizing the error introduced by the fitting procedure directly on the Coulomb integrals rather than on the (co)densities. In particular, by minimizing the Coulomb self-interaction ( δ ρ || δ ρ ) of the residual density δ ρ, 6 ACS Paragon Plus Environment

Page 7 of 33

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


where the notation ( · || · ) indicates a Coulomb integral, 4 one obtains the so-called Coulomb-fitting approach (RI-V) with: cβ = ∑[V −1 ]β β 0 ( β 0 || ρ ),

(4)

β0

where V −1 is now the inverse of the auxiliary basis Coulomb interaction matrix with coefficients [V ]β β 0 = ( β || β 0 ). Once two ρ and ρ 0 (co)densities optimal fitting coefficients have been obtained, the ERI can be approximated by: ( ρ || ρ 0 )

RI−SV S

'

∑0h ρ | β i [S−1V S−1]β β 0 h β 0 | ρ 0 i

(5)

∑0( ρ || β ) [V −1VV −1]β β 0 ( β 0 || ρ 0 )

(6)

ββ RI−V

'

ββ

(7)

While the RI-V expression can be clearly simplified (with V −1V = I), the notations above emphasize the formal similarities between the two approaches. This similarity also appears within the Fourier representation of the fitting problem, where the density and Coulomb equations only differ through the weights of the fitting coefficients. 35 In the case of an auxiliary basis composed, e.g., of decaying Gaussian orbitals, and due to the non-locality of the Coulomb operator, sparsity can be much less exploited in the Coulomb-fitting approach than in the density-fitting scheme, leading to a much larger number of coefficients to be stored in the RI-V approach, with an O(N 2 ) scaling of the number of non-zero ( β || ρ ) coefficients, instead of O(N) non-zero h β | ρ i coefficients in the density-fitting scheme. However, it has been shown in many studies that the RI-V approach is far superior to the usual RI-SVS density-fitting scheme in terms of accuracy, as emphasized below.

2.2

A generalized RI approach

An alternative interpretation of these two approaches is that one is inserting a projector Iβ onto the auxiliary basis in between the (co)densities and the Coulomb operator. Formally, Iβ is a linear retraction of the (co)densities functional space, encompassing both AO direct products and ABF 7 ACS Paragon Plus Environment


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

Page 8 of 33

functions, into the functional space supported by the ABF. As such it is the identity map (i.e., a closure relation) of the ABF space, and is idempotent. Following this definition, a construction of Iβ must thus satisfy, inter alia, the following property:

Iβ | β i = | β i

(8)

The definition of Iβ as a linear operator requires the choice of a set of linear forms h γ | that provides enough distinct measures to discern all elements of the {β } set. These linear forms defines a particular dual space that we associate to the ABF space, and by extension we will call γ a dual function. A generic expression for Iβ is then: Iβ = ∑[A]β γ | β ih γ | βγ

(9)

A = (T † Γ)−1 T † ,

[Γ]γβ = h γ | β i

with T any matrix so that T † Γ can be inverted. If {β } and {γ} share the same cardinality, namely |{β }| = |{γ}| in our notations, Γ and T are square matrices and the above expression simplifies with (T † Γ)−1 T † = Γ−1 to give: Iβ = ∑[Γ−1 ]β γ | β ih γ |

(10)

βγ

In this particular case, the dual function set {γ} defines the ensemble of measures that are left invariant under application of the operator Iβ over an arbitrary (co)density: h γ | ρ i = h γ | Iβ | ρ i. However, this properties does not hold exactly when one considers a supernumerary {γ} set, leading in particular to rectangular Γ and T matrices. In this case, enforcing a minimal deviation of the fitted measures with respect to the original ones leads to the choice T = Γ, recovering the


Page 9 of 33

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


expression of a standard least square estimator: ∂ ∂T ⇒

∑h γ | ρ − Iβ ρ i2 = 0 γ

(11)

Iβ = ∑[(Γ† Γ)−1 Γ† ]β γ | β ih γ | βγ

We will refer to criterion (11) as the minimal deviation condition, as it plays a central role in the following development. The minimization of the dual measure variations clearly outlines the idea that changing the dual basis leads to a different "measure" of the quality of the fit. As an example, one can consider the case |{β }| = |{γ}| and introduce Iβ on both sides of the Coulomb operator v, leading to: ( ρ || ρ 0 ) = h ρ | v | ρ 0 i ' h ρ |I†β v Iβ | ρ 0 i

(12)

' ∑h ρ | γ i [Γ−1†V Γ−1 ]γγ 0 h γ 0 | ρ 0 i γγ 0

By identification of expressions (5) and (6) to expression (12), we verify in particular that both the RI-SVS and RI-V approaches differ only through the choice of a different set of dual functions {γ}: in the RI-SVS approach, the choice {γ(r)} = {β (r)} is made, while in the RI-V approach, one chooses {γ(r)} = { dr0 β (r0 )/|r − r0 |}. As expected, the RI-V approach favors a projection R

that preserves the Coulomb interactions of the fitted density with element of the auxiliary basis: ( β || ρ ) = ( β || Iβ ρ ), recovering the robust fitting variational criterion on the Coulomb energy approximation described for example in Ref. 23.Contrary to the Coulomb-fitting RI-V approach, the RI-V0 approach cannot be considered robust, namely leading to an error on the Coulomb integrals that is in principle not quadratic with respect to the error on the fitted density. As shown below however, the induced error remains extremely limited, so that "robust corrections" would destroy the computational gain for a marginal accuracy increase.



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

2.2.1

Page 10 of 33

Adjusting the linear forms: the RI-V0 scheme

Here we momentarily restrict ourselves to the case |{β }| = |{γ}|, namely the dual and auxiliary spaces have identical size, so that the Coulomb integrals are approximated through expression (12). Beyond the two canonical RI-SVS and RI-V techniques, different RI schemes can be explored by modifying the dual functions as discussed above. Besides the choices {γ(r)} = {β (r)} or {γ(r)} = { dr0 β (r0 )/|r − r0 |}, let’s mention also for sake of illustration the choice of {γ(r)} = R

{ dr0 β (r0 ) erfc(−µ|r − r0 |)/|r − r0 |} where the complementary error function (erfc) allows to R

tune the accuracy and sparsity between the two limiting RI-SVS and RI-V schemes (see Ref. 36), but three-indices "attenuated" Coulomb integrals must be calculated. The dual-space formulation also shed another light on the reasons for the better performance of the RI-V approach with respect to the RI-SVS: since the dual in the RI-SVS approach is the auxiliary basis itself, the fitting procedure is blind to any residual which is not contained within the auxiliary basis set itself, even if it may contribute significantly to the Coulomb integrals. As an example, the monopolar component of a density will be truncated to keep only its projection onto the auxiliary basis. Exploiting the flexibility of the dual-space formulation, and the above-mentioned analysis that the dual basis defines the measures that are invariant during the fitting process, we introduce the following "hybrid" set of test functions:    ( r || β ) for l = 0 β components γ(r) =   β (r) otherwise

(13)

We use the notation RI-V0 for this scheme where the Coulomb-fitting like γ(r) = ( r || β ) test functions are only generated from the (l=0) angular-momentum components (e.g., s orbitals) of the Cartesian (CAB) or spherical (SAB) auxiliary {β } basis set. The {γ} dual function set is then completed directly with the remaining β (r) functions. In particular, we draw the reader attention


Page 11 of 33

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


on the fact that the quantities preserved through such a fitting procedure are:    ( β || ρ ) = ( β || I ρ ) for l = 0 β components β   h β | ρ i = h β | Iβ ρ i otherwise

(14)

so that the Coulomb interaction of the densities with monopolar components of the auxiliary basis set is preserved. We will see in consequence that this approach leads to a dramatic improvement over the standard RI-SVS approach. Such a scheme, intermediate between the density-fitting and the Coulomb-fitting approaches, presents the advantage of being straightforward to implement in standard codes. Let us stress out once again that this scheme implies the calculation of (αα 0 ||β ) three-center Coulomb integrals, where {α} is the AO basis used to express the molecular orbitals, only for the auxiliary basis (l=0) components. For the sake of illustration, this represent 8 orbitals out of 81 in the case of the carbon element cc-pVTZ-RI basis, 15 leading to a drop of the O(N 2 ) prefactor by ca. an order of magnitude. Moreover, augmenting the auxiliary basis with diffuse orbitals only brings one more s component, having thus very little impact on this prefactor. Thus, while the number of (αα 0 ||β ) three-center Coulomb integrals still scales as O(N 2 ) within the RI-V0 approach, a scaling identical to that of the standard RI-V, the advantageous associated prefactor makes it a computationally interesting alternative, even with an augmented aug-cc-pVTZ-RI basis. This is clearly illustrated in Figure 1 where we compare the number of non-zero three-center integrals associated with various approaches. Finally, it is worth mentionning that in the case of the linear RI-SVS approach, computation will eventually be be bound by the calculation of 2-center Coulomb integrals that scales as O(N 2 ). The corresponding numbers of integrals are also reported in Figure 1 (dashed line), showing clearly however that they remain marginal as compared to 3-center integrals up to rather large system sizes.



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

'

Figure 1: Number of three-center integrals involved through the different fitting schemes over the acene series, for a non-zero integral cut-off value set to 10−9 Hartree. We compare here the Coulomb (RI-V) and density (RI-SVS) fitting schemes using the Cartesian cc-pVTZ-RI auxiliary basis, with the standard and constrained (up to fourth order momenta µ4 ) RI-SVS scheme and the RI-V0 scheme using the augmented Cartesian aug-cc-pVTZ-RI auxiliary basis. The corresponding number of 2-center Coulomb integrals between elements of the Cartesian cc-pVTZ-RI auxiliary basis is also reported. 2.2.2

Supplementing linear forms: the constrained approach

A second possibility to be explored in order to enhance the quality of the fit is to expand the set of dual functions, so as to enforce preservation of important properties of the density that are not described through the measures provided by the usual {γ} sets. For instance, preservation of monopole, and later on monopole and dipole, of the fitted density were suggested successively by Baerends 5 and Van Alsenoy. 9 Preservation of low order momenta of the density proves indeed to be of great importance when tackling Coulomb interactions. 37 However, equation (11) tells us that such quantities cannot be strictly preserved by simply adding the corresponding measures to the dual function set. The strict respect of the density low order momenta, or any other measure over the density, can be however obtained through the inclusion of constraints imposed during the fitting procedure. In the case of Baerends and Van Alsenoy, this was done through the addition of a Lagrange multiplier to the fitting least square equation. 5,9 It turns out that the closure relation Iβ interpretation provides a general framework of which Baerends and Van Alsenoy procedures are specific cases. In the


Page 12 of 33

Page 13 of 33

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


present formulation, the constraints appear in the extended relation set that the closure relation must satisfy: Iβ | β i = | β i ,

h µ | Iβ = h µ |

(15)

where the set of dual functions {µ} is chosen here so as to enforce the preservation of the low order momenta of the density, i.e.,

µi jk (r) = xi y j zk

(16)

Our choice of low order momenta measures over low order multipoles comes from the fact that the former provides us with a set of measures which is somewhat invariant with respect of the chosen origin. Considering now the extended generic expression for Iβ :

Iβ = ∑[A]β γ | β ih γ | + ∑[B]β µ | β ih µ |, βγ

(17)

βµ

the coefficient [A]β γ and [B]β µ are set by the fulfillment of relations (15): Iβ | β i = | β i ⇒ ∑[A]β γ h γ | β 0 i + ∑[B]β µ h µ | β 0 i = δβ β 0 γ

h µ | Iβ = h µ | ⇒

µ

     ∑h µ | β i[A]β γ = 0

(18)

β

    ∑h µ | β i[B]β µ 0 = δµ µ 0 β

Writing M the matrix of coefficient [M]µβ = h µ | β i, the conditions (18) can be rewritten as a set or three linear equations in matrix form: AΓ + BM = I MA = 0 MB = I


(19)


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

Page 14 of 33

which leads to a general solution of the form: B = T † (MT † )−1 A = I − T † (MT † )−1 M Γ−1

(20)

where T can be any matrix that allows the evaluation (M · T † )−1 . In particular, one can show that Van Alsenoy procedure correspond to a specific choice of T , namely writing the fitting problem Lagrangian as: Z

L(cβ , λµ ) =

2 dr ρ(r) − ∑ cβ β (r) β

(21)

− ∑ λµ h µ | ρ − ∑ cβ β i µ

β

is equivalent to choosing (see appendix A for demonstration):

T = M(Γ† Γ)−α with α =

1 2

(22)

The columns of the matrix T † can be seen as the directions of the fit, i.e., the corrections brought to the fitted coefficients in the auxiliary basis, that are adjusted in order to enforce the constraints: they are removed from the fit through the A component and re-added explicitly through the B ones, given the exact measure h µ | ρ i over the density being fitted. A particular choice of T determines thus the deviation of the {h γ | Iβ ρ i} measure set with respect to the original {h γ | ρ i} one. The minimal deviation condition over h γ | ρ i measures (11) now leads to a different choice of T than Baerends and Alsenoy ones, namely (see appendix B for demonstration): ∂ ∂T ⇒

∑h γ | ρ − Iβ ρ i2 = 0 γ

(23)

T = M(Γ† Γ)−α with α = 1

We devote the next section to explore the performances of both choice (22) and (23) by assessing


Page 15 of 33

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


the error resulting from the fitting procedure using T = M(Γ† Γ)−α as a function of the α parameter, tested on a set of 28 molecules designed by Thiel and coworkers. 34,38–42 Of course, this expression is not the only possible choice for T , and one could ask in principle for an optimal solution. The answer depends strongly on the (co)densities test ensemble. It turns out that in absence of further considerations on the fitted (co)densities, the minimal deviation condition over h γ | ρ i measures provides a sensible α = 1 condition that brings in practice lower errors than the original Alsenoy’s α = 1/2 choice.

2.2.3

Validation of the α = 1 scheme

1

Figure 2: Total energy Mean Average Error (MAE) per atom (not counting H atoms) for selfconsistent Hartree-Fock calculation over Thiel’s set, with respect to exact calculations, using the cc-pVTZ Kohn-Sham basis and Cartesian cc-pVTZ-RI ABF coupled with the constrained RI-SVS scheme. µi denotes the enforcement of density momenta up to order i.

In order to demonstrate the best choice for α within the functional form (23), we focus here on self-consistent Hartree-Fock total energy calculations performed with the RI-SVS fitting procedure. As it can be expected, this scheme is the one that benefits the most from the enforcement of low order momentum constraints. Exact Hartree-Fock energies, namely without any RI approximations, are provided by the NWChem package. 30–33 Following the NWChem convention, the cc-pVTZ atomic orbitals basis is taken to be spherical while the corresponding cc-pvtz-RI is taken to be Cartesian. 30–33



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

As discussed above, constraints are applied to the (co)density momenta h xi y j zk | ρ i. The notation µl signifies that all momenta up to order i + j + k = l were enforced. Also, although it was not found to be the best normalization scheme, 43 for this specific case study we use auxiliary and dual basis sets that are both L2-normed. This facilitate indeed the exercise by making Γ symmetric, simplifying thus the computation of Γ−α , without significant loss of accuracy on the resulting total energies. We represent in Figure 2 the mean absolute error (MAE) obtained for the Hartree-Fock energy over Thiel’s set of molecules within the RI-SVS scheme combined with the enforcement of (co)density momenta as a function of the α coefficient. It appears clearly that α = 1 should be close to the optimal choice, independently of the maximum constraint order enforced. 44 The deterioration of the fit for low α values is very important when raising the maximum constraint order, and Alsenoy’s choice (α = 1/2) does not behave monotonically when increasing the constraint order. On the contrary, at α = 1, the error was found to decrease systematically when increasing the maximum constraint order. Based on these observation, we choose to organize the next Sections focusing on the result obtained with the (α = 1) choice resulting form the minimal deviation condition (23). As a final remark, we note that the number of constraints enforced is relatively independent of the total number N of AOs, and the additional number of two- and three-center integrals to be calculated, namely h µ | β i and h µ | αα 0 i, scales thus linearly with N. As a result, and as seen in Figure 1 where we compare the constrained and unconstrained RI-SVS approaches, the overall additional cost in terms of three-center integrals is negligible. We also emphasize that the constrained approach is of course compatible with the RI-V0 one, and the results obtained with mixed schemes are presented and discussed in the next Sections.


Page 16 of 33

Page 17 of 33

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


3

Applications

Our calculations are performed with a dedicated C++ implementation of the presented formalism. For validation and comparison, reference calculations, including in particular calculations without any RI approximation, were performed with the NWChem package. 45 We focus here below on self-consistent Hartree-Fock and perturbative MP2 total energy calculations, as well as on the expectation values of the bare exchange operator over the occupied and virtual canonical Hartree-Fock molecular orbitals. All calculations have been performed with the cc-pVTZ AO basis set. 30–33 Our different approximations are tested on the so-called Thiel’s set of molecules that contains a large set of small to medium size molecules, including unsaturated aliphatic hydrocarbons, aldehydes, ketones, amides and nucleobases. 34,38–42 Important technical details concerning the normalization of the auxiliary basis and the inversion threshold on the singular values of the pseudo inverse Γ−1 can be found in 43. In order to simplify the reading of the following sections, we will outline three specific procedures: (a) the standard RI-V based on cc-pVTZ-RI CAB that serves as a reference; (b) the constrained-momenta (up to fourth order) RI-SVS(µ4 ) approach with augmented aug-cc-pVTZRI CAB, and (c) the hybrid RI-V0 scheme with augmented aug-cc-pVTZ-RI CAB. As discussed in the following, the choice of the augmentation in the two later cases turns out to be important in order to provide an excellent accuracy, while the computational cost of both method remains favorable as compared to the one of the reference RI-V scheme (see Figure 1).

3.1

Self-Consistent Hartree-Fock Calculation

Statistics over total energies error for self-consistent Hartree Fock calculations with respect to the exact results are provided in Table 1, for both unconstrained and constrained fitting schemes combined with either SAB or CAB sets. We plot further in Figure 3 the error obtained for individual molecule (labeled from 1 to 28 following Thiel’s papers ordering). To clearly evidence the interest of the constrained schemes, we first emphasize that when left



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

Figure 3: Hartree-Fock total energy absolute errors per atom (excluding H atoms) for the 28 molecules in Thiel’s set. Data for calculations performed with the cc-pVTZ-RI Cartesian auxiliary set, with augmentation when indicated. unconstrained, the RI-SVS scheme is either unable to converge properly or leads to errors that are beyond the mHa per atom, even in the most favorable case of augmented CAB set. It is indeed known that robust fitting schemes 23–26,28 should be used when performing self-consistent calculations within the RI-SVS approach, but such schemes imply also the knowledge of the Coulomb three-center integrals, loosing thus the advantage of the simple RI-SVS approximation. The largest errors are indeed pathological of the self-consistent cycle: in particular, an analysis of the total fitted density shows clearly that it tends to move toward regions were it gets truncated, leading to a severe underestimation of the electron-electron repulsion energy. In that sense, both the augmentation of the auxiliary basis and the use of Cartesian orbitals instead of spherical ones have a strong impact on the RI-SVS fitting accuracy. Adding constraints on the RI-SVS fit improves drastically the situation, dropping the meanabsolute-error (MAE) down to few tenth of mHa per atom in the case of cc-pvtz-RI CAB. The use of an augmented auxiliary basis improves slightly the results, and both cases provide MAE in the same range as the reference RI-V one. The detailed analysis of the total energy error for each molecule of the set presented in Figure 3 confirms indeed that the accuracy of the constrained RI-SVS approach with augmented CAB is equivalent to the one of the reference RI-V scheme, but


Page 18 of 33

Page 19 of 33

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


Table 1: Self-consistent Hartree-Fock total energy mean absolute errors (MAE) in mHa/atom (excluding H atoms) with respect to exact calculations performed on Thiel’s set, as a function of the auxiliary basis set used. The parenthesis indicates the maximum order of the constraint that has been applied if any. µk indicates that momentum constraints up to nx + ny + nz = k have been used. The results of the three selected techniques have been highlighted. RI-SVS

RI-V0

RI-V

cc-pvtz-RI Sph.

>1Ha 0.879(µ4 )

0.319 0.107(µ4 )

0.022 0.085(µ4 )

cc-pvtz-RI Car.

>1Ha 0.074(µ4 )

0.061 0.015(µ4 )

0.040 0.026(µ4 )

aug-cc-pvtz-RI Sph.

26.78 0.233(µ4 )

0.044 0.031(µ4 )

0.016 0.018(µ4 )

aug-cc-pvtz-RI Car.

4.999 0.029(µ4 )

0.017 0.015(µ4 )

0.036 0.036(µ4 )

at a much lesser computational cost. Interestingly, the unconstrained RI-V0 scheme precision is already under the mHa per atom. As such, it can be already used as it is for most applications. When using simple cc-pVTZ-RI SAB and CAB, the accuracy can still be slightly improved through the application of constraints. All in all the RI-V0 scheme accuracy is similar to that of the RI-V one in the case of either CAB or augmented set. On the other hand, concerning the standard RI-V scheme, we can see that neither the use of SAB or CAB, nor the use of constraints, are able to significantly improve the already excellent accuracy. This confirms that the Coulomb fitting scheme efficiently preserves the loworder momenta of the exact (co)densities. In the RI-V0 case, where only the (l=0) auxiliary basis orbitals are combined with the Coulomb metric, the use of higher order constraints still appears to be beneficial, even though marginally as compared to the RI-SVS scheme.



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

3.2

MP2 total energy

Figure 4: MP2 total energy absolute errors per atom (excluding H atoms) for the 28 molecules in Thiel’s set. Data for calculations performed with the cc-pVTZ-RI Cartesian auxiliary set, with augmentation when indicated. We now explore the case of correlated MP2 calculations. Mean absolute errors (MAE) as compared to exact MP2 calculations are compiled in Table 2. The detail of the unconstrained approach errors for each molecule of the set are provided in Figure 4. As expected, the largest error for unconstrained fits is obtained with the RI-SVS scheme combined with the SAB, but interestingly the corresponding value can already be considered small, of the order of a few tenths of mHa per atom. Further, the RI-SVS scheme can be improved easily to excellent accuracy, either by imposing low order constraints, or by slightly increasing the auxiliary basis size, either switching from a spherical to a Cartesian representation, or by using augmentation. With such an accuracy of the simple RI-SVS scheme, the RI-V0 approach does not bring significant improvement. The effect of (co)density momenta enforcement is more subtle than for the previous HartreeFock calculations: for cc-pVTZ-RI auxiliary basis, it may bring an improvement at low order, but for higher µl values, the fit degrades and the error may actually increases. We interpret this outcome as the consequence of imposing constraints of too high order with respect to the quality of the auxiliary basis that is used. The error arises then for (co)densities that are poorly represented 20 ACS Paragon Plus Environment

Page 20 of 33

Page 21 of 33

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


Table 2: MP2 total energy mean absolute errors in mHa/atom (not counting H atoms) with respect to exact calculations performed on Thiel’s set, as a function of the auxiliary basis set used. The parenthesis indicates the maximum order of the constraint that has been applied if any. µk indicates that momentum constraints up to nx + ny + nz = k have been used. The results associated with the three selected techniques have been highlighted. RI-SVS

RI-V0

RI-V

cc-pvtz-RI Sph.

0.143 0.049(µ2 ) 0.061(µ3 ) 0.107(µ4 )

0.081 0.056(µ2 ) 0.114(µ3 ) 0.247(µ4 )

0.023 0.016(µ1 ) 0.102(µ2 ) 0.182(µ3 )

cc-pvtz-RI Car.

0.013 0.014(µ3 ) 0.023(µ4 ) 0.239(µ5 )

0.013 0.011(µ3 ) 0.031(µ4 ) 0.432(µ5 )

0.008 0.006(µ2 ) 0.012(µ3 ) 0.136(µ4 )

aug-cc-pvtz-RI Sph.

0.022 0.005(µ3 ) 0.004(µ4 ) 0.004(µ5 )

0.012 0.004(µ3 ) 0.005(µ4 ) 0.004(µ5 )

0.015 0.015(µ3 ) 0.014(µ4 ) 0.013(µ5 )

aug-cc-pvtz-RI Car.

0.003 0.002(µ3 ) 0.002(µ4 ) 0.002(µ5 )

0.007 0.003(µ3 ) 0.003(µ4 ) 0.003(µ5 )

0.007 0.007(µ3 ) 0.007(µ4 ) 0.007(µ5 )

within the auxiliary basis, typically induced by the products of canonical virtual orbital that appear within the MP2 formalism. In particular, improving the quality of the auxiliary basis, i.e., going from SAB to CAB, pushes further the maximum momenta accessible before degradation of the accuracy. Using augmented ABF sets completely cures this degradation, at least up to the maximum momenta tested here. This interpretation is also corroborated by the difference in accuracy between occupied and virtual canonical orbital exchange energies calculation, as reported in the following section. These results show that the RI-SVS scheme with augmented ABF sets can be considered a 21 ACS Paragon Plus Environment


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

Page 22 of 33

scheme of choice for MP2 total energy calculations, possibly combined with low order constraints, since its accuracy compares well to the RI-V one while retaining an O(N) complexity. Increasing the quality of the auxiliary basis set (i.e., the number of basis functions) within the RI-SVS scheme appears as a good strategy to achieve better accuracy while keeping reasonable computational costs.

3.3

Bare exchange contribution to molecular orbital energies

We explore now the accuracy of the various RI implementations in the calculation of the exactexchange contribution h ψn | ΣX | ψn i to the Hartree-Fock canonical orbital eigenvalues, namely: X

h ψn | Σ | ψn i =

Z

∑

drdr0

i ∈ occ.

ψn∗ (r)ψi∗ (r)ψn (r0 )ψi (r0 ) |r − r0 |

(24)

where the ψn/i are the exact Hartree-Fock canonical orbitals as provided by the NWChem package. Such a quantity is typical of perturbative correlated techniques, such as the GW formalism, 46 aiming at calculating accurate molecular orbital (HOMO, LUMO, etc.) energies. 47,48 Such techniques require the calculation of the expectation value of an exchange-correlation non-local self-energy ΣXC operator, expressed, e.g., in terms of the bare and screened Coulomb operators, over KohnSham or Hartree-Fock eigenstates. The statistics presented account for all molecules in the set of Thiel. The MAEs are first calculated over the manifold of all valence states (Table 3), and secondly over an equal number of empty states, taking the lowest ones (Table 4). Occupied states - As shown in Table 3, the RI-SVS scheme in combination with SAB leads to the largest error (∼25 meV). Again, the accuracy can be significantly improved by imposing constraints, reducing the error to less than a meV at the fourth-order (µ4 ) level. Further improvements can be obtained by increasing the auxiliary basis size, in particular by using augmentation. The favorable scaling in number of three-center integrals with auxiliary basis size for the density-fitting scheme clearly suggests that the constrained RI-SVS scheme with augmented basis is the method


Page 23 of 33

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


Table 3: Mean absolute errors (meV) for the bare exchange contribution to occupied molecular orbital energies. The parenthesis indicates the maximum order of the constraint that has been applied, if any. The results associated with the three selected techniques have been highlighted. RI-SVS

RI-V0

RI-V

cc-pvtz-RI Sph.

23.12 1.920(µ3 ) 0.932(µ4 )

4.225 1.192(µ3 ) 1.139(µ4 )

0.351 0.329(µ3 ) 0.754(µ4 )

cc-pvtz-RI Car.

3.224 0.341(µ3 ) 0.377(µ4 )

0.534 0.308(µ3 ) 0.360(µ4 )

0.213 0.186(µ3 ) 0.274(µ4 )

aug-cc-pvtz-RI Sph.

8.289 0.229(µ4 ) 0.219(µ5 )

0.308 0.198(µ4 ) 0.235(µ5 )

0.253 0.253(µ4 ) 0.254(µ5 )

aug-cc-pvtz-RI Car.

1.894 0.155(µ4 ) 0.149(µ5 )

0.159 0.151(µ4 ) 0.147(µ5 )

0.186 0.187(µ4 ) 0.187(µ5 )

of choice, bringing an accuracy comparable to the RI-V scheme with much reduced computer time. In particular, while RI-V0 offers a clear improvement over the standard unconstrained RI-SVS approach, the use of constraints dramatically reduces the differences between the two techniques introduced in the present study. Unoccupied states - Table 4 shows that, in the case of cc-pVTZ-RI ABF, error statistics on bare exchange orbital energies obtained for unoccupied states are much larger than their occupied states analogs, even in the case of the RI-V scheme. While the RI-V MAE remains under 10 meV independently of the use of SAB or CAB, only the constrained RI-SVS and RI-V0 MAE drop under 10 meV. The behavior of the error as a function of the applied constraint maximum order is similar to the one described for MP2 calculations. Again, the maximum momenta accessible before degradation of the fit increases when going from SAB to CAB.



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

Page 24 of 33

Table 4: Mean absolute errors (meV) for the bare exchange contribution to unoccupied molecular orbital energies. The parenthesis indicates the maximum order of the constraint that has been applied, if any. The results associated with the three selected techniques have been highlighted. RI-SVS

RI-V0

RI-V

cc-pvtz-RI Sph.

25.54 9.734(µ2 ) 12.72(µ3 )

28.09 14.16(µ2 ) 26.19(µ3 )

2.286 2.083(µ1 ) 15.29(µ2 )

cc-pvtz-RI Car.

13.44 5.677(µ3 ) 12.80(µ4 )

10.53 6.134(µ3 ) 15.66(µ4 )

1.175 3.044(µ2 ) 12.80(µ3 )

aug-cc-pvtz-RI Sph.

0.584 0.181(µ4 ) 0.171(µ5 )

0.456 0.187(µ4 ) 0.211(µ5 )

0.225 0.194(µ4 ) 0.204(µ5 )

aug-cc-pvtz-RI Car.

0.223 0.116(µ4 ) 0.112(µ5 )

0.174 0.120(µ4 ) 0.120(µ5 )

0.161 0.158(µ4 ) 0.153(µ5 )

The use of the augmented ABF turns out to be crucial int the case of virtual canonical orbitals, even for the standard RI-V scheme. To start with, it brings back the three method results on the occupied results footing. Also, in this case, even the unconstrained RI-SVS fit behaves accurately. We interpret this as the consequence of fitting neutral (co)densities, i.e., with no monopole, as implied with virtual canonical orbitals exchange energies calculation. Due to this fact, the virtual canonical orbital case is the one that benefits the less from the constrained approach. As in the MP2 case, the use of augmented ABF also cures the degradation of the fit that arises with high order momenta constraints. Finally, a good strategy for perturbative correlated techniques seems thus to consider augmented ABF sets, possibly used in conjunction with the constrained RI-SVS approach, in order to keep a low computational cost. This strategy provides indeed sub-meV accuracy while retaining


Page 25 of 33

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


overall an advantageous linear scaling of three-center overlap integral number.

4

Conclusions

We have explored novel resolution-of-identity (RI) strategies that provide similar accuracy, but with a much reduced number of needed integrals, as compared to the standard Coulomb-fitting (RI-V) approach. We have presented in particular a dual space formalism allowing to combine the density (RI-SVS), Coulomb (RI-V) and constrained fitting schemes, yielding flexible generalized approaches easy to implement in quantum chemistry codes. We also demonstrated that the originals Baerends and Van Alsenoy constrained density-fitting approaches can be significantly improved and stabilized. The schemes we introduced were tested by calculating the Hartree-Fock and MP2 total energies over a standard molecular set of medium size molecules, together with the corresponding bare exchange contribution to the Hartree-Fock occupied and virtual molecular orbitals. Our generalized constrained density-fitting RI-SVS formalism is shown to lead to significant improvements over the standard (unconstrained) density-fitting scheme, in particular used in conjunction with augmented auxiliary basis sets. The dramatically smaller number of needed integrals, as compared to the RI-V approach, and the favorable scaling with auxiliary basis size, leads to a situation where the constrained RI-SVS scheme with augmented auxiliary bases is as accurate as the standard RI-V approach, but with a dramatically reduced number of three-center integrals to be stored or re-calculated. Alternatively, an hybrid RI-V0 scheme, namely imposing the Coulomb metric for (l=0) auxiliary orbitals, and the density-fitting RI-SVS scheme for auxiliary basis orbitals with larger momentum, can be shown to provide an accuracy similar to the full RI-V approach upon slightly increasing the auxiliary basis size, either by going from spherical to Cartesian basis functions, or by including diffuse orbitals. While the RI-V0 scheme offers the same O(N 2 ) scaling of the number of three-center integrals as the full RI-V technique, the much reduced prefactor leads to



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

dramatically reduced number of three-centers integrals even when increasing the auxiliary basis size through the addition of diffuse functions. In all cases, the use of hybrid and/or constrained techniques with slightly increased auxiliary bases allows to obtain an accuracy similar to the standard RI-V Coulomb-fitting strategy but with a much reduced computational cost.

Acknowledgement The authors acknowledge Thierry Deutsch and Denis Jacquemin for a critical reading of the manuscript. The author are also indebted to Antoine Levitt for helpful comments and discussions. This project has received funding from the European Union Horizon 2020 research and innovation program under Grant agreement No 646176 (EXTMOS).

Supporting Information Available Supplementary computational details such as normalization schemes and corresponding pseudoinverse thresholds. This material is available free of charge via the Internet at http://pubs. acs.org.

References (1) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; McGraw-Hill Publishing Company: New York, 2000. (2) Helgaker, T.; Jorgensen, P.; Olsen, J. Molecular Electronic-Structure Theory; John Wiley and Sons, 1989. (3) Häser, M.; Ahlrichs, R. J. Comput. Chem. 1989, 10, 104–111. (4) Whitten, J. L. J. Chem. Phys. 1973, 58, 4496–4501. (5) Baerends, E.; Ellis, D.; Ros, P. Chem. Phys. 1973, 2, 41 – 51. 26 ACS Paragon Plus Environment

Page 26 of 33

Page 27 of 33

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


(6) Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. J. Chem. Phys. 1979, 71, 3396–3402. (7) Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. J. Chem. Phys. 1979, 71, 4993–4999. (8) Mintmire, J. W.; Dunlap, B. I. Phys. Rev. A 1982, 25, 88–95. (9) Van Alsenoy, C. J. Comput. Chem. 1988, 9, 620–626. (10) Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Chem. Phys. Lett. 1993, 208, 359 – 363. (11) Vahtras, O.; Almlof, J.; Feyereisen, M. Chem. Phys. Lett. 1993, 213, 514 – 518. (12) Klopper, W.; Samson, C. C. M. J. Chem. Phys. 2002, 116, 6397–6410. (13) Eichkorn, K.; Treutler, O.; Öhm, H.; H"aser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 240, 283 – 290. (14) Weigend, F.; H´’aser, M.; Patzelt, H.; Ahlrichs, R. Chem. Phys. Lett. 1998, 294, 143 – 152. (15) Weigend, F. Phys. Chem. Chem. Phys. 2002, 4, 4285–4291. (16) Weigend, F. Phys. Chem. Chem. Phys. 2006, 8, 1057–1065. (17) Hellweg, A.; Hättig, C.; Höfener, S.; Klopper, W. Theor. Chem. Acc. 2007, 117, 587–597. (18) Beebe, N. H. F.; Linderberg, J. Int. J. Quant. Chem. 1977, 12, 683–705. (19) Weigend, F.; Kattannek, M.; Ahlrichs, R. J. Chem. Phys. 2009, 130, . (20) Bernholdt, D. E.; Harrison, R. J. Chem. Phys. Lett. 1996, 250, 477 – 484. (21) Weigend, F.; Häser, M. Theor. Chem. Acc. 1997, 97, 331–340. (22) Ihrig, A. C.; Wieferink, J.; Zhang, I. Y.; Ropo, M.; Ren, X.; Rinke, P.; Scheffler, M.; Blum, V. New J. Phys. 2015, 17, 093020. (23) Dunlap, B. J. Mol. Struct.: {THEOCHEM} 2000, 501 - 502, 221 – 228.



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

(24) Dunlap, B. J. Mol. Struc: {THEOCHEM} 2000, 529, 37 – 40. (25) Reine, S.; Tellgren, E.; Krapp, A.; Kjaergaard, T.; Helgaker, T.; Jansik, B.; Höst, S.; Salek, P. J. Chem. Phys. 2008, 129, . (26) Domínguez-Soria, V. D.; Geudtner, G.; Morales, J. L.; Calaminici, P.; Köster, A. M. J. Chem. Phys. 2009, 131, . (27) Dunlap, B. I.; Rösch, N.; Trickey, S. Mol. Phys. 2010, 108, 3167–3180. (28) Mejía-Rodríguez, D.; Köster, A. M. J. Chem. Phys. 2014, 141, . (29) Ren, X.; Rinke, P.; Blum, V.; Wieferink, J.; Tkatchenko, A.; Sanfilippo, A.; Reuter, K.; Scheffler, M. New J. Phys. 2012, 14, 053020. (30) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007–1023. (31) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796–6806. (32) Feller, D. J. Comput. Chem. 1996, 17, 1571–1586. (33) Schuchardt, K. L.; Didier, B. T.; Elsethagen, T.; Sun, L.; Gurumoorthi, V.; Chase, J.; Li, J.; Windus, T. L. J. Chem. Inf. Model. 2007, 47, 1045–1052. (34) Schreiber, M.; Silva-Junior, M. R.; Sauer, S. P. A.; Thiel, W. J. Chem. Phys. 2008, 128, 134110. (35) Gill, P. M. W. The Journal of Physical Chemistry 1996, 100, 15421–15427. (36) Jung, Y.; Sodt, A.; Gill, P. M. W.; Head-Gordon, M. Proc. Natl. Acad. Sci. USA 2005, 102, 6692–6697. (37) Genovese, L.; Deutsch, T. Phys. Chem. Chem. Phys. 2015, 17, 31582–31591. (38) Silva-Junior, M. R.; Schreiber, M.; Sauer, S. P. A.; Thiel, W. J. Chem. Phys. 2008, 129, 104103. 28 ACS Paragon Plus Environment

Page 28 of 33

Page 29 of 33

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


(39) Sauer, S. P. A.; Schreiber, M.; Silva-Junior, M. R.; Thiel, W. J. Chem. Theory Comput. 2009, 5, 555–564. (40) Silva-Junior, M. R.; Thiel, W. J. Chem. Theory Comput. 2010, 6, 1546–1564. (41) Silva-Junior, M. R.; Sauer, S. P.; Schreiber, M.; Thiel, W. Mol. Phys. 2010, 108, 453–465. (42) Silva-Junior, M. R.; Schreiber, M.; Sauer, S. P. A.; Thiel, W. J. Chem. Phys. 2010, 133, 174318. (43) Computation details on auxiliary and dual basis normalization, as well as pseudo inverse threshold can be found in the Supporting Information. (44) It is striking that odd constraint orders do not bring much accuracy. This last fact however does not hold when changing the auxiliary basis normalization scheme. (45) Valiev, M.; Bylaska, E.; Govind, N.; Kowalski, K.; Straatsma, T.; Dam, H. V.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T.; de Jong, W. Comput. Phys. Comm. 2010, 181, 1477 – 1489. (46) Hedin, L. Phys. Rev. 1965, 139, A796–A823. (47) Blase, X.; Attaccalite, C.; Olevano, V. Phys. Rev. B 2011, 83, 115103. (48) Faber, C.; Attaccalite, C.; Olevano, V.; Runge, E.; Blase, X. Phys. Rev. B 2011, 83, 115123.



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

A

Page 30 of 33

Link to Van Alsenoy Lagrangian

In order to impose constraints on the fitted densities, Van Alsenoy 9 used the following Lagrangian, composed of the standard fitting Lagrangian completed with a Lagrange multiplier for each constraint: 1 2

2 dr ρ(r) − ∑ cβ β (r) − ∑ λµ h µ | ρ − ∑ cβ β i

Z

χ

β

µ

(25)

β

Here χ(r, r0 ) denotes the chosen metric, i.e. χ(r, r0 ) = δ (r, r0 ) for RI-SVS and χ(r, r0 ) = v(r, r0 ) for the RI-V scheme. The derivation condition with respect to cβ coefficient leads to:

∑0 cβ 0 h β | χ | β 0 i − h β | χ | ρ i + ∑ λµ h µ | β i = 0

(26)

µ

β

The derivation condition with respect to λµ coefficients simply gives:

h µ | ρ − ∑ cβ β i

(27)

β

For simplicity, we denotes vβ the vector which coefficients are h β | χ | ρ i. In the same manner, the vector of coefficients h µ | ρ i is written vµ , and Λ denotes the vector of coefficient λµ . Again, we write Γ the matrix which coefficients are h β | χ | β 0 i and M the matrix formed by h µ | β i coefficients. The problem can be rewritten under matrix form as: vµ − M · c = 0

(28)

†

Γ · c − vβ + M · Λ = 0 Reordering the second line leads to: vµ − M · c = 0 Γ

−1

†

· (vβ − M · Λ) = c


(29)

Page 31 of 33

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


Combining with the first line brings then: vµ − M · Γ−1 · (vβ − M † · Λ) = 0 ⇒Λ = (M · Γ

−1

† −1

·M )

· (M · Γ

−1

(30)

· vβ − vµ )

Finally, writing T = M · Γ−1 and substituting Λ leads to: c = Γ−1 · (vβ − M † · (T · M † )−1 · (T · vβ − vµ )) = Γ−1 · (I − M † · (T · M † )−1 · T ) · vβ

(31)

+ Γ−1 · M † · (T · M † )−1 · vµ By identification with components of Eq. (17) and considering that in this particular case Γ = Γ† , we can deduce that: A = Γ−1 · (I − M † · (T · M † )−1 · T ) (32) † −1

†

= (I − T · (M · T )

· M) · Γ

−1

B = Γ−1 · M † · (T · M † )−1 (33) † −1

†

= T · (M · T )

We fall back on the generic expression of A and B given in Eq. (20) by choosing to write T = 1

M · Γ−1 = M · (Γ† Γ)− 2 . Van Alsenoy expression corresponds thus to the choice α = 21 .

B α = 1 as optimal fitting criterion In this Appendix, we demonstrate that the α = 1 choice corresponds to a wise choice with respect to a minimal disturbance of the dual measures h γ | ρ i, without further knowledge over the set of test densities ρ. Let’s start by recalling the generic expression of the constrained fitting operator



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

Page 32 of 33

Iβ given by Eq. (17): Iβ = ∑[A]β γ | β ih γ | + ∑[B]β µ | β ih µ | βγ

βµ

with A and B generic expressions of Eq. (20): B = T † (MT † )−1 A = I − T † (MT † )−1 M Γ−1 Writing vρ the vector of h γ | ρ i coefficients and wρ the one with h µ | ρ i coefficients, we can express the fact that we request a minimal perturbation of the dual measures of a particular density ρ as: 0=

d dT

2 h γ | ρ i − h γ | I | ρ i ∑ β γ

2 d (I − ΓA)v − ΓBw ρ ρ dT 2 d † = Γ T (MT † )−1 MΓ−1 vρ − wρ dT

=

(34)

Now remains just to verify the property (34) for an ensemble {ρ} of test densities. Here, one could choose the Kohn-Sham or Hartree-Fock molecular orbitals product basis, or its restriction to a subset pertinent in terms of the (co)densities that we want to fit. However without considerations on the molecular orbitals, and thus no further knowledge on the test set, we assume here that property (34) is valid for any test density. Writing for clarity cρ = MΓ−1 vρ − wρ , expression


Page 33 of 33

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


(34) can be further developed as: 2 d † † −1 0= Γ T (MT ) cρ dT d = tr (T M † )−1 T Γ† Γ T † (MT † )−1 cρ c†ρ dT † −1 † † −1 † † † −1 = 2 (MT ) cρ cρ (T M ) T Γ Γ I − T (MT ) M

(35)

The condition (35) is then fulfilled for all cρ when: T Γ† Γ I − T † (MT † )−1 M = 0 ⇒T Γ† Γ = M

(36)

⇒T = M( Γ† Γ)−1 We see that the choice of T = M( Γ† Γ)−α with α = 1 corresponds to the choice that minimizes the correction brought to the dual measures of the (co)densities, in absence of any consideration on the (co)density being measured.


Hybrid and Constrained Resolution-of-Identity ... - ACS Publications

Recommend Documents