New Transcorrelated Method Improving the Feasibility of Explicitly Correlated Calculations

Abstract

We recently developed an explicitly correlated method using the transcorrelated Hamiltonian, which is preliminarily parameterized in such a way that the Coulomb repulsion is compensated at short inter-electronic distances. The extra part of the effective Hamiltonian features short-ranged, size-consistent, and state-universal. The localized and frozen nature of the correlation factor makes the enormous three-body interaction less important and enables us to bypass the complex nonlinear optimization. We review the basic strategy of the method mainly focusing on the applications to single-reference many electron theories using modified Møller-Plesset partitioning and biorthogonal orbitals. Benchmark calculations are performed for 10-electron systems with a series of basis sets.

1. INTRODUCTION

In modern ab initio quantum chemistry, the convergence of calculations with the size of one-electronic basis functions is crucial for predicting reliable energetics and molecular properties. It has been, however, recognized that an enormous basis with high angular momentum functions is required to achieve chemical accuracy of a few kcal/mol in total energies. The slow convergence of the dynamic correlation effects is a direct consequence of the inability of describing the correlation cusp [1] with one-electronic basis. Various alternatives incorporating the inter-electronic coordinate explicitly have been suggested to ameliorate this feature. With the objectives of many electronic systems, the Gaussian-type geminals (GTGs) [2,3] have been used in many places. Especially, Szalewicz and coworkers have used GTG successfully for accurate calculations of correlation energies in their series of papers [4,5]. Although the function form is not suitable from the viewpoint of the cusp condition, the methods with GTGs have demonstrated the efficiency [6,7] to give correlation energies to μHartree accuracy, as shown in the recent result of Cencek and Rychlewski [8]. This situation is comprehensible by evaluating the motion of pair-electrons in three-dimension, i.e. the gain from the cusp around a fixed electron scales as $4\pi r_{12}^2$ while the Coulomb repulsion factorize no more than $r_{12}^{-1}$.

The R12 method [9,10,11,12] developed by Klopper, Kutzelnigg and coworkers is another original scheme, which incorporates the linear $r_{12}$ behavior to describe the dynamic correlation effects. The method introduces systematic approximations, which bypass the explicit treatment of the so-called “difficult integrals”. The approximations become more reliable as the one-electronic expansion set, which is formally assumed to be complete, gets sufficiently large. The R12 ansatz has provided highly-reliable correlated methods. Especially, the R12 coupled-cluster (CC) method developed by Noga et al. [13,14,15] achieves the chemical accuracy both in the quality of the Hilbert space to describe the electronic motions and the sophistication of the correlated method. One disadvantage of the method is that the theoretical construction is intricate because of the wave operator form and the approximation to the various difficult integrals. The long-range nature of the linear $r_{12}$ behavior is another element that makes a scalable treatment complicated. Persson and Taylor developed another method related with the R12 theory [16]. They fit the linear $r_{12}$ function using GTG to avoid the nonlinear optimization for GTG. The approach aims at calculations of modest accuracy in comparison with the previous GTG theories, i.e. sub-mHartree accuracy for valence correlations within the framework of the second order Møller-Plesset perturbation (MP2) theory.

All of the above-mentioned methods use functions with inter-electronic coordinate and orthogonal projectors. In other words, the explicitly correlated wave operators are not commutable with the Coulomb interactions in the complete basis limit. This property is the main source complicating the explicitly correlated methods. The transcorrelated method [17,18] developed by Boys and Handy is on the basis of symmetric correlation factor commutable with the Coulomb interaction. As a result, the effective Hamiltonian, the so-called transcorrelated Hamiltonian, includes at most three-body effective interactions. In recent publications, we developed a method which uses the transcorrelated Hamiltonian especially for pair-wise short-range collisions [19,20,21]. The method does not include any nonlinear optimization of the correlation factor, which is present in the original transcorrelated method. Instead, we deal with the long-range correlation in terms of the usual configuration interaction (CI)-type expansion with one-electronic functions.

In this paper, we briefly review the basic ideas of our transcorrelated method in the next section. The manipulations of molecular integrals are explained in Sec. 3. We formulate the application to single-reference methods including the linearized CC (LCC) method in Sec. 4. Some benchmark calculations of 10-electronic systems are given in Sec. 5. Conclusions are depicted in Sec. 6.

2. TRANSCORRELATED HAMILTONIAN

In our transcorrelated method, the exact wave function is expressed as a product of a symmetric correlation factor, $\text{exp}(F)$, and an anti-symmetrized CI-type wave function, $\psi$ [19],

$$\Psi =\text{exp}(F)\psi ,$$

(2.1)

where $F$ is assumed to be a sum of two-electron functions (geminals),

$$F={\displaystyle \sum _{i>j}f(r_i,r_j)},$$

(2.2)

We aim at improving the convergence of CI treating the cusp behavior in the vicinity, $r_{ij}=0$, with the factor. Since the vicinity is dominated by the Coulomb repulsion, the simplest choice for the purpose is to employ a geminal, which is also spherically symmetric,

$$f(r_i,r_j)=f(|r_{ij}|)\equiv f_{ij}.$$

(2.3)

The similarity transformed effective Hamiltonian terminating at the double commutator [17,18],

$$\begin{gathered} \tilde{H}=\text{exp}(-F)H\text{exp}(F) \\ =H+[H,F]+\frac{1}{2}[[H,F],F], \end{gathered}$$

(2.4)

