Introducing explicitly correlated coupled cluster approaches into the world of astrophysics

Pierre Valiron’s contribution in the field of explicitly correlated coupled cluster approaches is briefly overviewed. Beside parallelization of the DIRCCR12-OS code, Pierre was involved in its development for open shell systems. His contribution towards a systematic development of special basis sets suited to explicitly correlated calculations within the R12 framework is appreciated in particular.


Introduction
In 1995 I met Pierre Valiron for the first time on a Summer School on Electron Correlation in Atoms and Molecules in a lovely Hungarian city of Debrecen, where we both lectured on our own topics.I about introducing the explicitly correlated R12 ansatz into the coupled cluster (CC) wave function expansion and Pierre about astrophysical issues.It turned out that our interests were complementary -he was looking for highly accurate methods and I was looking for a field of application for such approaches known for their not too favorable computational scaling with respect to the size of the investigated system.Little later in 1998, our real collaboration started and it was by far not restricted to astrophysical applications.In this short overview, I'll address the main issues related to our common developmental work in the direction towards a larger applicability of the CC-R12 approach.

The R1theory
In molecular electronic structure calculations it is known for a long time that the convergence to the exact solution of the Schrödinger equation with extension of the one electron basis is extremely slow within the standard electron correlation treatments which are based on wave function expansion in terms of antisymmetrized products of one-electron functions (Slater determinants) [1].The main reason is that the Coulomb hole cannot be correctly described within such approaches.For small interelectronic distance r 12 the wave function should behave as: which implies lim ( ∂Ψ ∂r 12 ) Relation ( 2) is known as Kato's electron-electron cusp condition [2].A configuration interaction (CI) expansion in orbital products does not reproduce Eq. (2).To satisfy this condition, it is, however, a e-mail: jozef.noga@savba.skEPJ Web of Conferences sufficient to introduce terms linear in r 12 to an orbital product expansion, which was the basic idea behind the R12 theory [3].The exact many electron wave function (Ψ ) in this theory is described as: Here, Φ is a reference single determinant that usually follows from the Hartree-Fock solution.Second term corresponds to the usual (conventional) configuration interaction type expansion.The operator R is related to the correlation factor [ f (r 12 )] that was originally represented by the linear r 12 .The content of the first term of (3) that overlaps with the second one has to be outprojected.In that case the many particle integrals that would arise in the working equations can be factorized into products of integrals over two particles using resolution of identity either in the computational or an auxiliary basis under the condition that the pertinent basis is saturated at the level of 3L occ , where L occ is the highest angular momentum function involved in the orbital product representation of Φ.At the same time, the convergence pattern towards the basis set limit changes from a very slowly decreasing error ∝(L+1) −3  for the conventional (CI) type expansion to ∝(L+1) −7 for an expansion as Eq. ( 3).Here, L is the highest angular momentum involved in the one particle basis set used to construct the configuration space [4].While a decade ago the r 12 as the correlation factor was almost exclusively used in R12 based theories [5,6], in the last decade, the use of Slater type geminal [7] with the length-scale parameter γ, represents the main stream [8].
are normal ordered n-body replacement operators and are matrix elements represented by the two-particle integrals over f (r 12 ).κ, λ, µ, ν formally run over the (imaginary) complete one particle basis, whereas "i, j, k, l" are reserved for occupied orbitals.Acting of FN on Φ results in: where α, β denote orbitals from a complete virtual space.The second term can be neglected if the basis set is close to the Hartree-Fock limit one, which was one of the R12 theory constraints until quite recently.The recent development is described elsewhere [9].

R12 within the coupled cluster ansatz
Evidently, the replacement operators remaining on the r.h.s. of Eq. ( 9) are excitation operators that commute with the conventional global excitation operator for n-electron system (a denote virtual orbitals within the given one-particle basis):  Based on this fact, we suggested a theory that combines the CC ansatz [10] with the R12 approach, giving rise to the CC-R12 theory [11].At the same time, we suggested diagrammatic technique to derive the pertinent working equations.
The exact wave function in the CC-R12 is expressed via exponential cluster expansion: where the operator R is defined as Evidently, R is derived from F , but via the "c"-parameters the wave function possesses more flexibility.As well, this ansatz leads to invariance with respect to rotations within the occupied orbital subspace [12].I shall skip other details here, and refer the reader to the aforementioned original literature.

Towards the open shell systems.
As I was visiting LAOG for the first time in 1998, relatively efficient sequential production code with direct calculation of two-electron integrals (DIRCCR12) [13] was available for closed shell systems.
With Pierre, we started working in two directions: i) parallelization of the code and ii) introducing the method for open shell systems.Soon after, we reported implementation of the R12 CC theory with single and double excitations (CCSD-R12) for open shell systems within the UHF formalism [14].With the aid of Pierre's expertise the code was then partially parallelized and extended with inclusion of non-iterative triple excitations for open shell systems.As well, the new version allowed using of restricted open shell reference.First applications were devoted to a hard test -calculation of atomization energies [15].Dramatic improvement of the convergence towards the basis set limit is demonstrated on Fig. 1, showing the normal distributions of errors for 7 tested molecules both for conventional CCSD(T) and CCSD(T)-R12 results.As it is seen, with a basis set saturated at spd f orbitals for a non-Hydrogen atom and spd for the Hydrogen, one essentially obtains the basis set limit values with very little uncertainty.Conventional calculations with 19s14p8d6 f 4g3h/9s6p4d3 f basis still provide an average error about 5 kJ mol −1 and a relatively wide distribution of errors.

