Laplacetransformed multireference secondorder perturbation theories in the atomic and active molecular orbital basis
Abstract
In the present article, we show how to formulate the partially contracted nelectron valence second order perturbation theory (NEVPT2) energies in the atomic and active molecular orbital basis by employing the Laplace transformation of orbitalenergy denominators (OED). As atomicorbital (AO) basis functions are inherently localized and the number of active orbitals is comparatively small, our formulation is particularly suited for a linearlyscaling NEVPT2 implementation. Some of the NEVPT2 energy contributions can be formulated completely in the AO basis as singlereference secondorder Møller–Plesset perturbation theory and benefit from sparse activepseudo density matrices — particularly if the active molecular orbitals are localized only in parts of a molecule. Furthermore, we show that for multireference perturbation theories it is particularly challenging to find optimal parameters of the numerical Laplace transformation as the fit range may vary among the 8 different OEDs by many orders of magnitude. Selecting the number of quadrature points for each OED separately according to an accuracybased criterion allows us to control the errors in the NEVPT2 energies reliably.
I Introduction
Many applications of quantum chemical simulations require a multireference (MR) wave function to provide an at least qualitatively correct description of the atomic or molecular system. Even the simplest molecule — the hydrogen molecule — becomes a MR case in the process of separating the two hydrogen atoms. MR theories are also indispensable when two potential energy surfaces of the same symmetry come close in energy, which is frequently encountered when studying the time evolution of excited electrons in molecules — in particular by surfacehopping excitedstate dynamics.Martínez (2006); *Subotnik2016 Beyond a reliable treatment of bond breaking and potential energy surfaces, MR theories are often considered in quantum chemical studies that involve transition metals, lanthanides, and actinides to cope with open shells in the electronic structure calculation. Alternatively, KohnSham density functional theory (DFT) can often provide meaningful energies and properties for transition metal complexes with low computational costs even when standard exchange correlation functionals are employed. However, DFT still lacks an ansatz that can improve systematically on either static or dynamic correlation effects in a computationally affordable manner.
In comparison to their singlereference (SR) analogues, wave functionbased MR electronic structure methods are conceptually and computationally much more demanding. At the multiconfigurational (MC) selfconsistent field (SCF) level, the energy needs to be minimized with respect to orbital rotations and configuration coefficients, which in many situations requires sophisticated optimizers with quadraticLengsfield (1980) or higherWerner and Meyer (1980); *Werner1981; *Werner1985 convergence rates. On the one hand, such algorithms can be incredibly efficient in terms of number of iterations; on the other hand, a transformation of the twoelectron integrals from the atomic (AO) to the molecular orbital (MO) basis is necessary, which scales with the system size . Moreover, the MCSCF wave function is expanded in configuration state functions (CSF) that can be generated by a complete active spaceRoos, Taylor, and Siegbahn (1980) (CAS). The latter includes all possible occupations of active MOs by active electrons usually denoted by CAS (,). With conventional determinantbased full configuration interaction (FCI) implementationsKnowles and Handy (1984); Olsen, Jørgensen, and Simons (1990) the computational costs grow exponentially and becomes soon a bottleneck in the MCSCF calculation. Usually, one is restricted to active spaces in the ballpark of a CAS (10,30) if no restrictedOlsen et al. (1988) (RAS) or generalized active spacesFleig, Olsen, and Marian (2001) (GAS) are desired.
There are recent developments in MCSCF that improve on both computational bottlenecks of the algorithms: the twoelectron integral transformation and the size limit of the active space. The latter can be lifted if the CAS CI secular equations are solved by modern quantum Monte CarloBooth, Thom, and Alavi (2009); Thomas et al. (2015) or density matrix renormalization groupWhite (1992); *whit93 (DMRG) algorithms. This facilitates CASSCF calculations with much larger active spaces, e. g. .Zgid and Nooijen (2008a); *zgid08b; *ghos08; *yana09a; *kurashige2011; *wouters2014; *wout15; *mayi16c; *LiManni2016; Knecht et al. (2016); Sun, Yang, and Chan (2017) As shown lately by Hohenstein et al.Hohenstein et al. (2015) for an approximate secondorder optimizerChaban, Schmidt, and Gordon (1997), the costly integral transformations can be avoided by initially computing Coulomb and exchangelike matrices in the sparse AO basis with costs and later transforming to Focklike intermediates and orbital gradients in the MO basis. Just recently, a truly secondorder optimizer that is also AOdriven and can handle a large number of AOs and active MOs was proposed in Ref. 26
The computationally simplest approach to account for dynamic correlation in molecules that demand for multiple reference determinants is secondorder multireference perturbation theory (MRPT2). MRPT2 methods made substantial progress over the last 25 years and are frequently used in quantum chemical investigations. In particular, the CASPT2 methodAndersson et al. (1990); *Andersson1992 initially developed by Andersson et al. has evolved into a reliable tool for the investiagation of molecules with a complicated electronic structure. Geometries can be conveniently optimized due to availability of analytic CASPT2 gradients,Celani and Werner (2003); *Shiozaki2011; *Gyoerffy2013 which were developed by Werner and his coworkers. Accurate CASPT2 energies and gradients near the complete basis set limit can be obtained with an explicitly correlated F12 variant proposed by Shiozaki et al.Shiozaki and Werner (2010); *Shiozaki2013 For calculations involving heavy elements, spinorbit coupling (SOC) can be treated either within a SOCCI formalismRoos and Malmqvist (2004) or by using the fourcomponent DiracCoulomb Hamiltonian.Abe, Nakajima, and Hirao (2006) Concerning studies of excited states, the erratic mixing of valence and Rydberg states at the CASSCF level can be cured by using quasidegenerate MRPT2.Finley et al. (1998); Angeli et al. (2004)
Despite its popularity and success in numerous application, CASPT2 has two major drawbacks: (i) It is not sizeconsistent; (ii) It is often plagued by intruder states, which complicates the iterative determination of the firstorder wavefunction coefficients. In practice, the sizeinconsistency errors are rather decent — even for larger moleculesMenezes, Kats, and Werner (2016) — and the intruderstate problem can be cured partially by level shifts.Roos and Andersson (1995); *Forsberg1997 However, since level shifts alter the secondorder energy unpredictably, CASPT2 should then be considered as a rather empirical method.
The remedies of CASPT2 are cured altogether within a different formulation of MRPT: the nelectron valence (NEV) PT of Angeli, Cimiraglia, Malrieu and their coworkers.Angeli, Cimiraglia, and Malrieu (2000); Angeli et al. (2001); Angeli, Cimiraglia, and Malrieu (2001, 2002) NEVPT2 differs from CASPT2 mainly in the definition in the zerothorder Hamiltonian. In NEVPT the Fock operator is replaced by the Dyall HamiltonianDyall (1995), which also accounts for explicit twoelectron interactions within the active orbitals. Initially, two variants of NEVPT were introduced that differ in the level of socalled outerspace contraction:Angeli et al. (2001) the partially contracted (PC) and strongly contracted (SC) NEVPT2. Both variants are fully internally contracted (IC) and require the computation of the fourparticle reduced density matrix (RDM), which scales with . With the rise of DMRG and FCIQMC algorithms in the last two decades, relatively large active spaces are now feasible for CASSCF, but become prohibitive in the subsequent ICMRPT calculation. A rigorous approach to this dilemma would be to use completely uncontracted firstorder wavefunctions, which was initially considered as computationally not feasible. Recent developmentsSharma and Chan (2014); Sokolov and Chan (2016); Giner et al. (2017) reveal, however, the clear potential of fully uncontracted MRPT for calculations that demand large active spaces.
Most systems with MR character that have a practical relevance, e. g. metalorganic complexes, enzymes, or organic radicals, have a rather small number of active orbitals compared to the total number of MOs. In such situations, the computational costs of MRPT2 are dominated by the scaling integral transformations as for SR secondorder Møller–Plesset PT (MP2). A substantial reduction of the timings can be achieved by employing a CD of twoelectron integralsAquilante et al. (2008); *Freitag2017, in particular, if such an integral decomposition technique is combined with a truncated or frozen natural orbital (FNO) virtuals basis.Aquilante et al. (2009) However, for FNO or CDMRPT2 only the prefactor of the computational costs is reduced but the scaling is still as for conventional implementations. Thus, calculation that exploit those approximations are still limited to mediumsized molecules in the ballpark of 100 atoms. Much larger systems can be calculated if the firstorder wave function is expanded in a truncated set of pair natural orbitals (PNO), which was explored recently by the NeeseGuo et al. (2016) and Werner groupsMenezes, Kats, and Werner (2016) for PNONEVPT2 and PNOCASPT2, respectively. Both implementations scale linearly with the system size and allow to perform MRPT2 calculations of molecules with up to 300 atoms.
In the present article, we pursue an alternative approach to reach linearly scaling MRPT2 implementations eventually. It is based on the Laplace transformation of orbitalenergy denominators,Almlöf (1991) which was introduced initially by Almlöf in the context of MP PT.Häser and Almlöf (1992); Häser (1993) In our formulation, all intermediates are given solely in the AO and active MO basis, which is in line with recent algorithmic developments in MCSCFHohenstein et al. (2015); Sun, Yang, and Chan (2017). On the one hand, the number of active MOs can still be considered as small in MR calculations because even with modern DMRG and FCIQMC implementations solving the CAS CI secular equations may easily become a computational bottleneck due to the formal exponential cost scaling with . On the other hand, localized AO basis functions are inherently suited for effective screening. It was shown by Ayala and ScuseriaAyala and Scuseria (1999) that SR Laplacetransformed AOMP2 energies show the physically correct decay with the distance between charge distributions. This was exploited later by Ochsenfeld and his coworkers,Lambrecht, Doser, and Ochsenfeld (2005); Doser et al. (2009); *Maurer2013b who introduced integral estimates with the proper dependence to develop linearly scaling MP2 implementations. With those AObased implementations MP2 calculations on molecules with more than 1000 atoms and more than 10000 basis functions were reported.Doser et al. (2009) However, since the initial AOMP2 implementation of HäserHäser (1993), AObased wave function methods were criticized for their large computational overhead that increases substantially for larger basis sets. This problems can be solved now by a Cholesky decomposition (CD) of (pseudo)densitiesZienau et al. (2009), which leads to intermediates that have the same dimension as the occupied or virtual MO coefficients but retain the sparsity of AO (pseudo)density matrices.Kussmann et al. (2015); Luenser, Schurkus, and Ochsenfeld () An extension of Laplacetransformed AObased MP2 for relativistic allelectron calculations was introduced recently by one of the present authors for twocomponent Hamiltonians and the Kramersrestricted formalism.HelmichParis, Repisky, and Visscher (2016)
In this work, we show how to formulate NEVPT2 in the atomic and active molecular orbital basis by means of the numerical Laplace transformation. We start in Sec. II with a résúme of present stateoftheart ICMRPT2 theories and argue why only the partiallycontracted NEVPT2 is suitable for such a reformulation. Moreover, we outline how to obtain PCNEVPT2 energies in the AO and active MO basis in the spinfree formalism and provide working equations for each energy contribution. In Sec. IV, the working equations are discussed in context of a potentially efficient implementation. Furthermore, we focus on how to control the error of the numerical quadrature of the Laplace transformation for each of the individual energy contributions. Finally, our procedure is validated for some typical applications of MRPT methods that involve bond dissociation and excited valence and Rydberg states.
Ii Theory
ii.1 Résúme of CASPT2 and NEVPT2
The zerothorder wavefunction is given by a linear combination
(1) 
of all configuration state functions . The wavefunction parameters, that is orbital rotations and CI coefficients , can be optimized either in a MCSCF or openshell HF and subsequent CI calculation. The firstorder correction to the wave function
(2) 
is expanded through all possible determinants that are generated by twoelectron excitations from the complete zerothorder wavefunction. This approximation is termed internal contraction (IC) and is referred to as internally contracted configurationsCelani and Werner (2000) (ICC). To determine the PT energy through second order,
(3) 
the expansion coefficients of the firstorder wave function are determined by solving the linear system of equations
(4) 
In Eq. (II.1) we assume that diagonalization of the zerothorder Hamiltonian,
(5) 
gives the MCSCF or CI eigenvalues and zerothorder wave function . As projection manifold we choose bra ICCs in Eq. (II.1) such that they keep the core and virtual part of the ICC overlaps orthogonalPulay, Sæbø, and Meyer (1984); *Helgaker2000biortho because simplified arithmetic expressions are obtained.
To build , the mutually orthogonal ICCs are generated by exciting two electrons out of the core or active into the virtual or active molecular orbitals (MO) of the zerothorder wavefunction . There are 8 different classes of ICCs if the reference determinants were built from a complete active space (CAS), which are compiled in Tab. 3. In the present work, we use spinfree operators that are composed of singlet oneparticle excitation operators
(6) 
We follow the convention that core orbitals are denoted by , active orbitals by , virtual orbitals by , and general orbitals by . The different classes of ICCs in Tab. 3 are labeled by symbols that became standard notation in NEVPT: denotes attachment of electron(s) into the active space while means ionization of electron(s) from the active space; a prime designates a singleparticle excitation within the active orbital space. In contrast to other ICCs, the one of the class are composed of two kinds of excited determinants, and , that are mutually nonorthogonal. If RAS or GAS were employed in the underlying MCSCF calculation, the expansion of in terms of ICCs needs to be augmented with a subset of the CAS space determinants, i. e. , to keep the zeroth and firstorder wavefunction mutually orthogonal.Celani and Werner (2000) In the present work, we only consider CAS reference wavefunctions.
In MRPT the definition of the zerothorder Hamiltonian is not unique and it is the choice of in which CASPT2 and PCNEVPT2 mainly differ. In CASPT2 the zerothorder Hamiltonian
(7) 
is defined by projection of the Fock operator
(8)  
(9)  
(10)  
(11) 
onto the space of reference determinants and its complement .Celani and Werner (2000) The inactive and active Fock matrixHelgaker, Jørgensen, and Olsen (2000b) are built from the one and twoelectron integrals , as they appear in the electronic Hamiltonian
(12)  
(13)  
(14) 
and the singlet oneparticle reduced density matrix (1RDM) (see Tab. 2). The integrals in the AO basis and are transformed to the corresponding MO integrals and by MO coefficients that are available from the preceding MCSCF or openshell HF calculation. The CI coefficients appear only in the nRDMs, as given in Tab. 2, for IC PTs, which are our focus here.
Due to the nondiagonal block structure of in CASPT2, can couple different classes of ICCs in . This has two unpleasant implications: first, Eq. (II.1) must be determined iteratively and a larger computational overhead compared to a direct method like singlereference MP2 is expected; second, CASPT2 may not be size extensivevan Dam, van Lenthe, and Pulay (1998); *VanDam1999 what, in principle, disqualifies the method to be applied to large molecular systems. Furthermore, the usage of the Fock operator not only for core and virtual but also for the active orbitals may lead to (near) singularities in Eq. (II.1). In particular, ICC blocks in that involve three active orbitals, and , are most likely to become zero.Dyall (1995); Angeli, Cimiraglia, and Malrieu (2000) This is the root of the notorious CASPT2 intruder state problem, which can be cured either by including more orbitals into the active space or in an empirical fashion by introducing level shifts in the zerothorder HamiltonianRoos and Andersson (1995); *Forsberg1997. Furthermore, near singularities in would impede its Laplace transform as many quadrature points in the numerical integration procedure would be required then.
The NEVPT2 zerothorder Hamiltonian features only diagonal ICC blocks,Dyall (1995)
(15) 
which can guarantee size extensivity.van Dam, van Lenthe, and Pulay (1998); *VanDam1999 Furthermore, Dyall’s Hamiltonian is employed in Eq. (15) that chooses the core and valence part differently,
(16)  
(17)  
(18) 
and accounts for twoelectron interactions among the active electrons.Dyall (1995); Angeli, Cimiraglia, and Malrieu (2000, 2002) The constant is chosen such that Eq. (5) is fulfilled.
The bielectronic valence part of Dyall’s Hamiltonian in the commutator of Eq. (II.1) leads to the socalled Koopmans matrices. For example, the Koopmans matrix is given as
(19)  
(20) 
and represents a singleionization potential.Morrell, Parr, and Levy (1975); Smith and Day (1975); Angeli et al. (2001) In Eq. (20) is the singlet 2RDM, which is defined in Tab. 2. The term Koopmans matrices for classes denoted with a prime is chosen for reasons of notational convenience rather than physical rationality. Explicit expressions for all 7 kinds of Koopmans matrices are presented in the spinfree formalism in Ref. Angeli, Cimiraglia, and Malrieu, 2002 and not repeated here. For large active spaces the computation of the and Koopmans matrices can be demanding since it requires prior calculation of the 4RDM (see Tab. 2).
The iterative solution of the linear system of equations (II.1) can be avoided in the partially contracted (PC) NEVPT2 variant by solving a generalized eigenvalue problem (GEP),
(21) 
for each excitation class separately. Eq. (21) involves the Koopmans matrices and the active part of the ICC overlap , which is denoted here as metric matrix. The metric matrices for all excitation classes are compiled in Tab. 2 and formulated in terms of singlet RDMs. If a CAS zerothorder wave function is employed, both and are Hermitian and, thus, the eigenvalues are necessarily real for each of the 7 classes.Angeli, Cimiraglia, and Malrieu (2002) Furthermore, the Koopmans matrices are either positive definite for those classes representing active space electron attachment or excitation or negative definite for the active electron ionization classes. Consequently, in the firstorder equation (II.1), can never become singular, that is, there are no intruder states in NEVPT2. This facilitates a direct calculation of the PCNEVPT2 energy with a formulation that is suited for a Laplace transformation of denominators involving Fock and Koopmans matrix eigenvalues, as we elaborate on in Sec. II.2.
Apart from choosing a IC firstorder wave function in Eq. (2), as done in CASPT2 and PCNEVPT2, a further level of contraction can be introduced by choosing an alternative excitation operator basis for the ICCs. For the strongly contracted (SC) NEVPT2, the active orbitals of the two particle excitations operators are contracted with parts of the perturbation operatorAngeli et al. (2001) as shown in Tab. 3 for the bra ICCs. This leads to intermediates in denominators of the energy expressions that have the same dimension as the SC ICCs, e. g. one obtains an effective quasihole energy
(22) 
for the class. Such effective energies, as in Eq. (22), are highly unsuited for a Laplace transformation of the denominators since, unlike the MP denominators, they cannot be partitioned into individual orbital contributions. Thus, it is not worth to pursue SCNEVPT2 further in the present work.
ii.2 Reformulation of PCNEVPT2 in the atomic and active orbital basis
In the following we will present a reformulation of PCNEVPT2 energies in the atomic and active molecular orbital basis. Only a few working equations of PCNEVPT2 in the MO basis are presented to discuss the reformulation in sufficient detail. A complete presentation of the conventional MObased NEVPT2 can be found in Ref. Angeli, Cimiraglia, and Malrieu, 2002.
The energy term is identical to the MP2 energy that depends only on the core and virtual orbitals:
(23) 
It can be expressed in terms of orbitalenergy denominators (OED) if a canonical MO basis is employed, which keeps the the corecore and virtualvirtual block of the Fock matrix in Eq. (9) diagonal,
(24)  
(25) 
Such kind of canonical core and virtual MOs are used in the remainder of this article.
The Laplace transformation of OEDs ,
(26) 
is employed to factorize the energy term.Almlöf (1991); Häser and Almlöf (1992) The parameters of the numerical quadrature in (26) are usually obtained either by leastsquare minimizationHäser and Almlöf (1992) or by using the minimax approximation.Takatsuka, Tenno, and Hackbusch (2008); HelmichParis and Visscher (2016) Eq. (26) is the foundation of a reformulation of (or MP2) energies
(27) 
in terms of intermediates in the AO basis that are obtained from two consecutive oneindex transformation of the twoelectron integrals
(28)  
(29) 
with the core and virtual pseudodensity matrices
(30)  
(31) 
The energy and wavefunction coefficients of the contribution in the MO basis is given by
(32)  
(33) 
Solving the linear system of equations (33) is avoided if the Koopmans and metric matrix are diagonalized in the GEP
(34) 
This leads to a direct expression for the energy contribution
(35)  
(36) 
with an OED that contains Koopmans matrix eigenvalues. The number of eigenvalues in the equations above may be smaller than the number of active orbitals if becomes singular and some pairs of singular values and vectors of need to be removed for reasons of numerical stability (vide infra). To obtain an energy expressions solely in the AO and active MO basis, like for the (or MP2) term, core and virtual MO coefficients and Fock matrix eigenvalues in Eqs. (24) and (25) can be incorporated into core (Eq. (30)) and virtual (Eq. (31)) pseudodensity matrices. Concerning the activeorbital part in Eq. (35), it is convenient to summarize all those quantities that appear in GEP (34) into a single intermediate
(37) 
which we refer to as Koopmans matrix pseudoexponential. By means of , the energy can be formulated in terms of intermediates that depend only on the AOs and active MOs:
(38) 
Alternatively, a pure AObased formulation, which is more inline with standard AOMP2 formulations,
(39)  
(40) 
can be chosen by incorporating the Koopmans matrix pseudoexponential into an active pseudodensity matrix
(41) 
Likewise, the energy term (Tab. 1) can be formulated completely in the AO basis by replacing the virtual pseudodensity matrix in one of the halftransformed integrals with an active pseudodensity matrix (Tab. 1) that includes the Koopmans matrix eigenvalues describing singleelectron attachment into the active space.
The remaining energy contributions are compiled in Tab. 1 along with their intermediates. For each class, Koopmans matrix pseudoexponentials of the form
(42)  
(43) 
are computed. In Eq. (42) is the number of active orbital pairs in metric and the sign with that the Koopmans matrix eigenvalues enter the OEDs. In analogy to AOMP2, where the core and virtual pseudodensity matrices can be computed without a preceding diagonalization of the Fock matrix in the SCF procedure,Surján (2005) one could consider a formulation of the Koopmans matrix pseudoexponential where solving the GEP (21) is avoided,
(44) 
However, in case of the pseudodensity matrices in AOMP2, the inverse of the metric matrix, that is the AO overlap , is readily available in form of the core and virtual density matrix :Ayala and Scuseria (1999); Surján (2005)
(45) 
In Eq. (44) the inverted metric matrix is initially unknown and must be determined by diagonalization and subsequent removal of linearly dependent eigenvectors. Consequently, it seems that there is no gain in computational performance if one tries to circumvent solving the GEP (21).
The explicit expressions for the secondorder energy contributions in Tab. 2 show a rather surprising resemblance with the recently proposed timedependent (t) NEVPT2 of Sokolov and Chan.Sokolov and Chan (2016) Apart from an AObased formulation for the core and virtual orbitals, the two formulations merely differ in the active orbitalbased intermediates. To obtain explicit expressions for our Laplacetransformed PCNEVPT2 formulation from those of tNEVPT2, the timeordered mparticle 1time Green’s functions with (m=13),Sokolov and Chan (2016) e. g. for the space
(46) 
needs to be replaced by the corresponding timedependent Koopmans matrix pseudo exponentials
(47) 
if we consider the analytic form of the Laplace transformation in Eq. (26).
Iii Computational details
We used the CASSCF procedure as implemented in the MOLCAS packageAquilante et al. (2016) to obtain the zerothorder MOs. Subsequently, Koopmans matrices were generated by the relmrpt2 module, a locally modified version of a DMRGNEVPT2 module Freitag et al. (2017), from nRDMs obtained by means of DMRG calculations with the QCMaquis packageKeller et al. (2015); Keller and Reiher (2016); Knecht et al. (2016).
A pilotimplementation of the Laplacetransformed PCNEVPT2 in the atomic and active molecular orbital basis was integrated into the Kramersrestricted twocomponent AOMP2HelmichParis, Repisky, and Visscher (2016) implementation, which is part of a development version of the DIRAC package.DIR () The parameters of the numerical quadrature were obtained by using an implementation of the minimax approximationTakatsuka, Tenno, and Hackbusch (2008); HelmichParis and Visscher (2016), which is provided as external opensource libraryHelmichParis (). The twoelectron integrals needed for NEVPT2 were calculated with the InteRest library.Repisky (2013)
Except for the and class, the GEP (21) is plagued by numerically instabilities caused by singularities in the metric matrix . To guarantee numerical stability when solving the GEP (21), the eigenvalues of are computed first and then compared to a threshold to remove singular valuevector pairs by a canonical orthogonalization procedure.Szabo and Ostlund (1996) The linear dependency threshold is set to for all classes as in the implementation described in Ref. Angeli, Cimiraglia, and Malrieu, 2001; Aidas et al., 2014. Pointgroup (PG) symmetry cannot be exploited for our NEVPT2 implementation at the current stage. Consequently, C PG symmetry was enforced for all CASSCF and DMRG calculations, though all molecules investigated in the this work have at least two symmetry elements.
All results presented here were obtained by employing the full Fock matrix in the definition of the core part of the Dyall HamiltonianAngeli, Borini, and Cimiraglia (2004); Havenith et al. (2004) rather than only the diagonal elementsAngeli et al. (2001) to guarantee rotational invariance within all three orbital spaces. However, we omitted the offdiagonal elements in the Fock matrix when validating the correctness of the working equations in Tab. 1 and our implementation by comparing our results with the ones from the NEVPT2 implementation in DALTON.Angeli, Cimiraglia, and Malrieu (2001); Aidas et al. (2014)
Iv Results and discussion
iv.1 Discussion of the working equations
In the previous section, we outlined the derivation of the Laplacetransformed PCNEVPT2 equations formulated in the AO and active MO basis and presented explicit expressions in Tab. 2. We proposed that the energy contributions that depend on integrals with either no or a single active MO, that is the , , and term, can be formulated completely in the AO basis similar to the SR AObased MP2Häser (1993). Those three energy terms differ merely in the pseudodensity matrices; the MP2like term requires only the core and virtual pseudodensity matrices, one active pseudodensity matrix replaces one core and one virtual for the and energy term, respectively.
A Laplacetransformed AOdriven implementation of the , , and terms increases the computational costs noticeably because, first, the number of basis functions is much larger than the number of core and active orbitals and, second, the timedetermining steps have to be repeated for every quadrature point of Laplace transformed OEDs. However, all intermediates are formulated in the AO basis and become sparse for sufficiently large systems. Furthermore, it was shown for SR AOMP2 that linearly scaling implementations with an early onsetDoser et al. (2009); *Maurer2013b can be obtained as the contributions of the AO quadruple to the direct or Coulomb MP2 energy decay as with respect to the separation of electronic charge distributions .Ayala and Scuseria (1999) This rapid decay could be attributed to the vanishing overlap of core and virtual pseudodensity matrices,
(48) 
in the multipole expansion of the integrals in orders of This leads to a leading dipoledipole integral term for the direct MP2 energy. Since the overlap of core and active respectively virtual and active pseudodensity matrices vanishes as well,
(49) 
due to the mutual orthogonality of the three separate orbital spaces, the and terms reveal the same decay of the direct energy contributions as AOMP2. This beneficial longrange behavior of the contributions has been exploited already by Guo et al. in their linearly scaling PNONEVPT2 implementation at the level of electron pair prescreening.Guo et al. (2016); Pinski et al. (2015)
Though the , and energy terms show the same interelectronic decay behavior, their magnitude may differ substantially and is sensitive to the number of electrons and orbitals in the active space. As and depend linearly on their respective active pseudodensity matrix, we show contour plots of and along with and for the polyene CH (PE32) and polyene glycol biradical (PEG32) in Figs. 2 and 3, respectively.si () The calculations were performed with the SVP basis setSchäfer, Horn, and Ahlrichs (1992) and the orbitals of the C and O atoms were kept frozen. For PE32 the CAS is composed of 8 electrons that are distributed amongst 8 active valence MOs. The latter are delocalized over the entire conjugated pisystem of the polyene. As can be see from Fig. 2, the active pseudodensity matrices and of the polyene are very similar and have their maxima along the diagonal of the CC blocks. The elements of the HH and HC blocks of are much smaller in . All type basis function do not contribute to the piorbitals in the CAS for symmetry reasons and, therefore, is much sparser than and . For the PEG32 biradical we build a CAS from the 10 electrons and 6 orbitals of the two O atoms. Those active orbitals are localized at the ends of PEG32 and have natural occupation numbers of . As can be seen from Fig. 3, the core and virtual pseudodensity matrices and have their largest elements at the diagonal shell pairs and their nearest neighbors whereas the largest elements of the active pseudodensities are located at the O atoms. As expected from , has much more significant elements than . Concerning , only those basis functions are important that have sizable coefficients for the two active MOs with .
The test molecules that we investigated have different physical characteristics of their valence orbitals in the active space. While for polyene the valence MOs are delocalized over nearly the whole molecule, for the polyethylene glycol biradical the valence MO are localized at the ends of the linear molecule. For both systems, we observed that the number of significant elements in the active pseudodensity matrix is much smaller than in the core and virtual pseudodensity matrix. In combination with the discussed decay, this should pave the way for a reduced or even linearly scaling implementation where the and energies are formulated completely in the AObasis. The computational overhead of an entirely AObased formulation of the and energies can be reduced drastically when performing a CDZienau et al. (2009); Kussmann et al. (2015); Luenser, Schurkus, and Ochsenfeld () of the active pseudodensity matrices. Such a CDDbased algorithm would combine successfully the sparsity with the low rank of active pseudodensity matrices.
iv.2 Accuracy of the numerical quadrature
The presence of three different orbital spaces results in eight different contributions to the IC firstorder wavefunction contribution with eight different OEDs. Each of the 8 OEDs is approximated by its own numerical quadrature of its Laplace transform, which is determined by minimization the maximum absolute error (MAE) of the error distribution functionTakatsuka, Tenno, and Hackbusch (2008)
(50) 
This procedure is known as the minimax algorithmTakatsuka, Tenno, and Hackbusch (2008); HelmichParis and Visscher (2016) and can be performed in interval with by scaling the parameters . This facilitates pretabulation of along with their MAEs for some discrete . These parameters can be used as start values in the iterative optimization procedure.Hackbusch (2007) We choose such that
(51) 
where is the pretabulated range closest to . The accuracy of the numerical quadrature depends on a single usergiven threshold that should correlate with the resultant error in each correlation energy contribution. Alternatively, one could simply set to the same value for each of the 8 OEDs based on experience. However, by doing so, the errors in cannot be easily controlled for each class. This can be clearly seen from Fig. 1 where the 8 different errors differ from each other even by more than 4 orders of magnitude for some . Conversely, the same errors converge much more uniformly with respect to the MAE , as shown in Fig. 1.
In the following, we investigate the errors when is selected according to (51) for the ground state of the fluorine and chlorine dimer. The results for different computational setups are compiled in Tab. 4. For F and Cl (a multiple of) the equilibrium bond distances were taken from experimental data.Huber and Herzberg (2015) For F the ccpVTZDunning (1989) and augccpVTZKendall, Dunning, and Harrison (1992) basis sets were employed; for Cl we used the ccpwCVTZWoon and Jr. (1993) basis set. The Laplace accuracy threshold was set to . We note that for all calculations in Tab. 4 the part of the is zero and is not shown.
For the first two calculation of F in Tab. 4 we chose a CAS (10,6) that includes all orbitals and electrons of F. The first calculation differs from the second only in correlation of the 1s core orbitals. The 8 OED ranges of the allelectron (AE) calculation varies roughly by a factor of two only and similar , i. e. 7 and 8, are selected by our procedure. The absolute errors in are all below . By freezing the orbitals the OED ranges of those classes that involve core orbitals are decreased significantly while the and contributions remain unchanged. The latter require the least (7) in the AE calculation but the most in the frozencore (FC) calculation. For the FC calculation the OED ranges differ by almost 1 order of magnitude. Nevertheless, all absolute errors in are smaller than . Augmenting the ccpVTZ basis by one diffuse basis function for each angular momentumKendall, Dunning, and Harrison (1992) lowers the eigenvalues of the virtualvirtual block of the Fock matrix. We can observe for the third calculation in Tab. 4 that the OED ranges of those classes that involve virtual orbitals are increased approximately by a factor of 2 while those of the space are almost unaffected when adding the diffuse “aug” functions to the basis set. The absolute errors in are still reliably below .
For Cl atom in Cl we correlate the orbitals and keep the frozen. For the fourth calculation in Tab. 4 we put the all orbitals and electrons of Cl into the active space, that is CAS (10,6). For this calculation the OED ranges become much larger than those of the FC F calculations, which can be attributed to the additional core electrons and the larger basis set for Cl. The largest OED range in the fourth calculation is 114.76, while the largest in F calculation with the augccpVTZ basis set is 24.32. Though the OED ranges vary by more than 1 order of magnitude for the Cl calculation, all absolute errors in are less than and seems to be selected properly. In the fifth calculation in Tab. 4 the orbitals and electrons of the Cl atoms are included in the valence space, i. e. CAS (14,8). This lowers substantially the OED ranges of all those classes that excite from core electrons. The smallest () and the largest OED range () differ by almost 2 orders of magnitude but is selected properly as the absolute error in is always smaller . The convergence of the absolute error with respect to and the MAE is shown in Fig. 1 and has been discussed already in the beginning of this section.
The previous example calculations of F and Cl had clearly singlereference character as deduced from their natural occupation numbers ( or ) and could have been performed more easily with closedshell singlereference MP2 or coupled cluster. Therefore, we performed the CAS (14,8) calculation also on a stretched Cl with bond distance. Now there are 2 open shells with and 6 closed shells with . As we investigate the state of Cl, at least 2 CSF are required to describe the system at least qualitatively correct. Thus, he have a MR case by definition. The MR character does not seem to change the OEDs much and we obtain absolute errors in that are all smaller than . It is noteworthy that the contribution to the correlation energy almost vanishes since Cl can be considered as nearly dissociated at and it is impossible to attach 2 electrons into a CAS (5,3) or (7,4) valence space of a single Cl atom.
A further field of application of MRPT methods are excited states. With our statespecific (SS) Laplacetransformed PCNEVPT2 method we investigate the OED ranges and absolute errors in of the ground state 1 A and two lowest vertically excited singlet states 1 A () and 1 B () in formaldehyde, which are shown in Tab. 5. Those and other valence and Rydberg excited states have been investigated rigorously by Merchán and RoosMerchán and Roos (1995) with SSCASPT2 and by Angeli et al. with SSNEVPT2Angeli, Borini, and Cimiraglia (2004). The geometry and the ANO(S)Ryd(S) basis set were taken from their works. To describe also Rydberg states, a contracted s, p, and d function, denoted by Ryd(S), were located altogether in the geometric center of HCO. In contrast to Refs. Merchán and Roos, 1995 and Angeli, Borini, and Cimiraglia, 2004 we chose a CAS (4,4) that includes the orbitals to optimize the orbitals and CI coefficients of the 1 A, 1 A, and 1 B states simultaneously in the SACASSCF procedure because, currently, we cannot exploit PG symmetry for the NEVPT2 calculations with our pilot implementation. Nevertheless, the excitation energies to the 1 A (4.10 eV) and 1 B state (7.18 eV) differ only by 0.07 and +0.10 eV, respectively, though an active space was employed that differs from the ones of Ref. Angeli, Borini, and Cimiraglia, 2004.
As observed already for the halogen dimers, those classes that attach electrons into the active space have the smallest OED ranges and demand the least to reach a desired accuracy for all 3 states in Tab. 5; while active space electron ionization classes have larger OED ranges and need the most . The combined single electron ionization and excitation class has the largest OED range for each of the 3 states but does not contribute with more than 7 % to the total energy. For the Rydberg excited state 1 B the OED range of the class is only a factor of two larger than for both the ground state 1 A and first valence excited state 1 A. Furthermore, only 1 additional quadrature point is required for in 1 B to reach the requested target accuracy. Therefore, OED ranges of contributions to Rydberg excited states are larger than for valence excited states; but even if diffuse basis functions are included in the basis set, the additional number of quadrature points is still rather decent. In fact, it is the less effective screening of shell contributions in the AO basis rather than the additional number of quadrature points that makes an AObased formulation less powerful for basis sets with diffuse basis functions in a productionlevel implementation. This we will study in future work.
V Conclusions
In the present article we showed how to formulate the partially contracted nelectron valence second order perturbation theory (PCNEVPT2) energies in the atomic and active molecular orbital basis by employing the Laplace transformation of orbitalenergy denominators (OED). As the number of active orbitals is comparatively small and AO basis functions are inherently localized, our formulation is particularly suited for a linearlyscaling NEVPT2 implementation. Some of the NEVPT2 energy contributions can be formulated completely in the AO basis as singlereference secondorder Møller–Plesset perturbation theory and benefit from sparse activepseudo density matrices — particularly if the active molecular orbitals are localized only in parts of a molecule. Furthermore, we could show that finding optimal parameters of the numerical Laplace transformation is particularly challenging for multireference perturbation theories as the fit range may vary among the 8 different OEDs by many orders of magnitude. Selecting the number of quadrature points for each OED separately according to an accuracybased criterion allows us to control the errors in the NEVPT2 energies reliably. Currently, we work on efficient loworder scaling implementations of NEVPT2 and relativistic twocomponent MP2. The extension of our formalism to quasidegenerate NEVPT2Angeli et al. (2004) or the uncontracted timedependent NEVPT2Sokolov and Chan (2016) will be presented in the future.
Vi Acknowledgments
This work was supported by the the German research foundation DFG through a research fellowship (Grant No. HE 7427/11) and the Dutch science foundation NWO through a VENI fellowship (Grant No. 722.016.011). B. H.P. would like to thank Markus Reiher for granting access to the local computing facilities at ETH Zürich and the NWO for computer time at the Dutch national supercomputers. The authors also thank Hans Jørgen Aa. Jensen for sharing his experiences in multireference electronic structure methods.
References
 Martínez (2006) T. J. Martínez, Acc. Chem. Res. 39, 119 (2006).
 Subotnik et al. (2016) J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi, Annu. Rev. Phys. Chem. 67, 387 (2016).
 Lengsfield (1980) B. H. Lengsfield, III., J. Chem. Phys. 73, 382 (1980).
 Werner and Meyer (1980) H. Werner and W. Meyer, J. Chem. Phys. 73, 2342 (1980).
 Werner and Meyer (1981) H. Werner and W. Meyer, J. Chem. Phys. 74, 5794 (1981).
 Werner and Knowles (1985) H. Werner and P. J. Knowles, J. Chem. Phys. 82, 5053 (1985).
 Roos, Taylor, and Siegbahn (1980) B. O. Roos, P. R. Taylor, and P. E. Siegbahn, Chem. Phys. 48, 157 (1980).
 Knowles and Handy (1984) P. Knowles and N. Handy, Chem. Phys. Lett. 111, 315 (1984).
 Olsen, Jørgensen, and Simons (1990) J. Olsen, P. Jørgensen, and J. Simons, Chem. Phys. Lett. 169, 463 (1990).
 Olsen et al. (1988) J. Olsen, B. O. Roos, P. Jørgensen, and H. J. A. Jensen, J. Chem. Phys. 89, 2185 (1988).
 Fleig, Olsen, and Marian (2001) T. Fleig, J. Olsen, and C. M. Marian, J. Chem. Phys. 114, 4775 (2001).
 Booth, Thom, and Alavi (2009) G. H. Booth, A. J. W. Thom, and A. Alavi, J. Chem. Phys. 131, 054106 (2009).
 Thomas et al. (2015) R. E. Thomas, Q. Sun, A. Alavi, and G. H. Booth, J. Chem. Theory Comput. 11, 5316 (2015).
 White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
 White (1993) S. R. White, Phys. Rev. B 48, 10345 (1993).
 Zgid and Nooijen (2008a) D. Zgid and M. Nooijen, J. Chem. Phys. 128, 144116 (2008a).
 Zgid and Nooijen (2008b) D. Zgid and M. Nooijen, J. Chem. Phys. 128, 144115 (2008b).
 Ghosh et al. (2008) D. Ghosh, J. Hachmann, T. Yanai, and G. K.L. Chan, J. Chem. Phys. 128, 144117 (2008).
 Yanai et al. (2009) T. Yanai, Y. Kurashige, D. Ghosh, and G. K.L. Chan, Int. J. Quantum Chem. 109, 2178 (2009).
 Kurashige and Yanai (2011) Y. Kurashige and T. Yanai, J. Chem. Phys. 135, 094104 (2011).
 Wouters et al. (2014) S. Wouters, W. Poelmans, P. W. Ayers, and D. van Neck, Comput. Phys. Commun. 185, 1501 (2014).
 Wouters et al. (2015) S. Wouters, W. Poelmans, S. De Baerdemacker, P. W. Ayers, and D. Van Neck, Comput. Phys. Commun. 191, 235 (2015).
 Ma et al. (2016) Y. Ma, S. Knecht, S. Keller, and M. Reiher, ArXiv eprints (2016), preprint: arXiv:1611.05972.
 Li Manni, Smart, and Alavi (2016) G. Li Manni, S. D. Smart, and A. Alavi, J. Chem. Theory Comput. 12, 1245 (2016).
 Knecht et al. (2016) S. Knecht, E. D. Hedegård, S. Keller, A. Kovyrshin, Y. Ma, A. Muolo, C. J. Stein, and M. Reiher, Chimia 70, 244 (2016).
 Sun, Yang, and Chan (2017) Q. Sun, J. Yang, and G. K.L. Chan, Chem. Phys. Lett. (2017), in press.
 Hohenstein et al. (2015) E. G. Hohenstein, N. Luehr, I. S. Ufimtsev, and T. J. Martínez, J. Chem. Phys. 142, 224103 (2015).
 Chaban, Schmidt, and Gordon (1997) G. Chaban, M. W. Schmidt, and M. S. Gordon, Theor. Chem. Acc. 97, 88 (1997).
 Andersson et al. (1990) K. Andersson, P. Å. Malmqvist, B. O. Roos, A. J. Sadlej, and K. Wolinski, J. Phys. Chem. 94, 5483 (1990).
 Andersson, Malmqvist, and Roos (1992) K. Andersson, P.Å. Malmqvist, and B. O. Roos, J. Chem. Phys. 96, 1218 (1992).
 Celani and Werner (2003) P. Celani and H.J. Werner, J. Chem. Phys. 119, 5044 (2003).
 Shiozaki et al. (2011) T. Shiozaki, W. Győrffy, P. Celani, and H.J. Werner, J. Chem. Phys. 135, 081106 (2011).
 Győrffy et al. (2013) W. Győrffy, T. Shiozaki, G. Knizia, and H.J. Werner, J. Chem. Phys. 138, 104104 (2013).
 Shiozaki and Werner (2010) T. Shiozaki and H.J. Werner, J. Chem. Phys. 133, 141103 (2010).
 Shiozaki and Werner (2013) T. Shiozaki and H.J. Werner, Mol. Phys. 111, 607 (2013).
 Roos and Malmqvist (2004) B. O. Roos and P.Å. Malmqvist, Phys. Chem. Chem. Phys. 6, 2919 (2004).
 Abe, Nakajima, and Hirao (2006) M. Abe, T. Nakajima, and K. Hirao, J. Chem. Phys. 125, 234110 (2006).
 Finley et al. (1998) J. Finley, P.Å. Malmqvist, B. O. Roos, and L. SerranoAndrés, Chem. Phys. Lett. 288, 299 (1998).
 Angeli et al. (2004) C. Angeli, S. Borini, M. Cestari, and R. Cimiraglia, J. Chem. Phys. 121, 4043 (2004).
 Menezes, Kats, and Werner (2016) F. Menezes, D. Kats, and H.J. Werner, J. Chem. Phys. 145, 124115 (2016).
 Roos and Andersson (1995) B. O. Roos and K. Andersson, Chem. Phys. Lett. 245, 215 (1995).
 Forsberg and Malmqvist (1997) N. Forsberg and P.Å. Malmqvist, Chem. Phys. Lett. 274, 196 (1997).
 Angeli, Cimiraglia, and Malrieu (2000) C. Angeli, R. Cimiraglia, and J.P. Malrieu, Chem. Phys. Lett. 317, 472 (2000).
 Angeli et al. (2001) C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger, and J.P. Malrieu, J. Chem. Phys. 114, 10252 (2001).
 Angeli, Cimiraglia, and Malrieu (2001) C. Angeli, R. Cimiraglia, and J.P. Malrieu, Chem. Phys. Lett. 350, 297 (2001).
 Angeli, Cimiraglia, and Malrieu (2002) C. Angeli, R. Cimiraglia, and J.P. Malrieu, J. Chem. Phys. 117, 9138 (2002).
 Dyall (1995) K. G. Dyall, J. Chem. Phys. 102, 4909 (1995).
 Sharma and Chan (2014) S. Sharma and G. K.L. Chan, J. Chem. Phys. 141, 111101 (2014).
 Sokolov and Chan (2016) A. Y. Sokolov and G. K.L. Chan, J. Chem. Phys. 144, 064102 (2016).
 Giner et al. (2017) E. Giner, C. Angeli, Y. Garniron, A. Scemama, and J.P. Malrieu, ArXiv eprints (2017), arXiv:1702.03133 [physics.chemph] .
 Aquilante et al. (2008) F. Aquilante, P.Å. Malmqvist, T. B. Pedersen, A. Ghosh, and B. O. Roos, J. Chem. Theory Comput. 4, 694 (2008).
 Freitag et al. (2017) L. Freitag, S. Knecht, C. Angeli, and M. Reiher, J. Chem. Theory Comput. 13, 451 (2017).
 Aquilante et al. (2009) F. Aquilante, T. K. Todorova, L. Gagliardi, T. B. Pedersen, and B. O. Roos, J. Chem. Phys. 131, 034113 (2009).
 Guo et al. (2016) Y. Guo, K. Sivalingam, E. F. Valeev, and F. Neese, J. Chem. Phys. 144, 094111 (2016).
 Almlöf (1991) J. Almlöf, Chem. Phys. Lett. 181, 319 (1991).
 Häser and Almlöf (1992) M. Häser and J. Almlöf, J. Chem. Phys. 96, 489 (1992).
 Häser (1993) M. Häser, Theor. Chem. Acc. 87, 147 (1993).
 Ayala and Scuseria (1999) P. Y. Ayala and G. E. Scuseria, J. Chem. Phys. 110, 3660 (1999).
 Lambrecht, Doser, and Ochsenfeld (2005) D. S. Lambrecht, B. Doser, and C. Ochsenfeld, J. Chem. Phys. 123, 184102 (2005).
 Doser et al. (2009) B. Doser, D. Lambrecht, J. Kussmann, and C. Ochsenfeld, J. Chem. Phys. 130, 064107 (2009).
 Maurer et al. (2013) S. A. Maurer, D. S. Lambrecht, J. Kussmann, and C. Ochsenfeld, J. Chem. Phys. 138, 014101 (2013).
 Zienau et al. (2009) J. Zienau, L. Clin, B. Doser, and C. Ochsenfeld, J. Chem. Phys. 130, 204112 (2009).
 Kussmann et al. (2015) J. Kussmann, A. Luenser, M. Beer, and C. Ochsenfeld, J. Chem. Phys. 142, 094101 (2015).
 (64) A. Luenser, H. F. Schurkus, and C. Ochsenfeld, J. Chem. Theory Comput. http://dx.doi.org/10.1021/acs.jctc.6b01235 .
 HelmichParis, Repisky, and Visscher (2016) B. HelmichParis, M. Repisky, and L. Visscher, J. Chem. Phys. 145, 014107 (2016).
 Celani and Werner (2000) P. Celani and H.J. Werner, J. Chem. Phys. 112, 5546 (2000).
 Pulay, Sæbø, and Meyer (1984) P. Pulay, S. Sæbø, and W. Meyer, J. Chem. Phys. 81, 1901 (1984).
 Helgaker, Jørgensen, and Olsen (2000a) T. Helgaker, P. Jørgensen, and J. Olsen, “Molecular electronicstructure theory,” (Wiley, New York, 2000) Chap. 13, p. 692.
 Helgaker, Jørgensen, and Olsen (2000b) T. Helgaker, P. Jørgensen, and J. Olsen, “Molecular electronicstructure theory,” (Wiley, New York, 2000) Chap. 10, p. 483.
 van Dam, van Lenthe, and Pulay (1998) H. J. J. van Dam, J. H. van Lenthe, and P. Pulay, Mol. Phys. 93, 431 (1998).
 van Dam, van Lenthe, and Ruttink (1999) H. J. J. van Dam, J. H. van Lenthe, and P. J. A. Ruttink, Int. J. Quantum Chem. 72, 549 (1999).
 Morrell, Parr, and Levy (1975) M. M. Morrell, R. G. Parr, and M. Levy, J. Chem. Phys. 62, 549 (1975).
 Smith and Day (1975) D. W. Smith and O. W. Day, J. Chem. Phys. 62, 113 (1975).
 Takatsuka, Tenno, and Hackbusch (2008) A. Takatsuka, S. Tenno, and W. Hackbusch, J. Chem. Phys. 129, 044112 (2008).
 HelmichParis and Visscher (2016) B. HelmichParis and L. Visscher, J. Comput. Phys. 321, 927 (2016).
 Surján (2005) P. R. Surján, Chem. Phys. Lett. 406, 318 (2005).
 Aquilante et al. (2016) F. Aquilante, J. Autschbach, R. K. Carlson, L. F. Chibotaru, M. G. Delcey, L. De Vico, I. Fdez. Galván, N. Ferré, L. M. Frutos, L. Gagliardi, M. Garavelli, A. Giussani, C. E. Hoyer, G. Li Manni, H. Lischka, D. Ma, P. Å. Malmqvist, T. Müller, A. Nenov, M. Olivucci, T. B. Pedersen, D. Peng, F. Plasser, B. Pritchard, M. Reiher, I. Rivalta, I. Schapiro, J. SegarraMartí, M. Stenrup, D. G. Truhlar, L. Ungur, A. Valentini, S. Vancoillie, V. Veryazov, V. P. Vysotskiy, O. Weingart, F. Zapata, and R. Lindh, J. Comput. Chem. 37, 506 (2016).
 Keller et al. (2015) S. Keller, M. Dolfi, M. Troyer, and M. Reiher, J. Chem. Phys. 143, 244118 (2015).
 Keller and Reiher (2016) S. Keller and M. Reiher, J. Chem. Phys. 144, 134101 (2016).
 (80) DIRAC, a relativistic ab initio electronic structure program, Release DIRAC16 (2016), written by H. J. Aa. Jensen, R. Bast, T. Saue, and L. Visscher, with contributions from V. Bakken, K. G. Dyall, S. Dubillard, U. Ekström, E. Eliav, T. Enevoldsen, E. Faßhauer, T. Fleig, O. Fossgaard, A. S. P. Gomes, T. Helgaker, J. Henriksson, M. Iliaš, Ch. R. Jacob, S. Knecht, S. Komorovský, O. Kullie, J. K. Lærdahl, C. V. Larsen, Y. S. Lee, H. S. Nataraj, M. K. Nayak, P. Norman, G. Olejniczak, J. Olsen, Y. C. Park, J. K. Pedersen, M. Pernpointner, R. di Remigio, K. Ruud, P. Sałek, B. Schimmelpfennig, J. Sikkema, A. J. Thorvaldsen, J. Thyssen, J. van Stralen, S. Villaume, O. Visser, T. Winther, and S. Yamamoto (see http://www.diracprogram.org).
 (81) B. HelmichParis, laplaceminimax library release v1.5 available from https://github.com/bhelmichparis/laplaceminimax.git.
 Repisky (2013) M. Repisky, (2013), InteRest 2.0, An integral program for relativistic quantum chemistry.
 Szabo and Ostlund (1996) A. Szabo and N. S. Ostlund, “Modern quantum chemistry: Introduction to advanced electronic structure theory,” (Dover Publications, New York, USA, 1996) Chap. 3, p. 144.
 Aidas et al. (2014) K. Aidas, C. Angeli, K. L. Bak, V. Bakken, R. Bast, L. Boman, O. Christiansen, R. Cimiraglia, S. Coriani, P. Dahle, E. K. Dalskov, U. Ekström, T. Enevoldsen, J. J. Eriksen, P. Ettenhuber, B. Fernández, L. Ferrighi, H. Fliegl, L. Frediani, K. Hald, A. Halkier, C. Hättig, H. Heiberg, T. Helgaker, A. C. Hennum, H. Hettema, E. Hjertenæs, S. Høst, I.M. Høyvik, M. F. Iozzi, B. Jansík, H. J. A. Jensen, D. Jonsson, P. Jørgensen, J. Kauczor, S. Kirpekar, T. Kjærgaard, W. Klopper, S. Knecht, R. Kobayashi, H. Koch, J. Kongsted, A. Krapp, K. Kristensen, A. Ligabue, O. B. Lutnæs, J. I. Melo, K. V. Mikkelsen, R. H. Myhre, C. Neiss, C. B. Nielsen, P. Norman, J. Olsen, J. M. H. Olsen, A. Osted, M. J. Packer, F. Pawlowski, T. B. Pedersen, P. F. Provasi, S. Reine, Z. Rinkevicius, T. A. Ruden, K. Ruud, V. V. Rybkin, P. Sałek, C. C. M. Samson, A. S. de Merás, T. Saue, S. P. A. Sauer, B. Schimmelpfennig, K. Sneskov, A. H. Steindal, K. O. SylvesterHvid, P. R. Taylor, A. M. Teale, E. I. Tellgren, D. P. Tew, A. J. Thorvaldsen, L. Thøgersen, O. Vahtras, M. A. Watson, D. J. D. Wilson, M. Ziolkowski, and H. Ågren, WIREs Comput. Mol. Sci. 4, 269 (2014).
 Angeli, Borini, and Cimiraglia (2004) C. Angeli, S. Borini, and R. Cimiraglia, Theor. Chem. Acc. 111, 352 (2004).
 Havenith et al. (2004) R. W. A. Havenith, P. R. Taylor, C. Angeli, R. Cimiraglia, and K. Ruud, J. Chem. Phys. 120, 4619 (2004).
 Pinski et al. (2015) P. Pinski, C. Riplinger, E. F. Valeev, and F. Neese, J. Chem. Phys. 143, 034108 (2015).
 (88) Molecular structures of and are provided in the Supporting Information.
 Schäfer, Horn, and Ahlrichs (1992) A. Schäfer, H. Horn, and R. Ahlrichs, J. Chem. Phys. 97, 2571 (1992).

Hackbusch (2007)
W. Hackbusch, “
Pretabulated Laplace exponents from http://www.mis.mpg.de/scicomp/EXP_SUM/1_x/,” (2007).  Huber and Herzberg (2015) K. Huber and G. Herzberg, in NIST Chemistry WebBook, NIST Standard Reference Database Number 69, edited by P. Linstrom and W. Mallard (National Institute of Standards and Technology, Gaithersburg MD, 20899, 2015) data prepared by J.W. Gallagher and R.D. Johnson, III.
 Dunning (1989) T. H. Dunning, J. Chem. Phys. 90, 1007 (1989).
 Kendall, Dunning, and Harrison (1992) R. A. Kendall, T. H. Dunning, and R. J. Harrison, J. Chem. Phys. 96, 6796 (1992).
 Woon and Jr. (1993) D. E. Woon and T. H. D. Jr., J. Chem. Phys. 98, 1358 (1993).
 Merchán and Roos (1995) M. Merchán and B. O. Roos, Theor. Chem. Acc. 92, 227 (1995).