is a sum of the original Hamiltonian and the effective two- and three-electronic interactions,

$$\tilde{H}=H+K+L,$$

(2.5)

$$K=\frac{1}{2}{\displaystyle \sum _{pqrs}\langle pq|K_{12}|rs\rangle a_p^{+}{a}_q^{+}{a}_s{a}_r},$$

(2.6)

$$L=\frac{1}{6}{\displaystyle \sum _{pqrstu}\langle pqr|L|stu\rangle _{123}{a}_p^{+}{a}_q^{+}{a}_r^{+}{a}_u{a}_t{a}_s},$$

(2.7)

The two-electronic part consists of the terms linear and quadratic to the factor,

$$K_{12}=K_{12}^{(L)}+K_{12}^{(Q)},$$

(2.8)

$$K_{12}^{(L)}=-{\nabla }_1^2{f}_{12}-2{\nabla }_1{f}_{12}\cdot {\nabla }_1,$$

(2.9)

$$K_{12}^{(Q)}=-{\nabla }_1{f}_{12}\cdot {\nabla }_1{f}_{12},$$

(2.10)

and the three-electron operator is given by

$$L_{123}=-3{\nabla }_1{f}_{12}\cdot {\nabla }_1{f}_{13}.$$

(2.11)

The operators, $K_{12}^{(L)}$ and $K_{12}^{(Q)}$, are anti-symmetric and symmetric, respectively. The latter can be combined with the Coulomb one in a practical computer code. It is possible to replace the operators by the averaged ones,

$$K=\frac{1}{2}{\displaystyle \sum _{pqrs}\langle pq|{\overline{K}}_{12}|rs\rangle a_p^{+}{a}_q^{+}{a}_s{a}_r},$$

(2.12)

$$L=\frac{1}{6}{\displaystyle \sum _{pqrstu}\langle pqr|\overline{L}|stu\rangle _{123}{a}_p^{+}{a}_q^{+}{a}_r^{+}{a}_u{a}_t{a}_s},$$

(2.13)

$${\overline{K}}_{12}=\frac{1}{2}(K_{12}^{(L)}+K_{21}^{(L)})+K_{12}^{(Q)},$$

(2.14)

$${\overline{L}}_{123}=\frac{1}{3}(L_{123}+L_{231}+L_{312}).$$

(2.15)

This simplifies the algebraic manipulations, which appear in correlated wave function methods.

The geminal is parameterized under the following conditions. Firstly, the Coulomb singularity is compensated by the effective potential in $K_{12}$. It is known that the condition leads to the linear $r_{12}$ behavior of the geminal, which enables us to express the Coulomb cusp explicitly with a finite number of one-electronic functions. Secondly, we localize the geminal. The linear $r_{12}$ behavior is generally adequate only in the vicinity, $r_{12}=0$, where the Coulomb potential is dominating the inter-electronic shape. We maintain the usual CI picture in the long-range behavior of the correlated wave function. For practical applications, the geminal is expressed as a sum of Gaussian-type functions,

$$f_{12}=-{\displaystyle \sum _{G=1}^{N_G}{c}_G\text{exp}(-{\zeta }_G{r}_{12}^2)},$$

(2.16)

and the parameters are determined using a least square fit for the approximate equality,

$$r_{12}^{-1}w(r_{12})-{\nabla }_1^2{f}_{12}-{\nabla }_1{f}_{12}\cdot {\nabla }_1{f}_{12}\cong 0,$$

(2.17)

where $w(r_{12})$ is a short-range weight function [19]. We have used an additional Gaussian for the weight function. The above condition involves no more than the derivatives of $f_{12}$. An additional requirement is therefore the geminal converges to zero at the infinite distance. This is automatically fulfilled with the present use of Gaussian-type functions. The choice of the localized geminal is crucial because the number of integrals increases only linearly with the size of the molecule. Moreover, the three-body interaction, $L_{123}$, becomes less important because of the decreasing probability of finding three electrons in the geminal radius. This is consistent with the localized nature of the dynamic correlation effects.

We use the tight geminal consisting of 10 Gaussian-type functions throughout of this paper. The exponents are determined as an even-tempered sequence of the function with the exponents between 904000.0 and 0.12, and that of the weight Gaussian is chosen to be 20.0. Fig.1 shows the shape of the geminal. The geminal coincides with the linear $r_{12}$ behavior around the origin. It will be shown that the introduction of the short-range correlation factor significantly improves the description of the dynamic correlation effects in latter sections.

3. INTEGRALS IN THE TRANSCORRELATED METHOD

The three-electron integrals, $\langle pqr|{\nabla }_1{f}_{12}\cdot {\nabla }_1{f}_{13}|stu\rangle$, are expressed in closed form in terms of recurrence relations [20]. The explicit treatment of the integrals, however, prohibits the application of the transcorrelated Hamiltonian to large molecular systems despite the use of the short-range geminal. The best feature of the present approach is that the three-body contribution is less important to make

Figure 1. The shape of the geminal, $f(r_{12})$, used in this work (full line). The dotted line denotes the linear $r_{12}$ behavior, required to balance the Coulomb repulsion around the origin.

Figure 1. The shape of the geminal, $f(r_{12})$, used in this work (full line). The dotted line denotes the linear $r_{12}$ behavior, required to balance the Coulomb repulsion around the origin.

the transcorrelated method feasible by using completeness insertions. One of such expressions is in the form,