Improvement of the algorithm for the triple excitations
The most time consuming part of the CCSD(T)-R12 calculations is represented by the contribution from triply excited configurations.This procedure involves ∝N 7 algebraic operations, where N is the number of basis functions.During the parallelization of our code which has been finished in 2003 [16], we have discovered a tricky way to improve upon the established algorithm(s) for this contribution by a speedup factor of 1.6-1.8 in real calculations.This was accomplished by converting some of the ∝N 7 steps to ∝N 6 using essentially very simple algebra, though with not obvious consequences.For a matrix C created from M arbitrary matrices A m and B m as we can write If the second sum in 14 could be obtained cheaper than with (M − 1) A m B m matrix multiplications, one would gain.This was the case in our algorithm.In the algorithm for triples, C represents a sixdimensional matrix which is traditionally obtained by six matrix multiplications of four-dimensional matrices A m and B m over a single index (i.e. ∝N 7 steps).We realized that if we use Eq. ( 14), the second term always contributes to a six-dimensional C matrix with some indices repeated, i. e. effectively to a five-dimensional matrix.Hence, also the A m B n matrix multiplications scale as ∝N 6 [17].
4 Towards the universal R12-suited basis sets.
In order to be able to factorize the many particle integrals that arise in the R12 theory into products of two-particle integrals, the basis set has to be saturated at the level of 3L occ , which together with other constraints is used to be denoted as the "standard approximation" [4].In the calculations related to Fig. 1 we used a 19s14p8d6 f 4g3h basis set for non-Hydrogen atoms, whose 19s14p core corresponded to the 18s13p Hartree-Fock limit basis [18] augmented with one diffuse s and p functions by a mere logarithmic extrapolation.Polarization functions were taken from the aug-cc-pV6Z basis [19], while two tight d functions and one tight f function were added.The 9s6p4d3 f basis for the Hydrogen atom has been taken from the aug-cc-pV5Z set that comprises 9s5p4d3 f 2g uncontracted functions.The 9s5p4d3 f subset was further augmented with one tight p function (α=11.9).These sets were denoted as "parent".Although such basis sets worked well, there were good reasons to look for improvements.One could certainly assume that optimization at the conventional level necessarily interferes with some partial compensation for the poor description of the electronic cusp.Hence, the range of exponents for the Gaussian functions might not be optimal for R12 calculations.It was likely that specifically R12optimized bases could be either smaller for comparable quality, or would provide more accurate results than the aforementioned R12 variants of the conventional sets.Following the idea of Pierre, we looked for basis sets that would be capable of describing different target systems with high (and similar) accuracy.Appropriate target systems were positively charged, neutral and negatively charged atoms.In the case of Hydrogen the H 2 molecule has been taken instead of the positively charged atom.We started from the aforementioned sets, while d, f , g and h sets were iteratively optimized as even tempered sets for the non-Hydrogen atoms [20] and similarly p, d, f sets for the Hydrogen atom [21].Hence, in each iteration we performed a two parameter optimization (α l , β l ) for the set of n l functions in the given shell l.The i-th function was defined as: This optimization gave rise to an unexpected positive side effect.Using extended basis sets in the R12 calculations has been often accompanied by numerical instabilities that were not related to the linear dependences in the one-particle basis, but to a linear dependences in the geminal (R12-pair) basis created from the pertinent one-particle basis.During the optimization procedure, we revealed ranges of parameters that should be "forbidden", because with those parameters the numerical difficulties were highly probable [22].A spectacular evidence is given on Fig. 2 which contains the energy map for C − as a function of the two parameters for the d set.The black dot corresponds to approximate parameters related to the parent basis set.Evidently, these parameters are very close to the numerical instability region.Consequently, in some molecular systems where the carbon atom is associated with a partial negative charge (like e. g.CO) the calculations were numerically instable.New sets are designed to avoid these areas of numerical instability.

Towards astrophysical applications.
The joint effort with Pierre was aimed to develop a quantum chemical tool to obtain highly accurate potential energy surfaces needed in many astrophysical applications.Indeed, first applications of our new code and basis sets in calculations of vibrational states for H 3 O + [23] and NH 3 [24] proved this ambition.The code was used for calculation of the 9-dimensional potential energy surface of H 2 O-H 2 interaction [25,26], as well as for the H 2 -CO [27].Pierre had great plans to use it in many further applications which have been stopped by the merciless fortune.

Fig. 1 .
Fig. 1.Normal distribution ρ(∆ AE ) of the deviations ∆ AE (in kJ/mol) of the CCSD(T) atomization energies from the experimental values.Results correspond to a set of 7 molecules CH 2 , NH 3 , H 2 O, HF, N 2 , CO, and F 2 calculated using subsets of the 19s14p8d6 f 4g3h/9s6p4d3 f (non-H/H) basis.Based on the data from [15].

Fig. 2 .
Fig.2.CCSD(T)-R12 energy map E=E(α d , β d ) for C − .Contours correspond to differences (in µE h ) with respect to the minimum of the energy surface.Larger differences are denoted by the color scale on the r.h.s.A 19s14p8d6 f basis was used, where 8d and 6 f functions form even-tempered series.A "stormy-like" numerical instability region is seen in the upper right triangle.The black dot corresponds to approximate parameters for the conventional parent basis.