$$\begin{gathered} \langle pqr|{\nabla }_1{f}_{12}\cdot {\nabla }_1{f}_{13}|stu\rangle ={\displaystyle \sum _{vwx}\langle pqr|{\nabla }_1{f}_{12}|vwx\rangle \cdot \langle vwx|{\nabla }_1{f}_{13}|stu\rangle } \\ ={\displaystyle \sum _v\langle pq|{\nabla }_1{f}_{12}|vt\rangle \cdot \langle vr|{\nabla }_1{f}_{12}|su\rangle }. \end{gathered}$$

(3.1)

The expression requires all Cartesian components of the integrals, $\langle pq|{\nabla }_1{f}_{12}|vt\rangle$, to which the point group symmetry cannot be applied. If the three basis functions, p, q and t, are s-orbitals with a coincident function center, only the p-component of the v function survives. This means that expansion functions, whose angular momentums are higher than the atomic orbitals by 1, are necessary in equation (3.1). Alternatively, we fall back on the idea of the operator identity,

$$L=\frac{1}{2}{[[H,F],F]}_s=\frac{1}{2}{[K^{(L)},F]}_s,$$

(3.2)

where $K^{(L)}$ and $F$ are given in the occupation number representation by

$$K^{(L)}=\frac{1}{2}{\displaystyle \sum _{pqrs}\langle pq|K_{12}^{(L)}|rs\rangle a_p^{+}{a}_q^{+}{a}_s{a}_r},$$

(3.3)

$$F=\frac{1}{2}{\displaystyle \sum _{pqrs}\langle pq|f_{12}|rs\rangle a_p^{+}{a}_q^{+}{a}_s{a}_r},$$

(3.4)

and the subscript, $s$, at the commutator denotes single contractions. Using the fact that $K_{12}^{(L)}$ is commutable with $f_{12}$ in the complete basis limit about the second electronic coordinate, we obtain the approximate expression for the three-electron integrals,

$$\begin{gathered} \langle pqr|{\nabla }_1{f}_{12}\cdot {\nabla }_1{f}_{13}|stu\rangle \\ =\frac{1}{2}{\displaystyle \sum _v(\langle pr|K_{12}^{(L)}|vu\rangle \langle vq|f_{12}|st\rangle -\langle pq|f_{12}|vt\rangle \langle vr|K_{12}^{(L)}|su\rangle )}. \end{gathered}$$

(3.5)

This expression features that the integrals, $\langle pq|f_{12}|rs\rangle$, are only required besides the anti-symmetric part of $K_{12}$. The point group symmetry is applicable to the integrals. Since the operator, $f_{12}$, does not increase the angular momentum of orbitals, the convergence of the completeness insertion of equation (3.5) is superior to the previous expression, equation (3.1). It is also possible to derive another expression using the averaged operator, ${\overline{K}}_{12}^{(L)}$, rather than $K_{12}^{(L)}$ as

$$\begin{gathered} \langle pqr|{\nabla }_1{f}_{12}\cdot {\nabla }_1{f}_{13}|stu\rangle \\ ={\displaystyle \sum _v(\langle pr|{\overline{K}}_{12}^{(L)}|vu\rangle \langle vq|f_{12}|st\rangle -\langle pq|f_{12}|vt\rangle \langle vr|{\overline{K}}_{12}^{(L)}|su\rangle )}. \end{gathered}$$

(3.6)

Both of the expressions with $K_{12}^{(L)}$ and ${\overline{K}}_{12}^{(L)}$ are exact in the complete basis limit. We examine the first in this work.

The present transcorrelated method requires the three different kinds of two-electron integrals, $\langle pq|K_{12}^{(L)}|rs\rangle$, $\langle pq|r_{12}^{-1}+K_{12}^{(Q)}|rs\rangle$ and $\langle pq|f_{12}|rs\rangle$. Assuming the basis functions to be primitive Gaussian-type functions and using the Laplace transform of the Coulomb operator, the electron repulsion integrals are written in the form,

$$\langle pq|r_{12}^{-1}|rs\rangle ={\displaystyle \sum _{l=0}^L{d}_l{g}_l^{(c)}},$$

(3.7)

$$g_l^{(c)}=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^{\infty }du{\Theta }^l(u^2)\langle 00|{\chi }_{12}(u^2)|00\rangle },$$

(3.8)

where $d_l$ are coefficients dependent only on the exponents and positions of the four basis functions, the zero indices means s-type functions, ${\chi }_{12}(\zeta )$ is the Gaussian function with the exponent, $\zeta$,

$${\chi }_{12}(\zeta )=\text{exp}(-\zeta r_{12}^2),$$

(3.9)

the function, $\Theta (\zeta )$, is given by

$$\Theta (\zeta )=\frac{\zeta }{\rho +\zeta },$$

(3.10)

and $\rho$ is a parameter which is a function of the four primitive exponents. The index, $l$, runs up to the sum of maximum angular momentums of the basis functions. For each point, $u$, the coefficients for the increments of angular momentums are independent in each Cartesian axis. The polynomial, ${\displaystyle \sum _{l=0}^L{d}_l{\Theta }^l}$, is therefore expressed as a product of the Cartesian contributions,

$${\displaystyle \sum _{l=0}^L{d}_l{\Theta }^l}=I_x(\Theta )I_y(\Theta )I_z(\Theta ).$$

(3.11)

The Cartesian components, $I_k(\Theta ),(k=x,y,z)$, can be calculated using recurrence relations [22].

Similarly, the integrals, $\langle pq|f_{12}|rs\rangle$, are expressed as,

$$\langle pq|f_{12}|rs\rangle ={\displaystyle \sum _{l=0}^L{d}_l}{g}_l^{(f)},$$

(3.12)

$$g_l^{(f)}=-{\displaystyle \sum _G^{N_G}{c}_G{\Theta }^l({\zeta }_G)\langle 00|{\chi }_{12}({\zeta }_G)|00\rangle }.$$

(3.13)

The recurrence relation [22] shows that the integrals for the operators, $-{\nabla }_1^2{f}_{12}$ and $-{\nabla }_1{f}_{12}\cdot {\nabla }_1$, which are more complex, also reduce to the kind of integrals, $\langle pq|f_{12}|rs\rangle$. The maximum value, $L$, increases by 2 because of the differentiations. For the operator, $K_{12}^Q$, we introduce the integrals,

$$\langle pq|q_{12}|rs\rangle ={\displaystyle \sum _{l=0}^L{d}_l}{g}_l^{(q)}$$

(3.14)

$$g_l^{(q)}=-{\displaystyle \sum _G^{N_G}{\displaystyle \sum _{\overline{G}}^{N_G}{c}_G{c}_{\overline{G}}\frac{{\zeta }_G{\zeta }_{\overline{G}}}{Z_{G\overline{G}}}{\Theta }^l(Z_{G\overline{G}})\langle 00|{\chi }_{12}(Z_{G\overline{G}})|00\rangle }},$$

(3.15)

where $Z_{G\overline{G}}$ denote sums of exponents of the geminal, $Z_{G\overline{G}}={\zeta }_G+{\zeta }_{\overline{G}}$. Using the relations,

$${\nabla }_1{\chi }_{12}({\zeta }_G)\cdot {\nabla }_1{\chi }_{12}({\zeta }_{\overline{G}})=\frac{{\zeta }_G{\zeta }_{\overline{G}}}{Z_{G\overline{G}}}(\frac{{\nabla }_1^2}{Z_{G\overline{G}}}-6){\chi }_{12}(Z_{G\overline{G}}),$$

(3.16)

$$1-\frac{Z_{G\overline{G}}}{\rho +Z_{G\overline{G}}}=\frac{\rho }{Z_{G\overline{G}}}\frac{Z_{G\overline{G}}}{\rho +Z_{G\overline{G}}},$$

(3.17)

it is shown that the integrals for $K_{12}^Q$ reduce to the combination of the type of the integrals, $\langle pq|q_{12}|rs\rangle$. One notes that the coefficients, $d_l$, which are dependent only on the angular momentums and the positions of the orbital quartet, are shared by all types of the integrals.

We express the functions, $g_l^{(x)}(x=c,f,g)$, using polynomials as

$$g_l^{(x)}={\displaystyle \sum _{\alpha =0}^L{R}_{\alpha }^l{W}_{\alpha }^{(\kappa )}},(\kappa =c,f,q),$$

(3.18)

where $R_{\alpha }$ and $W_{\alpha }^{(\kappa )}$ are the root and the weight for the point, $\alpha$, which runs over the same range as $l$. The functions, $g_l^{(\kappa )}$ and $W_{\alpha }^{(\kappa )}$, are connected through a linear transformation predetermined for $R_{\alpha }$. In the Rys quadrature method [23,24], both of the roots and weights are chosen to be variables reducing the number of root points to be larger than $L/2$ rather than $L$. We do not employ this convention but use fixed roots for all operators, $\kappa =c,f,q$. Once the integrals are written in the discrete polynomials, the summation over $l$ is expressed as a product of Cartesian components. For instance, the electron repulsion integrals become in the form,

$$\langle pq|r_{12}^{-1}|rs\rangle ={\displaystyle \sum _{\alpha =0}^L{W}_{\alpha }^{(c)}{\displaystyle \sum _{l=0}^L{d}_l{R}_{\alpha }^l}}={\displaystyle \sum _{\alpha =0}^L{W}_{\alpha }^{(c)}{I}_x(R_{\alpha })I_y(R_{\alpha })I_z(R_{\alpha })}.$$

(3.19)

The above procedure enables us to calculate the integrals for different two-electronic operators simultaneously by changing the weight, $W_{\alpha }^{(\kappa )}$, after setting up the common quantities, $I_k(\Theta ),(k=x,y,z)$.

4. APPLICATION TO SINGLE-REFERENCE MANY ELECTRON THEORIES

Basically, the application of the transcorrelated method is straightforward; the original Hamiltonian, $H$, is replaced by the transcorrelated one, $\tilde{H}$, in the usual ab initio theories. The transcorrelated Hamiltonian is, however, nonhermite including weak three-body interactions. This property gives rise to a few variations in formulating single reference theories. Throughout this paper, we use the notations, $a,b,\dots$ and $i,j,\dots$ for virtual and occupied orbitals within a given basis set, respectively. Our first method employs a coincident vacuum determinant for the bra and ket states in the CC equations,

$$\langle {\psi }_0|\tilde{H}\text{exp}(T)|{\psi }_0\rangle =E,$$

(4.1)

$$\langle X_{\mu }^{+}{\psi }_0|\tilde{H}\text{exp}(T){|{\psi }_0\rangle }_c=0,\forall \mu \in S,D,\dots ,$$

(4.2)

where ${\psi }_0$ is a single-determinant vacuum, $\left\{X_{\mu }^{+}\right\}=\{a_a^{+}{a}_b^{+}\dots a_j{a}_i\}$ denotes particle-hole excitation operators with respect to the vacuum, and $c$ means the connected entity of the operators. A suitable partitioning of the Hamiltonian leads to an order-by-order expansion of the wave operator. Since our transcorrelated Hamiltonian is designed to coincide with the long-range behavior of the original Coulomb operator, the usual Hartree-Fock (HF) model Hamiltonian in the MP perturbation theory is a good starting point for the expansion [19],

$$\tilde{H}=H_0^{(MP)}+\tilde{V},$$

(4.3)

choosing ${\psi }_0$ to be the HF determinant for the original Hamiltonian. One notes that the new perturbation, $\tilde{V}$, is nearly free from the Coulomb singularity. We therefore can anticipate that the perturbational series converges quicker than the usual MP one.

The above argument is not the case for electrons around nuclei especially when we introduce an approximation to the CC method. The operator dependent on the momentum of the anti-symmetric part, $-{\nabla }_1{f}_{12}\cdot {\nabla }_1$, becomes significant for core electrons of heavy atomic elements. In the present application, we preliminarily treat the orbital relaxation with the pseudo-orbital equation, excluding the expensive but less important three-body contribution from the iterative loop,

$$\langle X_{\mu }^{+}{\psi }_0|(H+K)\text{exp}({\overline{T}}_1){|{\psi }_0\rangle }_c=0,\forall \mu \in S.$$

(4.4)

The subsequent connected singles and doubles are restored order-by-order approximately using the LCC equation known to be CEPA0 [25,26],

$$\begin{gathered} \langle X_{\mu }^{+}{\psi }_0|[H_0^{(MP)},\Delta T^{(n+1)}]|{\psi }_0\rangle \\ =-\langle X_{\mu }^{+}{\psi }_0|\widehat{\tilde{H}}+[\widehat{\tilde{H}},T^{(n)}]|{\psi }_0\rangle ,\forall \mu \in S,D, \end{gathered}$$

(4.5)

The similarity transformed Hamiltonian,

$$\widehat{\tilde{H}}=\text{exp}(-{\overline{T}}_1)\tilde{H}\text{exp}({\overline{T}}_1),$$

(4.6)

analogous to the one, which appeared in the treatment of CCSD amplitudes [27,28], retains the operator rank of the transcorrelated Hamiltonian. The procedure is adequate especially for a molecule consisting of light atomic elements.

Another option, which gives results similar to those of the above-mentioned modified MP partitioning with the pseudo-orbital equation but uses more optimized reference functions, is based on biorthogonal one-electronic basis sets [29,30]. For a while, we denote the orbitals in the one-electronic sets as $\{{\chi }_1,{\chi }_2,\dots \}$ and $\{{\phi }_1,{\phi }_2,\dots \}$, respectively. The biorthogonality condition is given by

$$\langle {\chi }_p|{\phi }_q\rangle ={\delta }_{pq}.$$

(4.7)

Creation and annihilation operators are defined using the field operators as

$${\Psi }^{+}(1)={\displaystyle \sum _p{\chi }_p(1)b_p^{+}},$$

(4.8)

$$\Psi (1)={\displaystyle \sum _q{\phi }_q(1)c_q},$$

(4.9)

Operators are expressed in the second quantized form of the biorthogonal set, like

$$K=\frac{1}{2}{\displaystyle \sum _{pqrs}\langle {\chi }_p{\chi }_q|K_{12}|{\phi }_r{\phi }_s\rangle b_p^{+}{b}_q^{+}{c}_s{c}_r}.$$

(4.10)

Since the anti-commutation relation holds for the new creation and annihilation operators,

$$\{b_p^{+},c_q\}={\delta }_{pq},$$

(4.11)

the Wick theorem and accordingly the usual diagrammatic technique can be used for manipulating operators in the biorthogonal form.

We introduce separate single-reference vacua consisting of $n$ biorthogonal orbitals, ${\rm X}=\Vert {\chi }_1{\chi }_2\dots {\chi }_n\Vert$ and $\Phi =\Vert {\phi }_1{\phi }_2\dots {\phi }_n\Vert$. The canonical orbitals are determined self-consistently using the approximate Fock operator [21],

$$\langle {\chi }_p|\tilde{f}|{\phi }_q\rangle =\langle {\chi }_p|h|{\phi }_q\rangle +{\displaystyle \sum _i\langle {\chi }_p{\chi }_i|r_{12}^{-1}+K_{12}{|{\phi }_q{\phi }_i\rangle }_A}={\delta }_{pq}{\epsilon }_p,$$

(4.12)

where $A$ means anti-symmetrized integrals. We solve the LCCSD equation,

$$\begin{gathered} \langle Z_{\mu }^{+}{\rm X}|[H_0^{(BO)},\Delta T_b^{(n+1)}]|\Phi \rangle \\ =-\langle Z_{\mu }^{+}{\rm X}|\tilde{H}+[\tilde{H},T_b^{(n)}]|\Phi \rangle ,\forall \mu \in S,D, \end{gathered}$$

(4.13)

iteratively with the zeroth order Hamiltonian, $H_0^{(BO)}={\displaystyle \sum _p{b}_p^{+}{c}_p{\epsilon }_p}$, where we define the excitation and deexciation operators with respect to the right vacuum, $\left\{Y_{\mu }^{+}\right\}=\{b_a^{+}{b}_b^{+}\dots c_j{c}_i\}$ and $\left\{Z_{\mu }\right\}=\{b_i^{+}{b}_j^{+}\dots c_b{c}_a\}$, and the cluster operator, $T_b={\displaystyle \sum _{\mu }{t}_{\mu }{Y}_{\mu }^{+}}$.

It should be noted that the present methods do not alter the original scaling of the correlated methods. The pseudo-orbital and biothogonal self-consistent field (SCF) methods, which exclude the contribution of $L$, are driven in the atomic orbital basis. This gives rise to the formal $N^4$ scaling as the usual SCF method. The subsequent LCCSD code is commonly used in both of the methods substituting the individual transformed integrals and energy denominators. The double contractions of $T_2$ and the two-electron operators in the transcorrelated Hamiltonian, $r_{12}^{-1}+K_{12}^{(Q)}$, $K_{12}^{(L)}$ and $f_{12}$, become the rate-determining steps at this level. These scale as at most $V^4{O}^2$ as the conventional CCSD and LCCSD methods for $V$ and $O$ denoting the numbers of virtual and occupied orbitals, respectively. The effect of the connected triples from $L$ is negligibly small owing to the localized nature of the correlation factor.

5. RESULT AND DISCUSSION

We examine the efficiency of the transcorrelated method in the calculations of the typical 10-electronic systems, CH₃, NH₃, H₂O, HF and Ne. We first check the convergence of the correlation energies for H₂O and Ne with the correlation consistent basis sets of Dunning [31], cc-pVXZ (X = D, T, Q and 5). The basis functions are used as primitives to increase the flexibility. Each basis set is further augmented with the d- f- and g-core functions in the corresponding cc-pCXZ sets [32], excluding the primitives in the range of the exponents in the parent set. This gives rise to the $9s4p1d$, $10s5p3d1f$, $12s6p5d3f1g$ and $14s8p7d5f3g$ sets. The corresponding primitive sets derived from cc-pVXZ are used for hydrogen with the exception that cc-pVQZ is used in the largest set for main elements. We approximate the connected double excitations correct to the first order from the secondary perturbation of the operator, $L$, in this particular work.

Table I shows the result for H₂O. At the second order level, the results are relatively sensitive to the choice of the reference functions; the differences of the correlation energies are as large as $7m{E}_h$ between the modified MP and biorthogonal theories. However, they decrease as the order of perturbation increases. Finally the differences converge to the LCCSD values, ca. $0.5m{E}_h$, with all of the basis sets. The basis set dependencies of the MBPT(2) and LCCSD energies are depicted in Fig. 2. For comparison, we show the R12-MBPT(2), CCSD and CCSD(T) energies for H₂O with 374 basis functions [33]. The difference between the conventional and transcorrelated LCCSD energies is $60m{E}_h$ with the $9s4p1d$ set. This means that the use of the present methods are particularly advantageous with such a small basis set. The energies of the transcorrelated methods seem almost

Table I. Basis set dependence of the correlation energies of water ^a)

Basis Set	Orderb)	Modified MP	Biorthogonal	Conventional
$9s4p1d$	2nd	325.38	317.76	257.74
	3rd	323.40	322.82	263.61
	4th	328.08	327.38	268.58
	$\infty$	329.34	329.02	270.20
$10s5p3d1f$	2nd	354.54	347.34	318.46
	3rd	352.19	351.37	323.38
	4th	357.78	356.97	329.65
	$\infty$	358.95	358.53	331.14
$12s6p5d3f1g$	2nd	362.20	355.47	341.90
	3rd	358.86	357.89	344.43
	4th	364.70	363.89	351.14
	$\infty$	365.85	365.38	352.55
$14s8p7d5f3g$	2nd	363.11	356.54	348.89
	3rd	359.24	358.17	350.14
	4th	365.18	364.38	357.11
	$\infty$	366.34	365.87	358.50

^a) In -mE_h unit.

^b) Order of the perturbation in the correlation energy within LCCSD.

^c) Nearly exact correlation energies are -361.486, -359.775 and -369.89 mE_h at the R12-MBPT(2)-A, R12-CCSD and R12-CCSD(T) levels of approximations with 374 Gaussian-type functions.

Table I. Basis set dependence of the correlation energies of water ^a)

Basis Set	Orderb)	Modified MP	Biorthogonal	Conventional
$9s4p1d$	2nd	325.38	317.76	257.74
	3rd	323.40	322.82	263.61
	4th	328.08	327.38	268.58
	$\infty$	329.34	329.02	270.20
$10s5p3d1f$	2nd	354.54	347.34	318.46
	3rd	352.19	351.37	323.38
	4th	357.78	356.97	329.65
	$\infty$	358.95	358.53	331.14
$12s6p5d3f1g$	2nd	362.20	355.47	341.90
	3rd	358.86	357.89	344.43
	4th	364.70	363.89	351.14
	$\infty$	365.85	365.38	352.55
$14s8p7d5f3g$	2nd	363.11	356.54	348.89
	3rd	359.24	358.17	350.14
	4th	365.18	364.38	357.11
	$\infty$	366.34	365.87	358.50

^a) In -mE_h unit.

^b) Order of the perturbation in the correlation energy within LCCSD.

^c) Nearly exact correlation energies are -361.486, -359.775 and -369.89 mE_h at the R12-MBPT(2)-A, R12-CCSD and R12-CCSD(T) levels of approximations with 374 Gaussian-type functions.

Figure 2. Convergence of the correlation energy of H₂O. The abbreviated notations, BO and MMP, denote the biorthogonal and modified MP methods, respectively.

Figure 2. Convergence of the correlation energy of H₂O. The abbreviated notations, BO and MMP, denote the biorthogonal and modified MP methods, respectively.

saturated in each order with the $12s6p5d3f1g$ set, while the conventional methods still gain the correlation energies of $5-7m{E}_h$ from $12s6p5d3f1g$ to $14s8p7d5f3g$ set. The absolute correlation energies of the transcorrelated LCCSD methods with the largest set are $6m{E}_h$ as low as that of R12-CCSD. The usual situation that the fourth order quadruples diagrams contribute positive [34] holds in the transcorrelated methods.

Let us proceed to the result of Ne in Table II. In comparison with the H₂O result, the differences between the modified MP and biorthogonal methods are larger by the amounts, ca. $10m{E}_h$ and $3m{E}_h$ at the MBPT(2) and LCCSD levels, respectively. Especially, the modifies MP series tends to overshoot the correlation energy for the system. The method is considered to become less adequate for Ne because of the increasing magnitude of $-{\nabla }_1{f}_{12}\cdot {\nabla }_1$ especially for the description of correlations including core electrons. The inclusions of higher-order corrections will be required to recover the difference. The convergence behaviors are shown in Fig. 3. In spite of the increase of the discrepancy between the two biorthogonal methods, the basis set dependence is similar to that for H₂O beyond the $10s5p3d1f$ set; the differences are sufficiently small within $1m{E}_h$ between the $12s6p5d3f1g$ to $14s8p7d5f3g$ sets. With the small sets, $9s4p1d$ and $10s5p3d1f$, the transcorrelated method of the biorthogonal reference reproduces 95.4% and 98.8% of the saturated energies at the second order

Table II. Basis set dependence of the correlation energies of neon^a)

Basis Set	Order b)	Modified MP	Biorthogonal	Conventional
$9s4p1d$	2nd	382.54	368.19	255.48
	3rd	371.75	366.67	255.97
	4th	371.73	368.91	258.20
	$\infty$	372.19	369.20	258.56
$10s5p3d1f$	2nd	392.39	381.17	330.53
	3rd	386.18	379.36	331.79
	4th	386.16	383.06	335.18
	$\infty$	386.60	383.23	335.54
$12s6p5d3f1g$	2nd	395.40	385.93	360.63
	3rd	391.21	383.65	361.12
	4th	391.28	388.09	364.93
	$\infty$	391.72	388.15	365.27
$14s8p7d5f3g$	2nd	394.62	385.94	370.69
	3rd	391.00	383.08	370.13
	4th	391.09	387.91	374.13
	$\infty$	391.52	387.89	374.46

^a) In -mE_h unit.

^b) Order of the perturbation in the correlation energy within LCCSD.

^c) Nearly exact correlation energies are -388.311, -383.823 and -390.508 mE_h at the R12-MBPT(2)-A, R12-CCSD and R12-CCSD(T) levels of approximations with 293 Gaussian-type functions.

Table II. Basis set dependence of the correlation energies of neon^a)

Basis Set	Order b)	Modified MP	Biorthogonal	Conventional
$9s4p1d$	2nd	382.54	368.19	255.48
	3rd	371.75	366.67	255.97
	4th	371.73	368.91	258.20
	$\infty$	372.19	369.20	258.56
$10s5p3d1f$	2nd	392.39	381.17	330.53
	3rd	386.18	379.36	331.79
	4th	386.16	383.06	335.18
	$\infty$	386.60	383.23	335.54
$12s6p5d3f1g$	2nd	395.40	385.93	360.63
	3rd	391.21	383.65	361.12
	4th	391.28	388.09	364.93
	$\infty$	391.72	388.15	365.27
$14s8p7d5f3g$	2nd	394.62	385.94	370.69
	3rd	391.00	383.08	370.13
	4th	391.09	387.91	374.13
	$\infty$	391.52	387.89	374.46

^a) In -mE_h unit.

^b) Order of the perturbation in the correlation energy within LCCSD.

^c) Nearly exact correlation energies are -388.311, -383.823 and -390.508 mE_h at the R12-MBPT(2)-A, R12-CCSD and R12-CCSD(T) levels of approximations with 293 Gaussian-type functions.

Figure 3. Convergence of the correlation energy of Ne. The notations are same as those in Figure 2.

Figure 3. Convergence of the correlation energy of Ne. The notations are same as those in Figure 2.

level, while the conventional method covers only 65.8% and 85.1% of the R12-MBPT(2)-A energy. Although the modified MP series slightly overshoots the correlation energies, the almost saturated biorthogonal LCCSD energy is in between the R12-CCSD and R12-CCSD(T) ones consistently with the previous results for H₂O.

We show the results of the entire 10-electronic systems with the $12s6p5d3f1g$ set in Table III. The biorthogonal SCF and pseudo orbital energies, $\langle {\rm X}|\tilde{H}|\Phi \rangle$ and $\langle {\psi }_0|\widehat{\tilde{H}}|{\psi }_0\rangle$, differ from each other by $5-40m{E}_h$ increasingly with the mass of the main elements. This is relevant to the descriptions concerning 1s core electrons as mentioned previously. The differences, however, reduce to $0.2-3.5m{E}_h$ at the LCCSD levels. The biorthogonal LCCSD results are always in the range between R12-CCSD and R12-CCSD(T). Especially the energies are consistently below the R12-CCSD energies by $5-7m{E}_h$, implying that the results are reliable in the basis set convergence. For further improvements of the absolute energies, inclusion of the primary fourth order diagrams of triple- and quadruple-excitations will be necessary.

6. CONCLUSION

We have seen that the use of the transcorrelated Hamiltonian is particularly effective for accelerating the convergence of correlation energies with the size of one electronic basis. The key features of the method are in the following. We employ the symmetric correlation factor, which is commutable with the Coulomb potential. This leads to the simple structure of the transcorrelated Hamiltonian. The explicitly correlated functions are used to deal with the short-range behavior of the many-electron wave functions. Instead, the residual correlation effects, which are long-ranged, are treated in terms of the usual CI-type expansion with one-electronic basis. Thus the number of the additional integrals

Table III. Correlation energies of Ne, HF, H₂O, NH₃ and CH₄ molecules with the [12s6p5d3f1g/5s3p2d1f] basis set^a)

	Ne	HF	H2O	NH3	CH4
HF	128543.47	100067.71	76064.89	56223.14	40216.34
Conv. MBPT(2)	360.63	360.58	341.90	307.07	260.33
Conv. LCCSD	365.27	366.36	352.55	325.52	287.61
Biorthogonal
SCF	184.14	138.06	102.13	74.70	53.77
2nd order	385.93	379.11	355.47	316.74	267.71
3rd order	383.65	376.99	357.89	326.77	285.68
4th order	388.09	382.95	363.89	332.56	291.82
LCCSD	388.15	383.63	365.38	334.62	294.39
Modified MP
Pseudo Orbital	222.47	158.84	112.55	79.55	55.84
2nd order	395.40	387.80	362.20	321.14	270.03
3rd order	391.21	381.04	358.86	326.33	285.11
4th order	391.28	384.66	364.70	332.80	291.70
LCCSD	391.72	385.47	365.85	334.41	294.02
R12b)
MBPT(2)-A	387.910	384.357	361.691	322.400	273.579
CCSD	382.945	378.392	359.312	327.829	288.557
CCSD(T)	389.456	387.260	369.228	337.248	295.948

^a) In -mE_h unit. ^b) The correlation energies with $[15s9p6d5f4g/6s3p2d]$ basis set.

Table III. Correlation energies of Ne, HF, H₂O, NH₃ and CH₄ molecules with the [12s6p5d3f1g/5s3p2d1f] basis set^a)

	Ne	HF	H2O	NH3	CH4
HF	128543.47	100067.71	76064.89	56223.14	40216.34
Conv. MBPT(2)	360.63	360.58	341.90	307.07	260.33
Conv. LCCSD	365.27	366.36	352.55	325.52	287.61
Biorthogonal
SCF	184.14	138.06	102.13	74.70	53.77
2nd order	385.93	379.11	355.47	316.74	267.71
3rd order	383.65	376.99	357.89	326.77	285.68
4th order	388.09	382.95	363.89	332.56	291.82
LCCSD	388.15	383.63	365.38	334.62	294.39
Modified MP
Pseudo Orbital	222.47	158.84	112.55	79.55	55.84
2nd order	395.40	387.80	362.20	321.14	270.03
3rd order	391.21	381.04	358.86	326.33	285.11
4th order	391.28	384.66	364.70	332.80	291.70
LCCSD	391.72	385.47	365.85	334.41	294.02
R12b)
MBPT(2)-A	387.910	384.357	361.691	322.400	273.579
CCSD	382.945	378.392	359.312	327.829	288.557
CCSD(T)	389.456	387.260	369.228	337.248	295.948

^a) In -mE_h unit. ^b) The correlation energies with $[15s9p6d5f4g/6s3p2d]$ basis set.

grows linearly to the system size, making the application to large molecular systems feasible. The short-range nature of the correlation factor is also effective to make the three-body contributions less important. This property allows us to use the resolution of identity appropriately for the quantities less significant than the ones like $r_{12}^{-1}{f}_{13}$, which appear in other theories. As a result, the original scaling of a correlated method is maintained even after the introduction of the transcorrelated Hamiltonian.

There are mass dependences of the reference energies of the transcorrelated Hamiltonian, $\langle {\rm X}|\tilde{H}|\Phi \rangle$ and $\langle {\psi }_0|\widehat{\tilde{H}}|{\psi }_0\rangle$, because the operator, $-{\nabla }_1{f}_{12}\cdot {\nabla }_1$, brings more significant contribution on the core electrons as the atomic mass increases. The correlation behaviors are different between the core and valence shells. The original transcorrelated method as well as most of the GTG methods use explicitly correlated functions, which are dependent on the positions of electrons. The introduction of such position dependence will enable more balanced treatment of different shells though the evaluation of integrals becomes complicated. The use of the spherically symmetric geminal is much more feasible. The present result implies that the ignorance of the position dependence does not spoil the accuracy at the LCCSD level to say the least for elements as heavy as Ne.

For calculations of excited states and bond breakings, multireference treatment is necessary. Such an application of the transcorrelated Hamiltonian is also straightforward. Extensions in this line are now under progress. Another application will be the calculation of static structure factors [35,36]. Since this quantity is a probe sensitive to the dynamic correlation effects, comparisons with experimental data will be an important measure of the explicitly correlated method.