Multicenter Molecular Integrals over Dirac Wave Functions for Several Fundamental Properties

Abstract

Multicenter molecular integrals over Dirac wave functions can be derived by using the Gaussian-transform for the Dirac wave function, which was derived by the author, for several fundamental properties; i.e., the overlap integral, the kinetic energy one, the nuclear attraction one for the point-like nucleus and for the finite one, and the electron-repulsion integral.

1. Introduction

Recently, Sun, et al. [1] derived the gauge invariant Dirac equation given by

$\left(\begin{array}{cc}{m}_e{c}^2+V& c\overrightarrow{\sigma }\cdot (\overrightarrow{P}+\overrightarrow{A})\\ c\overrightarrow{\sigma }\cdot (\overrightarrow{P}+\overrightarrow{A})& -m_e{c}^2+V\end{array}\right)\left(\begin{array}{c}{\Psi }^L\\ {\Psi }^S\end{array}\right)=\left(\begin{array}{c}{\Psi }^L\\ {\Psi }^S\end{array}\right)E_0$ (1.1)

where m_e is the electron rest mass, c is the speed of light, V is the scalar potential, $\overrightarrow{\sigma }$ is the Pauli spin matrices, $\overrightarrow{p}=-i\hslash \nabla$ is the momentum, $\overrightarrow{A}$ is the vector potential of the magnetic field due to the nuclear spin, ψL is the large component spinor, ψs is the small component spinor, and E₀ is the energy. We use the atomic units throughout the present article ( $m_e=1,e=1,\hslash =1,4\pi {\epsilon }_0=1,c=\text{137}.0\text{35999139}$ ). However, we describe m_e, e, and $\hslash$ explicitly, for the readers' convenience when one converts the units to the natural units. We subtract the rest- mass energy m_ec^e from E₀ to align the energy scale to that of the Schrödinger equation. So Eq. (1.1) can be modified to

$\left(\begin{array}{cc}V& c\overrightarrow{\sigma }\cdot (\overrightarrow{P}+\overrightarrow{A})\\ c\overrightarrow{\sigma }\cdot (\overrightarrow{P}+\overrightarrow{A})& -2m_e{c}^2+V\end{array}\right)\left(\begin{array}{c}{\Psi }^L\\ {\Psi }^S\end{array}\right)=\left(\begin{array}{c}{\Psi }^L\\ {\Psi }^S\end{array}\right)E$ (1.2)

Where E = E₀ − m_ec^e. To solve this Dirac equation, we may use a suitable basis set, {χ}. The large component spinor can be expressed as a linear combination in terms of these basis functions as given by

${\Psi }_i^L={\displaystyle {\sum }_{\mu }{C}_{i\mu }^L{\chi }_{\mu }}$ (1.3)

However, the small component spinor is in the variational collapse until using the restricted magnetic balance (RMB) [2] as given by

${\Psi }_i^S={\displaystyle {\sum }_{\mu }{C}_{i\mu }^S}\overrightarrow{\sigma }\cdot (\overrightarrow{P}+\overrightarrow{A}){\chi }_{\mu }$ (1.4)

Recently, Yoshizawa [3] derived the matrix Dirac equation using the RMB as given by

$\left(\begin{array}{cc}\overleftrightarrow{V}& \overleftrightarrow{T_m}\\ \overleftrightarrow{T_m}& \overleftrightarrow{W_m}-\overleftrightarrow{T_m}\end{array}\right)\left(\begin{array}{cc}\overleftrightarrow{C_{-}^L}& \overleftrightarrow{C_{+}^L}\\ \overleftrightarrow{C_{-}^S}& \overleftrightarrow{C_{+}^S}\end{array}\right)=\left(\begin{array}{cc}\overleftrightarrow{S}& \overleftrightarrow{0}\\ \overleftrightarrow{0}& \frac{1}{2m_e{c}^2}\overleftrightarrow{T_m}\end{array}\right)\left(\begin{array}{cc}\overleftrightarrow{C_{-}^L}& \overleftrightarrow{C_{+}^L}\\ \overleftrightarrow{C_{-}^S}& \overleftrightarrow{C_{+}^S}\end{array}\right)\left(\begin{array}{cc}\overleftrightarrow{{ϵ}_{-}}& \overleftrightarrow{0}\\ \overleftrightarrow{0}& \overleftrightarrow{{ϵ}_{+}}\end{array}\right)$ (1.5)

where $\overleftrightarrow{C_{-}^L}$ is the coefficient matrix of the large component spinor for the energy ∈−, $\overleftrightarrow{C_{+}^L}$ is that for ∈+, $\overleftrightarrow{C_{-}^S}$ and $\overleftrightarrow{C_{+}^S}$ are those of the small component spinor, $\overleftrightarrow{{ϵ}_{-}}$ and $\overleftrightarrow{{ϵ}_{+}}$ are the energy matrices, $\overleftrightarrow{0}$ is the zero matrix,

$V_{\mu \nu }=<{\chi }_{\mu }\left|V\right|{\chi }_{\nu }>$ (1.6)

${(T_m)}_{\mu \nu }=\frac{1}{2m_e}<{\chi }_{\mu }\left|\overrightarrow{\sigma }\cdot \left(\overrightarrow{P}+\overrightarrow{A}\right)\overrightarrow{\sigma }\cdot \left(\overrightarrow{P}+\overrightarrow{A}\right)\right|{\chi }_{\nu }>$ (1.7)

${(W_m)}_{\mu \nu }=\frac{1}{4m_e^2{c}^2}<{\chi }_{\mu }\left|\overrightarrow{\sigma }\cdot \left(\overrightarrow{P}+\overrightarrow{A}\right)V\overrightarrow{\sigma }\cdot \left(\overrightarrow{P}+\overrightarrow{A}\right)\right|{\chi }_{\nu }>$ (1.8)

and

$S_{\mu \nu }=<{\chi }_{\mu }|{\chi }_{\nu }>$ (1.9)

Many researchers extend the matrix Dirac equation to the molecule [2-17]. Especially, many are for relativistic calculations of NMR spectra [2,3,13-17]. It is natural to use the atomic Dirac wave function as the function among the basis set in order to solve the molecular matrix Dirac equation. However, it has not been used yet, because there are no molecular integral formulas. In a previous article [18], the author derived the Gaussian-transform for the Dirac wave function centered at A given by

$r_A^{-{\epsilon }_A}\exp \left(-\zeta r_A\right)=\frac{{\zeta }_A^{1+{\epsilon }_A}}{\Gamma \left(1+{\epsilon }_A\right)2{\pi }^{1/2}}{\displaystyle {\int }_0^{\infty }ds{s}^{-3/2}}\text{exp}[-s{r}_A^2]\left[\frac{{\zeta }_A^2}{2s}{\displaystyle {\int }_0^{1}dt\frac{{(1-t)}^{{\epsilon }_A}}{t^{4+{\epsilon }_A}}}-{\displaystyle {\int }_0^{1}dt\frac{{(1-t)}^{{\epsilon }_A}}{t^{2+{\epsilon }_A}}}\right]\text{exp}\left[-\frac{{\zeta }_A^2}{4s{t}^2}\right]$ (1.10)

where ${\epsilon }_A=1-\sqrt{1-{(Z_A\alpha )}^2}$ and ζ_A = Z_A in which α = 1/137.035999139 is the fine structure constant and $Z_{A}e$ is the nuclear charge. Equation (1.10) is the only formula to be able to evaluate the multicenter molecular integral over Dirac wave functions. No study treats the Dirac wave function in molecular systems. The author’s study is the first time to treat the Dirac wave function in molecular systems. Anyone can derive desired molecular integral by using Eq. (1.10). In the present article, we derive multicenter molecular integrals over Dirac wave functions for several fundamental properties; i.e., the overlap integral in the next section, the Kinetic Energy Integral (KEI) in the third section, the Nuclear Attraction Integral (NAI) for both the point-like nucleus and the finite one in the fourth section, and the Electron-Repulsion Integral (ERI) in the fifth section.

2. Overlap integral

The two-center overlap integral over Dirac wave functions can be given by

$OVL=N_A{N}_B{S}_{AB}$ (2.1)

where $N_A=\sqrt{\frac{{(2{\zeta }_A)}^{3-2{\epsilon }_A}}{4\pi \Gamma (3-2{\epsilon }_A)}}$ is the normalization constant and

$S_{AB}={{\displaystyle \int }}^{\text{}}d\overrightarrow{r}{r}_A^{-{\epsilon }_A}{r}_B^{-{\epsilon }_B}\text{exp}[-{\zeta }_A{r}_A-{\zeta }_B{r}_B]$ (2.2)

Using Eq. (1.10), we have

$S_{AB}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}}{4\pi \Gamma \left(1+{\epsilon }_A\right)\Gamma \left(1+{\epsilon }_B\right)}{\displaystyle {\int }_0^{\infty }d{S}_1}{\displaystyle {\int }_0^{\infty }d{S}_2}{(S_1{S}_2)}^{－３/2}exp\left[-S_1{r}_A^2-S_2{r}_B^2\right]\left[\frac{{\zeta }_A^2}{2S_1}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4S_1{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2S_2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4S_2{t}_2^2}\right]I_1^{(2)}$ (2.3)

where

$I_1^{(2)}={{\displaystyle \int }}^{\text{}}d\overrightarrow{r}\text{exp}[-S_{12}{r}_P^2]=\frac{{\pi }^{3/2}}{s_{12}^{3/2}}$ (2.4)

in which S₁₂ = S₁ + S₂. In the above derivation, we use the Gaussian product rule given by

$\exp \left[-S_1{r}_A^2-S_2{r}_B^2\right]=exp\left[-\frac{S_1{S}_2}{S_{12}}{\overline{AB}}^2-S_{12}{r}_P^2\right]$ (2.5)

where $\overrightarrow{P}=\frac{S_1}{S_{12}}\overrightarrow{A}+\frac{S_2}{S_{12}}\overrightarrow{B}$ . Let us change the integral variables as S₁₂ =z and $\frac{S_1}{S_{12}}=w$ . The Jacobian is $\frac{\partial (S_1,S_2)}{\partial (w,z)}=z$ . Thus, we have

$S_{AB}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}}{\Gamma \left(1+{\epsilon }_A\right)\Gamma \left(1+{\epsilon }_B\right)}\frac{\sqrt{\pi }}{4}{\displaystyle {\int }_0^{1}dw{\displaystyle {\int }_0^{\infty }dz}}{[w(1-w)]}^{-\frac{3}{2}}{z}^{-\frac{7}{2}}\text{exp}[-w\left(1-w\right)z{\overline{AB}}^2]\left[\frac{{\zeta }_A^2}{2wz}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4wz{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2\left(1-w\right)z}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4\left(1-w\right)z{t}_2^2}\right]$ (2.6)

We separate the integral over z as follows:

${\displaystyle {\int }_0^{\infty }dz}={\displaystyle {\int }_0^{a^2}dz}+{\displaystyle {\int }_{a^2}^{\infty }dz}$ (2.7)

Where a² can be chosen arbitrarily. We choose as a² = 4, here. Next, we change the integral variables as follows: In the first term of the right-hand side in Eq. (2.7), we do as z = ua². We do as $z=\frac{a^2}{u}$ in the last term in Eq. (2.7). Thus, we have

${\displaystyle {\int }_0^{\infty }dz}=a^2{\displaystyle {\int }_0^{1}du}+a^2{\displaystyle {\int }_0^{1}du}\frac{1}{u^2}$ (2.8)

Substituting Eq. (2.8) into Eq. (2.6), we have the final formula for the overlap integral over Dirac wave functions given by

$S_{AB}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}}{\Gamma \left(1+{\epsilon }_A\right)\Gamma \left(1+{\epsilon }_B\right)}\frac{\sqrt{\pi }}{4a^5}{\displaystyle {\int }_0^{1}dw}{[w(1-w)]}^{-\frac{3}{2}}\{{\displaystyle {\int }_0^{1}du}\frac{1}{u^{7/2}}\text{exp}[-w\left(1-w\right)u{a}^2{\overline{AB}}^2]\left[\frac{{\zeta }_A^2}{2wu{a}^2}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{{\displaystyle {\int }_0^{1}d{t}_1}}_1\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4wu{a}^2{t}_1^2}\right]$ $\left[\frac{{\zeta }_B^2}{2\left(1-w\right)u{a}^2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4\left(1-w\right)u{a}^2{t}_2^2}\right]+{\displaystyle {\int }_0^{1}du{u}^{3/2}}\text{exp}[-\frac{w\left(1-w\right)}{u}{a}^2{\overline{AB}}^2]$ $\left[\frac{u{\zeta }_A^2}{2w{a}^2}{{\displaystyle {\int }_0^{1}d{t}_1}}_1\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{{\displaystyle {\int }_0^{1}d{t}_1}}_1\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{u{\zeta }_A^2}{4w{a}^2{t}_1^2}\right]\left[\frac{u{\zeta }_B^2}{2\left(1-w\right)a^2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{u{\zeta }_B^2}{4\left(1-w\right)a^2{t}_2^2}\right]\}$ (2.9)

Integrals over w, u, t¹, and t² can be evaluated by the Gauss-Legendre quadrature. The 64 point-quadrature can give a good precision of 8 significant figures (SFs). The calculated value of the overlap integral is as OVL=0.29845563 with ζA = ζB = 1, $\overrightarrow{A}=(-\frac{\sqrt{8}}{3}, -\sqrt{\frac{8}{3}}, \frac{2}{3})$ , and $\overrightarrow{B}=(-\frac{\sqrt{8}}{3}, \sqrt{\frac{8}{3}}, \frac{2}{3})$ , which is the case of two hydrogen atoms at $\overrightarrow{A}$ and $\overrightarrow{B}$ . For the case of two carbon +5 cations, the 256-point quadrature is necessary to give the 7 SF precision. The value is as OVL=0.4570495(-6) with ζA = ζB = 6, and $\overrightarrow{A}$ and $\overrightarrow{B}$ are the same as the above.

3. Kinetic energy integral

The two-center KEI over Dirac wave functions is given by

$KEI=\frac{{\hslash }^2}{2m_e}{N}_A{N}_B{T}_{AB}$ (3.1)

Where

$T_{AB}={{\displaystyle \int }}^{\text{}}d\overrightarrow{r}{r}_A^{-{\epsilon }_A}\text{exp}(-{\zeta }_A{r}_A)(-{\nabla }^2)r_B^{-{\epsilon }_B}\text{exp}(-{\zeta }_B{r}_B)={{\displaystyle \int }}^{\text{}}d\overrightarrow{r}[\nabla r_A^{-{\epsilon }_A}\text{exp}(-{\zeta }_A{r}_A)]\cdot [\nabla r_B^{-{\epsilon }_B}\text{exp}(-{\zeta }_B{r}_B)]$ (3.2)

The Gaussian-transform for the derivative of the Dirac function can be derived in Appendix A(See Below) as given by

$\nabla r_A^{-{\epsilon }_A}\exp \left(-{\zeta }_A{r}_A\right)=\frac{-\overrightarrow{r_A}{\zeta }_A^{3+{\epsilon }_A}}{2\sqrt{\pi }\Gamma (2+{\epsilon }_A)}{\displaystyle {\int }_0^{\infty }d{S}_1{S}_1^{-3/2}{}_1}\text{exp}(-S_1{r}_A^2)\left[{\epsilon }_A{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}+{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{3+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4S_1{t}_1^2}\right]$ (3.3)

Using Eq. (3.3), we have

$T_{AB}=\frac{{\zeta }_A^{3+{\epsilon }_A}{\zeta }_B^{3+{\epsilon }_B}}{4\pi \Gamma (2+{\epsilon }_A)\Gamma (2+{\epsilon }_B)}{\displaystyle {\int }_0^{\infty }d{S}_1}{\displaystyle {\int }_0^{\infty }d{S}_2}{(S_1{S}_2)}^{-3/2}\left[{\epsilon }_A{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}+{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{3+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4S_1{t}_1^2}\right]\left[{\epsilon }_B{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}+{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{3+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4S_2{t}_2^2}\right]I_1^{(3)}$ (3.4)

Where

$I_1^{(3)}={{\displaystyle \int }}^{\text{}}d\overrightarrow{r}\overrightarrow{r_A}\cdot \overrightarrow{r_B}\text{exp}[-S_1{r}_A^2-S_2{r}_B^2]$ (3.5)

We use the Gaussian product rule, Eq. (2.5), for Eq. (3.5) and know that $\overrightarrow{r_A}\cdot \overrightarrow{r_B}=\left(\overrightarrow{r_P}+\overrightarrow{AP}\right)\cdot \left(\overrightarrow{r_P}+\overrightarrow{BP}\right)=r_P^2+\overrightarrow{r_P}\cdot \left(\overrightarrow{AP}+\overrightarrow{BP}\right)+\overrightarrow{AP}\cdot \overrightarrow{BP}$ .

Then we have

$I_1^{(3)}=exp\left[-\frac{S_1{S}_2}{S_{12}}{\overline{AB}}^2\right][I_2^{(3)}+I_3^{(3)}+I_4^{(3)}]$ (3.6)

Where

$I_2^{(3)}=4\pi {\displaystyle {\int }_0^{\infty }d{r}_P{r}_P^4}\text{exp}[-S_{12}{r}_P^2]=\frac{3{\pi }^{3/2}}{2S_{12}^{5/2}}$ (3.7)

$I_3^{(3)}=(\overrightarrow{AP}+\overrightarrow{BP})\cdot {\displaystyle {\int }_0^{\infty }d{r}_P{r}_P^2}\exp \left[-S_{12}{r}_P^2\right]{{\displaystyle \int }}^{\text{}}d\widehat{r_P}\overrightarrow{r_P}=0$ (3.8)

And

$I_4^{(3)}=4\pi \overrightarrow{AP}\cdot \overrightarrow{BP}{\displaystyle {\int }_0^{\infty }d{r}_P{r}_P^2}\text{exp}[-S_{12}{r}_P^2]=\frac{{\pi }^{3/2}{S}_1{S}_2{\overline{AB}}^2}{S_{12}^{7/2}}$ (3.9)

Let us change the integral variables as S₁₂ =z and $\frac{S_1}{S_{12}}=w$ . The Jacobian is $\frac{\partial (S_1,S_2)}{\partial (w,z)}=z$ . Thus, we have

$T_{AB}=\frac{{\zeta }_A^{3+{\epsilon }_A}{\zeta }_B^{3+{\epsilon }_B}}{\Gamma (2+{\epsilon }_A)\Gamma (2+{\epsilon }_B)}{\displaystyle {\int }_0^{1}dw}{\displaystyle {\int }_0^{\infty }dz}{[w\left(1-w\right)]}^{-3/2}\frac{\sqrt{\pi }}{8}\text{exp}[-w\left(1-w\right)z{\overline{AB}}^2]\frac{3-2w\left(1-w\right)z{\overline{AB}}^2}{z^{9/2}}$ $\left[{\epsilon }_A{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}+{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{3+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4wz{t}_1^2}\right]\left[{\epsilon }_B{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}+{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{3+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4\left(1-w\right)z{t}_2^2}\right]$ (3.10)

We separate the integral over z as given by Eq. (2.8). Then we have the final formula for the kinetic energy integral over Dirac wave functions as given by

$T_{AB}=\frac{{\zeta }_A^{3+{\epsilon }_A}{\zeta }_B^{3+{\epsilon }_B}}{\Gamma (2+{\epsilon }_A)\Gamma (2+{\epsilon }_B)}\frac{\sqrt{\pi }}{8a^7}{\displaystyle {\int }_0^{1}dw}{[w\left(1-w\right)]}^{-3/2}\{{\displaystyle {\int }_0^{1}du}\frac{3-2w\left(1-w\right)u{a}^2{\overline{AB}}^2}{u^{9/2}}\text{exp}[-w\left(1-w\right)u{a}^2{\overline{AB}}^2]\left[{\epsilon }_A{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}+{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{3+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4wu{a}^2{t}_1^2}\right]$ $\left[{\epsilon }_B{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}+{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{3+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4\left(1-w\right)u{a}^2{t}_2^2}\right]+{\displaystyle {\int }_0^{1}du{u}^{5/2}}\left[3-\frac{2w(1-w)}{u}{a}^2{\overline{AB}}^2\right]exp\left[-\frac{w(1-w)}{u}{a}^2{\overline{AB}}^2\right]$ $\left[{\epsilon }_A{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}+{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{3+{\epsilon }_A}}\right]exp\left[-\frac{u{\zeta }_A^2}{4w{a}^2{t}_1^2}\right]\left[{\epsilon }_B{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}+{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{3+{\epsilon }_B}}\right]exp\left[-\frac{u{\zeta }_B^2}{4\left(1-w\right)a^2{t}_2^2}\right]\}$ (3.11)

Integrals over w, u, t¹, and t² can be evaluated by the Gauss-Legendre quadrature. The 64 point-quadrature can give a good precision of 8 SFs. The calculated value for the kinetic energy integral is as $\frac{{\hslash }^2}{2m_e}\times 0.\text{271}0\text{6788}\left(-\text{1}\right)$ with, ${\zeta }_A={\zeta }_B=1,\overrightarrow{A}=(-\frac{\sqrt{8}}{3}, -\sqrt{\frac{8}{3}}, \frac{2}{3})$ , and $\overrightarrow{B}=(-\frac{\sqrt{8}}{3}, \sqrt{\frac{8}{3}}, \frac{2}{3})$ , which is the case of two hydrogen atoms at $\overrightarrow{A}$ and $\overrightarrow{B}$ . For the case of two carbon +5 cations, the 128-point quadrature is necessary to give the 7 SF precision. The value is as $\text{KEI}=-\frac{{\hslash }^2}{2m_e}\times 0.\text{1189531}\left(-\text{4}\right)$ with ${\zeta }_A={\zeta }_B=6$ , and $\overrightarrow{A}$ and $\overrightarrow{B}$ are the same as the above.

4. Nuclear attraction integral

4.1. Point-like nucleus

The three-center NAI over Dirac wave functions for the point-like nucleus is given by

$NA{I}_P=-Z_M{e}^2{N}_A{N}_B{U}_{AB}^{(P)}$ (4.1.1)

Where

$U_{AB}^{(P)}={{\displaystyle \int }}^{\text{}}d\overrightarrow{r_M}\frac{1}{r_M}{r}_A^{-{\epsilon }_A}{r}_B^{-{\epsilon }_B}\text{exp}[-{\zeta }_A{r}_A-{\zeta }_B{r}_B]$ (4.1.2)

where the nucleus is at $\overrightarrow{M}=(0,0,0)$ and ZMe is its nuclear charge. Using Eq. (1.10), we have

$U_{AB}^{(P)}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}}{4\pi \Gamma (1+{\epsilon }_A)\Gamma (1+{\epsilon }_B)}{\displaystyle {\int }_0^{\infty }d{S}_1}{\displaystyle {\int }_0^{\infty }d{S}_2}{(S_1{S}_2)}^{-3/2}\left[\frac{{\zeta }_A^2}{2S_1}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4S_1{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2S_2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4S_2{t}_2^2}\right]I_1^{(P)}$ (4.1.3)

where

$I_1^{(P)}={{\displaystyle \int }}^{\text{}}d\overrightarrow{r_M}\frac{1}{r_M}\text{exp}[-S_1{r}_A^2-S_2{r}_B^2]$ (4.1.4)

We use the Gaussian product rule, Eq. (2.5) and do the translation of Gaussian-type orbital (GTO) by Sack [19] given by

$\exp \left[-S_{12}{r}_P^2\right]=4\pi \text{exp}[-S_{12}{r}_M^2-S_{12}{\overline{MP}}^2]{\displaystyle \sum }_l{i}_l[2S_{12}\overline{MP}{r}_M]{\displaystyle \sum }_m{Y}_l^m(\widehat{MP})Y_l^m{(\widehat{r_M})}^{*}$ (4.1.5)

Where il(x) is the modified spherical Bessel function of the first kind and is the spherical harmonics. We use the Gaussian product rule again as given by

$exp\left[-\frac{S_1{S}_2}{S_{12}}{\overline{AB}}^2\right]\exp \left[-S_{12}{\overline{MP}}^2\right]=\text{exp}[-S_1{\overline{MA}}^2-S_2{\overline{MB}}^2]$ (4.1.6)

Then we have

$U_{AB}^{(P)}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}}{\Gamma (1+{\epsilon }_A)\Gamma (1+{\epsilon }_B)}{\displaystyle {\int }_0^{\infty }d{S}_1}{\displaystyle {\int }_0^{\infty }d{S}_2}{(S_1{S}_2)}^{-3/2}\text{exp}[-S{\overline{MA}}^2-S_2{\overline{MB}}^2]\left[\frac{{\zeta }_A^2}{2S_1}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4S_1{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2S_2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4S_2{t}_2^2}\right]I_2^{(P)}$ (4.1.7)

Where

$I_2^{(P)}={\displaystyle {\int }_0^{\infty }d{r}_M{r}_M}\exp \left[-S_{12}{r}_M^2\right]{\displaystyle {\sum }_l{i}_l}[2S_{12}\overline{MP}{r}_M]{{\displaystyle \int }}^{\text{}}d\widehat{r_M}{\displaystyle {\sum }_m{Y}_l^m}\left(\widehat{MP}\right)Y_l^m\left(\widehat{r_M}\right)*={\displaystyle {\int }_0^{\infty }d{r}_M{r}_M}\exp \left[-S_{12}{r}_M^2\right]i_0[2S_{12}\overline{MP}{r}_M]$ (4.1.8)

We use, in the above derivation, the following relation derived in a previous article [20]:

${{\displaystyle \int }}^{\text{}}d\widehat{r_M}{\displaystyle {\sum }_m{Y}_l^m}\left(\widehat{MP}\right)Y_l^m{(\widehat{r_M})}^{*}={\delta }_{l0}{\delta }_{m0}$ (4.1.9)

We know the power series of the modified spherical Bessel function as given by

$i_l\left(x\right)=\frac{x^l}{\left(2l+1\right)!!}{\displaystyle {\sum }_j\frac{{[x^2/4]}^j}{j!{(l+3/2)}_j}}$ (4.1.10)

Where (a)j = a(a+1) ⋅⋅(a+j−1) is the Pochhammer symbol. Using Eq. (4.1.10), we have

$I_2^{(P)}=\frac{1}{2}{\displaystyle {\sum }_j\frac{{[S_{12}^2{\overline{MP}}^2]}^j}{j!{(\frac{3}{2})}_j}\frac{\Gamma (j+1)}{S_{12}^{j+1}}=\frac{1}{2S_{12}}{}_{1}F{}_1[1;\frac{3}{2};S_{12}{\overline{MP}}^2]}=\frac{1}{2S_{12}}\exp \left[S_{12}{\overline{MP}}^2\right]{}_{1}F{}_1[\frac{1}{2};\frac{3}{2};-S_{12}{\overline{MP}}^2]=\frac{1}{2S_{12}}\exp \left[S_{12}{\overline{MP}}^2\right]F_0[S_{12}{\overline{MP}}^2]$ (4.1.11)

where $F_m\left(x\right)={\displaystyle {\int }_0^{1}dt{t}^{2m}}\text{exp}(-x{t}^2)$ is the molecular incomplete gamma function. Substituting Eq. (4.1.11) into Eq. (4.1.7), we have

$U_{AB}^{(P)}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}}{\Gamma (1+{\epsilon }_A)\Gamma (1+{\epsilon }_B)}{\displaystyle {\int }_0^{\infty }d{S}_1}{\displaystyle {\int }_0^{\infty }d{S}_2}\frac{F_0[S_{12}{\overline{MP}}^2]}{2S_{12}{(S_1{S}_2)}^{3/2}}\text{exp}[-\frac{S_1{S}_2}{S_{12}}{\overline{AB}}^2]\left[\frac{{\zeta }_A^2}{2S_1}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4S_1{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2S_2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4S_2{t}_2^2}\right]$ (4.1.12)

In the above derivation, we use the relation given by $\exp \left[-S_1{\overline{MA}}^2-S_2{\overline{MB}}^2+S_{12}{\overline{MP}}^2\right]=\text{exp}\left[-\frac{S_1{S}_2}{S_{12}}{\overline{AB}}^2\right]$ . Let us change the integral variables as S₁₂ =z and $\frac{S_1}{S_{12}}=w$ . The Jacobian is $\frac{\partial (S_1,S_2)}{\partial (w,z)}=z$ . Thus, we have

$U_{AB}^{(P)}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}}{\Gamma (1+{\epsilon }_A)\Gamma (1+{\epsilon }_B)}{\displaystyle {\int }_0^{1}dw}{[w\left(1-w\right)]}^{-3/2}{\displaystyle {\int }_0^{\infty }dz}\frac{1}{2z^3}\text{exp}[-w\left(1-w\right)z{\overline{AB}}^2]F_0(z{x}_0)$ $\left[\frac{{\zeta }_A^2}{2wz}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4wz{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2\left(1-w\right)z}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4\left(1-w\right)z{t}_2^2}\right]$ (4.1.13)

where $x_0=w^2{\overline{MA}}^2+{(1-w)}^2{\overline{MB}}^2+2w(1-w)\overrightarrow{MA}\cdot \overrightarrow{MB}$ . We separate the integral over z as given by Eq. (2.8). Then we have the final formula of the nuclear attraction integral over Dirac wave functions for the point-like nucleus as given by

$U_{AB}^{(P)}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}}{\Gamma (1+{\epsilon }_A)\Gamma (1+{\epsilon }_B)}\frac{1}{2a^4}{\displaystyle {\int }_0^{1}dw}{[w\left(1-w\right)]}^{-3/2}\{{\displaystyle {\int }_0^{1}du}\frac{1}{u^3}\text{exp}[-w\left(1-w\right)u{a}^2{\overline{AB}}^2]F_0(u{a}^2{x}_0)$ $\left[\frac{{\zeta }_A^2}{2wu{a}^2}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4wu{a}^2{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2\left(1-w\right)u{a}^2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4\left(1-w\right)u{a}^2{t}_2^2}\right]$ $+{\displaystyle {\int }_0^{1}duu} exp\left[-\frac{w(1-w)}{u}{a}^2{\overline{AB}}^2\right]F_0\left(\frac{a^2{x}_0}{u}\right)\left[\frac{u{\zeta }_A^2}{2w{a}^2}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{u{\zeta }_A^2}{4w{a}^2{t}_1^2}\right]\left[\frac{u{\zeta }_B^2}{2\left(1-w\right)a^2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{u{\zeta }_B^2}{4\left(1-w\right)a^2{t}_2^2}\right]\}$ (4.1.14)

Integrals over w, u, t¹, and t² can be evaluated by the Gauss-Legendre quadrature. The 64 point-quadrature can give a good precision of 8 SFs. The calculated value for the nuclear attraction integral is as NAIP = -ZMe2 × 0.16764413 with ζA = ζB = 1, ZM=1, $\overrightarrow{M}=\left(0, 0, 0\right),$ $\overrightarrow{A}=(-\frac{\sqrt{8}}{3}, -\sqrt{\frac{8}{3}}, \frac{2}{3})$ , and $\overrightarrow{B}=(-\frac{\sqrt{8}}{3}, \sqrt{\frac{8}{3}}, \frac{2}{3})$ , which is the case of three hydrogen atoms at $\overrightarrow{M}$ , $\overrightarrow{A}$ and $\overrightarrow{B}$ . For the case of three carbon +5 cations, the 128-point quadrature is necessary to give the 7 SF precision. The value is as NAIp = -e2 × 0.2062558(-5) with ζA = ζB = 6, ZM = 6, and $\overrightarrow{M}$ , $\overrightarrow{A}$ and $\overrightarrow{B}$ are the same as the above.

4.2. Gauss-type charge density distribution model

Some experiment shows the real nucleus is not a point-like one but a finite one [21]. Among the finite nucleus model, the Gauss-type charge density distribution (GCDD) model is frequently used [22-24]. For the GCDD model, the nuclear attraction potential is given by [21],

$V=-Z_M{e}^2\frac{\text{erf}\left(\sqrt{\xi }{r}_M\right)}{r_M}=-Z_M{e}^2\frac{2}{\sqrt{\pi }{r}_0}{F}_0\left(\frac{r_M^2}{r_0^2}\right)$ (4.2.1)

Where erf (x) is the error function and each of the size parameters ξ and r0 is relative to each other as $r_0=\sqrt{1/\xi }$ . The value of r0 is very small as r0 = 0.2169394461(−4) for the hydrogen atom. The three-center NAI over Dirac wave functions for the GCDD model is given by

$NA{I}_G=-Z_M{e}^2{N}_A{N}_B{U}_{AB}^{(G)}$ (4.2.2)

where

$U_{AB}^{(G)}=\frac{2}{\sqrt{\pi }{r}_0}{{\displaystyle \int }}^{\text{}}d\overrightarrow{r_M}{F}_0\left(\frac{r_M^2}{r_0^2}\right)r_A^{-{\epsilon }_A}{r}_B^{-{\epsilon }_B}\text{exp}[-{\zeta }_A{r}_A-{\zeta }_B{r}_B]$ (4.2.3)

Using Eq. (1.10), Gaussian product rule, Eq. (2.5), Sack’s translation of GTO, Eq. (4.1.5), Eq. (4.1.6), and Eq. (4.1.9), we have

$U_{AB}^{(G)}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}}{\Gamma (1+{\epsilon }_A)\Gamma (1+{\epsilon }_B)}{\displaystyle {\int }_0^{\infty }d{S}_1}{\displaystyle {\int }_0^{\infty }d{S}_2}{(S_1{S}_2)}^{-3/2}\text{exp}[-S_1{\overline{MA}}^2-S_2{\overline{MB}}^2]\left[\frac{{\zeta }_A^2}{2S_1}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4S_1{t}_1^2}\right]$ $\left[\frac{{\zeta }_B^2}{2S_2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4S_2{t}_2^2}\right]I_1^{(G)}$ (4.2.4)

where

$I_1^{(G)}=\frac{2}{\sqrt{\pi }{r}_0}{{\displaystyle {\int }_0^{\infty }d{r}_M{r}_M^{2}F}}_0\left(\frac{r_M^2}{r_0^2}\right)\exp \left[-S_{12}{r}_M^2\right]i_0[2S_{12}\overline{MP}{r}_M]=I_1^{(G)in}+I_1^{\left(G\right)out}$ (4.2.5)

in which

$I_1^{(G)in}=\frac{2}{\sqrt{\pi }{r}_0}{\displaystyle {\int }_0^{R_0}d{r}_M{r}_M^2}{F}_0\left(\frac{r_M^2}{r_0^2}\right)\exp \left[-S_{12}{r}_M^2\right]i_0[2S_{12}\overline{MP}{r}_M]$ (4.2.6)

and

$I_1^{(G)out}=\frac{2}{\sqrt{\pi }{r}_0}{\displaystyle {\int }_{R_0}^{\infty }d{r}_M{r}_M^2}{F}_0\left(\frac{r_M^2}{r_0^2}\right)\exp \left[-S_{12}{r}_M^2\right]i_0[2S_{12}\overline{MP}{r}_M]$ (4.2.7)

With R0 = bro. Here we choose b from the asymptotic expansion of the molecular incomplete gamma function as given by

$F_m\left(\frac{r_M^2}{r_0^2}\right)=\frac{\Gamma (m+\frac{1}{2})}{2}{\left(\frac{r_0}{r_M}\right)}^{2m+1} (r_M\to \infty )$ (4.2.8)

In a previous article [25], the author shows that $F_0\left(\frac{r_M^2}{r_0^2}\right)$ becomes its asymptotic value for $r_M^2>36r_0^2$ within 15 significant-figure precision, $F_1\left(\frac{r_M^2}{r_0^2}\right)$ does for $r_M^2>40r_0^2$ , and $F_2\left(\frac{r_M^2}{r_0^2}\right)$ does for $r_M^2>43r_0^2$ . We choose b = 7(b2 = 49), here. At rM > R0 =7r0, we know that $\frac{2}{\sqrt{\pi }{r}_0}{F}_0\left(\frac{r_M^2}{r_0^2}\right)=\frac{1}{r_M}$ and $\frac{4}{\sqrt{\pi }{r}_0^3}{F}_1\left(\frac{r_M^2}{r_0^2}\right)=\frac{1}{r_M^3}$ . Thus, for r > R0, the scalar potential is equal to that for the point-like nucleus as given by $-Z_M{e}^2\frac{2}{\sqrt{\pi }}{F}_0\left(\frac{r^2}{r_0^2}\right)=\frac{-Z_M{e}^2}{r}$ and the vector potential is also given by $\frac{4Z_{M}e}{\sqrt{\pi }{r}_0^3}{F}_1\left(\frac{r_M^2}{r_0^2}\right)\overrightarrow{\mu }\times \overrightarrow{r_M}=\frac{Z_{M}e}{r_M^3}\overrightarrow{\mu }\times \overrightarrow{r_M}$ . We can recognize rM > R0 is clearly the outer part of the finite nucleus of the GCDD model. We know the power series of Fm(x) given by

$F_m\left(x\right)=\frac{1}{2m+1}{}_{1}F{}_1\left(m+\frac{1}{2};m+\frac{3}{2};x\right)=\frac{1}{2m+1}{\displaystyle {\sum }_k\frac{x^k{(m+\frac{1}{2})}_k}{k!{(m+\frac{3}{2})}_k}}$ (4.2.9)

Using Eq. (4.1.10) and (4.2.9), we have

$I_1^{\left(G\right)in}=\frac{2}{\sqrt{\pi }{r}_0}{\displaystyle \sum }_j\frac{{[S_{12}^2{\overrightarrow{MP}}^2]}^j}{j!{(3/2)}_j}{\displaystyle {\sum }_k\frac{{(-1/r_0^2)}^k{(1/2)}_k}{k!{(3/2)}_k}{I}_2^{\left(G\right)in}}$ (4.2.10)

where

$=\frac{1}{2}{R}_0^{2j+2k+3}\frac{\Gamma (j+k+3/2)}{\Gamma (j+k+5/2)}{}_{1}F{}_1[j+k+\frac{3}{2};j+k+\frac{5}{2}; -S_{12}{R}_0^2]$ (4.2.11)

in which γ(x) is the incomplete gamma function of the first kind and 1F1(α; γ; x) is the confluent hypergeometric function. Substituting Eq. (4.2.11) into Eq. (4.2.10), we have

$I_1^{\left(G\right)in}=\frac{b^3{r}_0^2}{\sqrt{\pi }}{\displaystyle {\sum }_j\frac{{[S_{12}^2{\overline{MP}}^2{R}_0^2]}^j}{j!{(3/2)}_j}{\displaystyle {\sum }_k\frac{{(-b^2)}^k{(1/2)}_k}{k!{(3/2)}_k}}}\frac{\Gamma (j+k+3/2)}{\Gamma (j+k+5/2)}{}_{1}F{}_1[j+k+\frac{3}{2};j+k+\frac{5}{2}; -S_{12}{R}_0^2]$ $=\frac{b^3{r}_0^2}{\sqrt{\pi }}{\displaystyle {\sum }_k\frac{{\left(-b^2\right)}^k{(1/2)}_k}{k!{(3/2)}_k}\frac{\Gamma (k+3/2)}{\Gamma (k+5/2)}+O(R_0^4)}=\frac{b^3{r}_0^2}{\sqrt{\pi }}\frac{\Gamma \left(3/2\right)}{\Gamma \left(5/2\right)}{}_{1}F{}_1\left(\frac{1}{2};\frac{5}{2}; -b^2\right)+O\left(R_0^4\right)$ $=\frac{b^3{r}_0^2}{\sqrt{\pi }}\frac{\Gamma \left(3/2\right)}{\Gamma \left(5/2\right)}\frac{\Gamma (5/2)}{b}\left[1-\frac{1}{2b^2}\right]+O(R_0^4)=\frac{1}{2}{R}_0^2-\frac{r_0^2}{4}+O(R_0^4)$ (4.2.12)

Next, we evaluate Eq. (4.2.7). We use the asymptotic expansion Eq. (4.2.8) for Eq. (4.2.7) and have

$I_1^{\left(G\right)out}={\displaystyle {\int }_{R_0}^{\infty }d{r}_M r_M}\exp \left[-S_{12}{r}_M^2\right]i_0[2S_{12}\overline{MP}{r}_M]$ (4.2.13)

Using Eq. (4.1.10), we have

$I_1^{\left(G\right)out}={\displaystyle {\sum }_j\frac{{[S_{12}^2{\overline{MP}}^2]}^j}{j!{(3/2)}_j}}{\displaystyle {\int }_{R_0}^{\infty }d{r}_M r_M^{2j+1}}\text{exp}[-S_{12}{r}_M^2]=\frac{1}{2}{\displaystyle {\sum }_j\frac{{[S_{12}^2{\overline{MP}}^2]}^j}{j!{(3/2)}_j}}{\left(\frac{1}{S_{12}}\right)}^{j+1}\Gamma [j+1;S_{12}{R}_0^2]$ $=\frac{1}{2S_{12}}{\displaystyle {\sum }_j\frac{{[S_{12}{\overline{MP}}^2]}^j}{j!{(3/2)}_j}}[\Gamma \left(j+1\right)-\frac{{[S_{12}{R}_0^2]}^{j+1}}{j+1}+O(R_0^4)]=\frac{1}{2S_{12}}{}_{1}F{}_1\left(1;\frac{3}{2}; S_{12}{\stackrel{-}{MP}}^2\right)-\frac{1}{2}{R}_0^2+O(R_0^4)$ $=\frac{1}{2S_{12}}\exp \left[S_{12}{\overline{MP}}^2\right]F_0\left(S_{12}{\overline{MP}}^2\right)-\frac{1}{2}{R}_0^2+O(R_0^4)$ (4.2.14)

where is the incomplete gamma function of the second kind. It is easy to derive the following relation:

$\Gamma \left(\alpha ;x\right)=\Gamma \left(\alpha \right)-\frac{x^{\alpha }}{\alpha }+\frac{x^{\alpha +1}}{\alpha +1}+\cdots (x\ll 1)$ (4.2.15)

We use Eq. (4.2.15) in deriving Eq. (4.2.14). Substituting Eq. (4.2.12) and (4.2.14) into Eq. (4.2.5), we have

$I_1^{(G)}=\frac{1}{2S_{12}}\exp \left[S_{12}{\overline{MP}}^2\right]F_0\left(S_{12}{\overline{MP}}^2\right)-\frac{r_0^2}{4}+O(R_0^4)$ (4.2.16)

Substituting Eq. (4.2.16) into Eq. (4.2.4), we have

$U_{AB}^{(G)}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}}{\Gamma (1+{\epsilon }_A)\Gamma (1+{\epsilon }_B)}{\displaystyle {\int }_0^{\infty }d{S}_1}{\displaystyle {\int }_0^{\infty }d{S}_2}{(S_1{S}_2)}^{-3/2}\text{exp}[-S_1{\overline{MA}}^2-S_2{\overline{MB}}^2][\frac{1}{2S_{12}}\exp \left[S_{12}{\overline{MP}}^2\right]F_0\left(S_{12}{\overline{MP}}^2\right)-\frac{r_0^2}{4}]$ $\left[\frac{{\zeta }_A^2}{2S_1}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4S_1{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2S_2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4S_2{t}_2^2}\right]+O(R_0^4)$ (4.2.17)

Let us change the integral variables as S₁₂ =z and $\frac{S_1}{S_{12}}=w$ . The Jacobian is $\frac{\partial (S_1,S_2)}{\partial (w,z)}=z$ . Thus, we have

$U_{AB}^{(G)}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}}{\Gamma (1+{\epsilon }_A)\Gamma (1+{\epsilon }_B)}{\displaystyle {\int }_0^{1}dw}{[w\left(1-w\right)]}^{-3/2}{\displaystyle {\int }_0^{\infty }dz}\frac{1}{z^3}\text{exp}[-wz{\overline{MA}}^2-\left(1-w\right)z{\overline{MB}}^2][\frac{1}{2z}\exp \left[z{x}_0\right]F_0\left(z{x}_0\right)-\frac{r_0^2}{4}]$ $\left[\frac{{\zeta }_A^2}{2wz}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4wz{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2\left(1-w\right)z}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4\left(1-w\right)z{t}_2^2}\right]+O(R_0^4)$ (4.2.18)

We separate the integral over z as given by Eq. (2.8). Then we have the final formula of the nuclear attraction integral over Dirac wave functions for the GCDD model as given by

$U_{AB}^{(G)}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}}{\Gamma (1+{\epsilon }_A)\Gamma (1+{\epsilon }_B)}\frac{1}{2a^4}{\displaystyle {\int }_0^{1}dw}{[w\left(1-w\right)]}^{-3/2}\{{\displaystyle {\int }_0^{1}du}\frac{1}{u^3}\text{exp}[-wu{a}^2{\overline{MA}}^2-\left(1-w\right)u{a}^2{\overline{MB}}^2][\frac{1}{u{a}^2}\exp \left[u{a}^2{x}_0\right]F_0\left(u{a}^2{x}_0\right)-\frac{r_0^2}{2}]$ $\left[\frac{{\zeta }_A^2}{2wu{a}^2}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4wu{a}^2{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2\left(1-w\right)u{a}^2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4\left(1-w\right)u{a}^2{t}_2^2}\right]$ $+{\displaystyle {\int }_0^{1}duu} \text{exp}[-\frac{w}{u}{a}^2{\overline{MA}}^2-\frac{1-w}{u}{a}^2{\overline{MB}}^2]\text{ [}\frac{u}{a^2}\exp \left[\frac{a^2}{u}{x}_0\right]F_0\left(\frac{a^2}{u}{x}_0\right)-\frac{r_0^2}{2}]\left[\frac{u{\zeta }_A^2}{2w{a}^2}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{u{\zeta }_A^2}{4w{a}^2{t}_1^2}\right]$ $\left[\frac{u{\zeta }_B^2}{2\left(1-w\right)a^2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{u{\zeta }_B^2}{4\left(1-w\right)a^2{t}_2^2}\right]\}+O(R_0^4)$ (4.2.19)

The error term, $O(R_0^4)$ , is in the order of $R_0^4$ , of which value is as $R_0^4=0.53179747(-15)$ . Integrals over w, u, t¹, and t² can be evaluated by the Gauss-Legendre quadrature. The 64 point-quadrature can give a good precision of 8 SFs. The calculated value for the nuclear attraction integral is as NAIIG = -ZMe2 × 0.16764413 with ζA = ζB = 1, $\overrightarrow{M}=\left(0, 0, 0\right),$ , $\overrightarrow{A}=(-\frac{\sqrt{8}}{3}, -\sqrt{\frac{8}{3}}, \frac{2}{3})$ and $\overrightarrow{B}=(-\frac{\sqrt{8}}{3}, \sqrt{\frac{8}{3}}, \frac{2}{3})$ , which is the case of three hydrogen atoms at $\overrightarrow{M}$ , $\overrightarrow{A}$ and $\overrightarrow{B}$ . Thus, the value is the same as that of the pint-like nucleus within the 8 SF precision. The contribution from the term $-\frac{r_0^2}{2}$ in Eq. (4.2.19) is as ZMe2 × 0.9110420(-9). Using the 128 point-quadrature, we find the small difference as follows: The value is as NAIP = -ZMe2 × 0.167644127 for the point-like nucleus and as NAIG = -ZMe2 × 0.167644126 for the GCDD model. The difference is very small because that $r_0^2=0.\text{47}0\text{6272323}\left(-\text{9}\right)$ (for the case of the hydrogen atom) is very small. For the case of three carbon +5 cations, the 128-point quadrature is necessary to give the 7 SF precision. The value is as NAIG = -e2 × 0.2062558(-5) with ζA = ζB = 6, ZM = 6, and $\overrightarrow{M}$ , $\overrightarrow{A}$ and $\overrightarrow{B}$ are the same as the above. Thus, the value is also the same as that of the pint-like nucleus within the 7 SF precision, because that $r_0^2=0.146891404(-8)$ is also very small for the case of the carbon +5 cation. The contribution from the term $-\frac{r_0^2}{2}$ in Eq. (4.2.19) is e2 × 0.65044(-16). So, the difference between the value for the point-like nucleus and that for the GCDD model can not be found until that as in 11 SF precision. One can up the precision by upping the point number of the quadrature if desired.

5. Electron-repulsion integral

The four-center ERI over Dirac wave functions is given by

$ERI=e^2{N}_A{N}_B{N}_C{N}_D{V}_{ABCD}$ (5.1)

where

$V_{ABCD}={{\displaystyle \int }}^{\text{}}d\overrightarrow{r_1} d\overrightarrow{r_2}\frac{1}{r_{12}}{r}_{1A}^{-{\epsilon }_A}{r}_{1B}^{-{\epsilon }_B}{r}_{2C}^{-{\epsilon }_C}{r}_{2D}^{-{\epsilon }_D}\text{exp}[-{\zeta }_A{r}_{1A}-{\zeta }_B{r}_{1B}-{\zeta }_C{r}_{2C}-{\zeta }_D{r}_{2D}]$ (5.2)

Using Eq. (1.10). we have

$V_{ABCD}=\frac{1}{{(4\pi )}^2}\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}{\zeta }_C^{1+{\epsilon }_C}{\zeta }_D^{1+{\epsilon }_D}}{\Gamma \left(1+{\epsilon }_A\right)\Gamma \left(1+{\epsilon }_B\right)\Gamma \left(1+{\epsilon }_C\right)\Gamma \left(1+{\epsilon }_D\right)}{\displaystyle {\int }_0^{\infty }d{S}_1}{\displaystyle {\int }_0^{\infty }d{S}_2}{\displaystyle {\int }_0^{\infty }d{S}_3}{\displaystyle {\int }_0^{\infty }d{S}_4}{(S_1{S}_2{S}_3{S}_4)}^{-3/2}$ $\left[\frac{{\zeta }_A^2}{2S_1}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4S_1{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2S_2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4S_2{t}_2^2}\right]$ $\left[\frac{{\zeta }_C^2}{2S_3}{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{4+{\epsilon }_C}}-{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{2+{\epsilon }_C}}\right]exp\left[-\frac{{\zeta }_C^2}{4S_3{t}_3^2}\right]\left[\frac{{\zeta }_D^2}{2S_4}{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{4+{\epsilon }_D}}-{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{2+{\epsilon }_D}}\right]exp\left[-\frac{{\zeta }_D^2}{4S_4{t}_4^2}\right]I_1^{(5)}$ (5.3)

where

$I_1^{(5)}={{\displaystyle \int }}^{\text{}}d\overrightarrow{r_1} d\overrightarrow{r_2} \frac{1}{r_{12}}\text{exp}[-S_1{r}_{1A}^2-S_2{r}_{1B}^2-S_3{r}_{2C}^2-S_4{r}_{2D}^2]$ (5.4)

We use the Gaussian product rule, Eq. (2.5) and have

$I_1^{\left(5\right)}=exp\left[-\frac{S_1{S}_2}{S_{12}}{\overline{AB}}^2-\frac{S_3{S}_4}{S_{34}}C{D}^2\right]I_2^{(5)}$ (5.5)

where

$I_2^{(5)}={{\displaystyle \int }}^{\text{}}d\overrightarrow{r_1} d\overrightarrow{r_2} \frac{1}{r_{12}}\text{exp}[-S_{12}{r}_{1P}^2-S_{34}{r}_{2Q}^2]$ (5.6)

in which S34 = S3 + S4 and $\overrightarrow{Q}=\frac{S_3}{S_{34}}\overrightarrow{C}+\frac{S_4}{S_{34}}\overrightarrow{D}$ . We know the Fourier transform of 1/r12 given by

$\frac{1}{r_{12}}=\frac{1}{2{\pi }^2}{{\displaystyle \int }}^{\text{}}\overrightarrow{K}\frac{1}{K^2}\text{exp}[i\overrightarrow{K}\cdot \left(\overrightarrow{r_{2Q}}-\overrightarrow{r_{1P}}+\overrightarrow{PQ}\right)]$ (5.7)

We use the partial wave expansion of the plane wave as given by

$\exp \left[i\overrightarrow{K}\cdot \overrightarrow{r_{2Q}}\right]=4\pi {\displaystyle {\sum }_{{}_2}{i}^l}{j}_{l_2}(K{r}_{2Q}){\displaystyle {\sum }_{m_2}{Y}_{l_2}^{m_2}}(\widehat{K})Y_{l_2}^{m_2}{(\widehat{r_{2Q}})}^{*}$ (5.8)

$\exp \left[-i\overrightarrow{K}\cdot \overrightarrow{r_{1P}}\right]=4\pi {\displaystyle {\sum }_{l_1}{(-i)}^{l_1}}{j}_{l_1}(K{r}_{1P}){\displaystyle {\sum }_{m_1}{Y}_{l_1}^{m_1}}{(\widehat{K})}^{*}{Y}_{l_1}^{m_1}(\widehat{r_{1P}})$ (5.9)

and

$\exp \left[i\overrightarrow{K}\cdot \overrightarrow{PQ}\right]=4\pi {\displaystyle {\sum }_l{i}^l}{j}_l(K\overrightarrow{PQ}){\displaystyle {\sum }_m{Y}_l^m}(\widehat{K})Y_l^m{(\widehat{PQ})}^{*}$ (5.10)

Where jl (x) is the spherical Bessel function.

Using Eq. (5.7)-(5.10), we have

$I_2^{(5)}=\frac{{\left(4\pi \right)}^3}{2{\pi }^2}{\displaystyle {\sum }_{l, l_1,l_2}{i}^{l+l_2}}{(-i)}^{l_1}{{\displaystyle \int }}^{\text{}}d\overrightarrow{K}\frac{1}{K^2}{j}_l(K\overline{PQ}){\displaystyle {\sum }_m{Y}_l^m}(\widehat{K})Y_l^m{(\widehat{PQ})}^{*}{\displaystyle {\int }_0^{\infty }d{r}_{1P}{r}_{1P}^2}\text{exp}[-S_{12}{r}_{1P}^2]j_{l_1}(K{r}_{1P}){{\displaystyle \int }}^{\text{}}d\widehat{r_{1P}}{\displaystyle {\sum }_{m_1}{Y}_{l_1}^{m_1}}(\widehat{K})Y_{l_1}^{m_1}{(\widehat{r_{1P}})}^{*}$ ${\displaystyle {\int }_0^{\infty }d{r}_{2Q} r_{2Q}^2}\text{exp}[-S_{34}{r}_{2Q}^2]j_{l_2}(K{r}_{2Q}){{\displaystyle \int }}^{\text{}}d\widehat{r_{2Q}}{\displaystyle {\sum }_{m_2}{Y}_{l_2}^{m_2}}{(\widehat{K})}^{*}{Y}_{l_2}^{m_2}(\widehat{r_{2Q}})$ (5.11)

Using Eq. (4.1.9) for the angular parts in Eq. (5.11), we have

$I_2^{(5)}=32\pi {\displaystyle {\sum }_l{i}^l}{\displaystyle {\int }_0^{\infty }dK{j}_l}(K\overline{PQ}){{\displaystyle \int }}^{\text{}}d\widehat{K}{\displaystyle {\sum }_m{Y}_l^m}\left(\widehat{K}\right)Y_l^m{(\widehat{PQ})}^{*}{\displaystyle {\int }_0^{\infty }d{r}_{1P}{r}_{1P}^2}\exp \left[-S_{12}{r}_{1P}^2\right] j_0(K{r}_{1P}){\displaystyle {\int }_0^{\infty }d{r}_{2Q} r_{2Q}^2}\exp \left[-S_{34}{r}_{2Q}^2\right] j_0(K{r}_{2Q})$ (5.12)

In the textbook by Watson [26], we have

${\displaystyle {\int }_0^{\infty }dt{J}_{\nu }\left(at\right)}\exp \left(-p^2{t}^2\right) t^{\mu -1}=\frac{\Gamma (\frac{1}{2}\nu +\frac{1}{2}\mu ){(\frac{1}{2}a/p)}^{\nu }}{2p^{\mu }\Gamma \left(\nu +1\right)}{}_{1}F{}_1(\frac{1}{2}\nu +\frac{1}{2}\mu ; \nu +1; -\frac{a^2}{4p^2})$ (5.13)

Where jv (x) is the Bessel function. We know that $j_n\left(x\right)=\sqrt{\frac{\pi }{2x}}{J}_{n+\frac{1}{2}}(x)$ . Using Eq. (5.13) (with v = n+ 1/2 and µ = m + 1/2) and Eq. (4.1.9), we have

$I_2^{(5)}=32\pi \frac{\sqrt{\pi }}{4S_{12}^{3/2}}\frac{\sqrt{\pi }}{4S_{34}^{3/2}}{\displaystyle {\int }_0^{\infty }dK{j}_0}(\text{K}\stackrel{-}{PQ})\text{exp}[-\frac{K^2}{4S_{12}}-\frac{K^2}{4S_{34}}]=\frac{2{\pi }^2}{{(S_{12}{S}_{34})}^{3/2}}\frac{\Gamma \left(1/2\right)\sqrt{\pi }}{4\Gamma \left(3/2\right)}{\frac{2\sqrt{S_{12}{S}_{34}}}{\sqrt{S_{12}+S_{34}}}}_1{}_{1}F{}_1(\frac{1}{2}; \frac{3}{2}; -\frac{S_{12}{S}_{34}}{S_{12}+S_{34}}{\overline{PQ}}^2)=\frac{2{\pi }^{5/2}}{S_{12}{S}_{34}\sqrt{S_{12}+S_{34}}} F_0\left[\frac{S_{12}{S}_{34}}{S_{12}+S_{34}}{\overline{PQ}}^2\right]$ (5.14)

Substituting Eq. (5.14) into Eq. (5.5) and doing Eq. (5.5) into Eq. (5.3), we have

$V_{ABCD}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}{\zeta }_C^{1+{\epsilon }_C}{\zeta }_D^{1+{\epsilon }_D}}{\Gamma \left(1+{\epsilon }_A\right)\Gamma \left(1+{\epsilon }_B\right)\Gamma \left(1+{\epsilon }_C\right)\Gamma \left(1+{\epsilon }_D\right)}\frac{\sqrt{\pi }}{8}{\displaystyle {\int }_0^{\infty }d{S}_1}{\displaystyle {\int }_0^{\infty }d{S}_2}{\displaystyle {\int }_0^{\infty }d{S}_3}{\displaystyle {\int }_0^{\infty }d{S}_4}\frac{F_0\left(z_0\right)\text{exp}[-(S_1{S}_2/S_{12}){\overline{AB}}^2-(S_3{S}_4/S_{34}){\overline{CD}}^2]}{{(S_1{S}_2{S}_3{S}_4)}^{3/2}{S}_{12}{S}_{34}\sqrt{S_{12}+S_{34}}}$ $\left[\frac{{\zeta }_A^2}{2S_1}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4S_1{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2S_2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4S_2{t}_2^2}\right]$ $\left[\frac{{\zeta }_C^2}{2S_3}{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{4+{\epsilon }_C}}-\underset{0}{\overset{1}{{\displaystyle \int }}}d{t}_3\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{2+{\epsilon }_C}}\right]exp\left[-\frac{{\zeta }_C^2}{4S_3{t}_3^2}\right]\left[\frac{{\zeta }_D^2}{2S_4}{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{4+{\epsilon }_D}}-{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{2+{\epsilon }_D}}\right]exp\left[-\frac{{\zeta }_D^2}{4S_4{t}_4^2}\right]$ (5.15)

where $z_0=\frac{S_{12}{S}_{34}}{S_{12}+S_{34}}{\overline{PQ}}^2$ . Let us change the integral variables as S12 + S34 = z, $\frac{S_{12}}{S_{12}+S_{34}}=w$ , $\frac{S_1}{S_{12}}=u$ , and $\frac{S_3}{S_{34}}=v$ .

The Jacobian is $\frac{\partial (S_1,S_2,S_3,S_4)}{\partial (z, w, u, v)}=w(1-w)z^3$ . Then we have

$V_{ABCD}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}{\zeta }_C^{1+{\epsilon }_C}{\zeta }_D^{1+{\epsilon }_D}}{\Gamma \left(1+{\epsilon }_A\right)\Gamma \left(1+{\epsilon }_B\right)\Gamma \left(1+{\epsilon }_C\right)\Gamma \left(1+{\epsilon }_D\right)}\frac{\sqrt{\pi }}{8}{\displaystyle {\int }_0^{1}dw}{\displaystyle {\int }_0^{1}du}{\displaystyle {\int }_0^{1}dv}{\displaystyle {\int }_0^{\infty }dz}\frac{F_0(w\left(1-w\right)z{\overline{PQ}}^2)\text{exp}[-u\left(1-u\right)wz{\overline{AB}}^2-v\left(1-v\right)\left(1-w\right)z{\overline{CD}}^2]}{{[u\left(1-u\right)v(1-v)]}^{3/2}{[w(1-w)]}^3{z}^{11/2}}$ $\left[\frac{{\zeta }_A^2}{2uwz}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]\text{exp}\left[-\frac{{\zeta }_A^2}{4uwz{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2\left(1-u\right)wz}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]\text{exp}\left[-\frac{{\zeta }_B^2}{4\left(1-u\right)wz{t}_2^2}\right]$ $\left[\frac{{\zeta }_C^2}{2v\left(1-w\right)z}{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{4+{\epsilon }_C}}-{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{2+{\epsilon }_C}}\right]\text{exp}\left[-\frac{{\zeta }_C^2}{4v\left(1-w\right)z{t}_3^2}\right]\left[\frac{{\zeta }_D^2}{2\left(1-v\right)\left(1-w\right)z}{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{4+{\epsilon }_D}}-{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{2+{\epsilon }_D}}\right]\text{exp}\left[-\frac{{\zeta }_D^2}{4\left(1-v\right)\left(1-w\right)z{t}_4^2}\right]$ (5.16)

We separate the integral over z similarly to Eq. (2.8) as given by

${\displaystyle {\int }_0^{\infty }dz}=a^2{\displaystyle {\int }_0^{1}ds}+a^2{\displaystyle {\int }_0^{1}ds}\frac{1}{s^2}$ (5.17)

Then we have the final formula of the electron-repulsion integral over Dirac wave functions as given by

$V_{ABCD}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}{\zeta }_C^{1+{\epsilon }_C}{\zeta }_D^{1+{\epsilon }_D}}{\Gamma \left(1+{\epsilon }_A\right)\Gamma \left(1+{\epsilon }_B\right)\Gamma \left(1+{\epsilon }_C\right)\Gamma \left(1+{\epsilon }_D\right)}\frac{\sqrt{\pi }}{8a^9}{\displaystyle {\int }_0^{1}dw}{\displaystyle {\int }_0^{1}du}{\displaystyle {\int }_0^{1}dv}{[u\left(1-u\right)v(1-v)]}^{-3/2}{[w(1-w)]}^{-3}$ $\{{\displaystyle {\int }_0^{1}ds}\frac{1}{s^{11/2}}{F}_0\left[z_1\right]\text{exp}[-u\left(1-u\right)ws{a}^2{\overline{AB}}^2-v\left(1-v\right)(1-w)s{a}^2{\overline{CD}}^2]\left[\frac{{\zeta }_A^2}{2uws{a}^2}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]\text{exp}\left[-\frac{{\zeta }_A^2}{4uws{a}^2{t}_1^2}\right]$ $\left[\frac{{\zeta }_B^2}{2\left(1-u\right)ws{a}^2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]\text{exp}\left[-\frac{{\zeta }_B^2}{4\left(1-u\right)ws{a}^2{t}_2^2}\right]\left[\frac{{\zeta }_C^2}{2v\left(1-w\right)s{a}^2}{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{4+{\epsilon }_C}}-{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{2+{\epsilon }_C}}\right]\text{exp}\left[-\frac{{\zeta }_C^2}{4v\left(1-w\right)s{a}^2{t}_3^2}\right]$ $\left[\frac{{\zeta }_D^2}{2\left(1-v\right)\left(1-w\right)s{a}^2}{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{4+{\epsilon }_D}}-{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{2+{\epsilon }_D}}\right]\text{exp}\left[-\frac{{\zeta }_D^2}{4\left(1-v\right)\left(1-w\right)s{a}^2{t}_4^2}\right]+{\displaystyle {\int }_0^{1}ds{s}^{7/2}}{F}_0\left[z_2\right]\text{exp}[-\frac{u\left(1-u\right)w}{s}{a}^2{\overline{AB}}^2-\frac{v(1-v)(1-w)}{s}{a}^2{\overline{CD}}^2]$ $\left[\frac{s{\zeta }_A^2}{2uw{a}^2}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]\text{exp}\left[-\frac{s{\zeta }_A^2}{4uw{a}^2{t}_1^2}\right]\left[\frac{s{\zeta }_B^2}{2\left(1-u\right)w{a}^2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]\text{exp}\left[-\frac{s{\zeta }_B^2}{4\left(1-u\right)w{a}^2{t}_2^2}\right]$ $\left[\frac{s{\zeta }_C^2}{2v\left(1-w\right)a^2}{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{4+{\epsilon }_C}}-{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{2+{\epsilon }_C}}\right]\text{exp}\left[-\frac{s{\zeta }_C^2}{4v\left(1-w\right)a^2{t}_3^2}\right]\left[\frac{s{\zeta }_D^2}{2\left(1-v\right)\left(1-w\right)a^2}{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{4+{\epsilon }_D}}-{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{2+{\epsilon }_D}}\right]\text{exp}\left[-\frac{s{\zeta }_D^2}{4\left(1-v\right)\left(1-w\right)a^2{t}_4^2}\right]\}$ (5.18)

where $z_1=w\left(1-w\right)s{a}^2{\overline{PQ}}^2$ , $z_2=\frac{w(1-w)}{s}{a}^2{\overline{PQ}}^2$ , and ${\overline{PQ}}^2={\overline{BD}}^2+u^2{\overline{AB}}^2+v^2{\overline{CD}}^2+2u\overrightarrow{BD}\cdot \overrightarrow{AB}-2v\overrightarrow{BD}\cdot \overrightarrow{CD}-2uv\overrightarrow{AB}\overrightarrow{CD}$ .

Integrals over w, u, v, s, t¹, t², t3, and t4 can be evaluated by the Gauss-Legendre quadrature. The 64 point-quadrature can give a good precision of 7 SFs. The calculated value is ERI = 0.2343726 × e2 for the electron-repulsion integral with ζA = ζB = ζC = ζD = 1, $\overrightarrow{A}=(-\frac{\sqrt{8}}{3}, -\sqrt{\frac{8}{3}}, \frac{2}{3})$ , $\overrightarrow{B}=(-\frac{\sqrt{8}}{3}, \sqrt{\frac{8}{3}}, \frac{2}{3})$ , $\overrightarrow{C}=(\frac{\sqrt{24}}{3}, 0, \frac{2}{3})$ , and $\overrightarrow{D}=(0, 0, -2)$ , which is the case of four hydrogen atoms at $\overrightarrow{A}$ , $\overrightarrow{B}$ , $\overrightarrow{C}$ , $\overrightarrow{D}$ and $\overrightarrow{D}$ . For the case of four carbon +5 cations, the 128-point quadrature is necessary to give 5 SF precision. The value is as ERI = 0.61259(-13) × e2 with ζA = ζB =ζC = ζD = 6, and $\overrightarrow{A}$ , $\overrightarrow{B}$ , $\overrightarrow{C}$ and $\overrightarrow{D}$ are the same as in the case of the four hydrogen atoms. Because the absolute value is very small, the 5 SF precision may be satisfactory.

The author has been interested in whether the electron is not a point-like particle but a finite-sized one. If the charge density distribution of the finite-sized one is like the GCDD model of the nucleus, we can calculate the ERI of the finite-sized one. How the value of the ERI for the finite-sized electron is different from that of Eq. (5.18). Then we derive the ERI for the finite-sized electron in Appendix B (See Below). As seen in Appendix B, we find the value of the ERI for the finite-sized electron is the same as that for the point-like electron within 7 significant-figures precision for the case of four hydrogen atoms and within 5 significant-figure precision for four carbon +5 cations.

6. Conclusion

Using the Gaussian-transform for the Dirac wave function, which is derived in a previous article [18], we derive multicenter molecular integrals over Dirac wave functions for several fundamental properties; i.e., the overlap integral, the KEI, the NAI for both the point-like nucleus and the finite one, and the ERI. The value of the NAI for the finite nucleus with the GCDD model is the same as that for the point-like nucleus within 8 SF precision in the case of three hydrogen atoms and within 7 SF precision in the case of three carbon +5 cations. Also, the value of the ERI for the finite-sized electron (with the GCDD model) is the same as that for the point-like electron within 7 SF precision in the case of four hydrogen atoms and within 5 SF precision in the case of the four carbon +5 cations.

For solving the molecular matrix Dirac equation, we need more molecular integrals than those for the present fundamental properties. Such a project is in progress.

7. Acknowledgement

The author thanks a reviewer for many helpful suggestions to improve the explanation of the article.

Appendix

Appendix A. Gaussian-transform for the derivative of the Dirac wave function<

The derivative of the Dirac wave function centered at B is given by

$\nabla r_B^{-{\epsilon }_B}\exp \left[-{\zeta }_B{r}_B\right]=-\overrightarrow{r_B}\left(\frac{{\epsilon }_B}{r_B^{2+{\epsilon }_B}}+\frac{{\zeta }_B}{r_B^{1+{\epsilon }_B}}\right)\text{exp}[-{\zeta }_B{r}_B]$ (A1)

We know the identity given by [27]

$\frac{\text{exp}(-\beta )}{{\beta }^{\mu }}=\frac{1}{\Gamma (\mu )}{\displaystyle {\int }_0^{1}dt}\frac{{(1-t)}^{\mu -1}}{t^{\mu +1}}exp\left(-\frac{\beta }{t}\right)$ (A2)

Using Eq. (A2) with $\beta ={\zeta }_B{r}_B$ , we have

$\nabla r_B^{-{\epsilon }_B}\exp \left[-{\zeta }_B{r}_B\right]=-\overrightarrow{r_B}\frac{{\epsilon }_B{\zeta }_B^{2+{\epsilon }_B}}{\Gamma \left(2+{\epsilon }_B\right)}{\displaystyle {\int }_0^{1}dt}\frac{{(1-t)}^{1+{\epsilon }_B}}{t^{3+{\epsilon }_B}}exp\left(-\frac{{\zeta }_B{r}_B}{t}\right)-\overrightarrow{r_B}\frac{{\zeta }_B^{2+{\epsilon }_B}}{\Gamma (1+{\epsilon }_B)}{\displaystyle {\int }_0^{1}dt}\frac{{(1-t)}^{{\epsilon }_B}}{t^{2+{\epsilon }_B}}exp\left(-\frac{{\zeta }_B{r}_B}{t}\right)$ (A3)

We know the Gaussian-transform of 1s Slater-type orbital (STO) given by [28]

$\exp \left[-\zeta r\right]=\frac{\zeta }{2\sqrt{\pi }}{\displaystyle {\int }_0^{\infty }dS} S^{-3/2}exp\left(-S{r}^2-\frac{{\zeta }^2}{4S}\right)$ (A4)

Using Eq. (A4) with $\zeta =\frac{{\zeta }_B}{t}$ , we have the final formula of the Gaussian-transform for the derivative of the Dirac function as given by

$\nabla r_B^{-{\epsilon }_B}\exp \left[-{\zeta }_B{r}_B\right]=-\overrightarrow{r_B}\frac{{\zeta }_B^{3+{\epsilon }_B}}{2\sqrt{\pi }\Gamma \left(2+{\epsilon }_B\right)}{\displaystyle {\int }_0^{\infty }dS} S^{-\frac{3}{2}}\text{exp}[-S{r}_B^2]\left[{\epsilon }_B{\displaystyle {\int }_0^{1}dt}\frac{{(1-t)}^{{\epsilon }_B}}{t^{4+{\epsilon }_B}}+{\displaystyle {\int }_0^{1}dt}\frac{{(1-t)}^{{\epsilon }_B}}{t^{3+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4S{t}^2}\right]$ (A5)

where we use the relation given by

$\Gamma \left(2+{\epsilon }_B\right)=(1+{\epsilon }_B)\Gamma (1+{\epsilon }_B)$ (A6)

and

$\frac{(1-t){\epsilon }_B}{\left(1+{\epsilon }_B\right) t}+1=\frac{{\epsilon }_B+t}{\left(1+{\epsilon }_B\right)t}$ (A7)

Appendix B. Electron-repulsion integral for the finite-sized electron

We consider here in the case if the electron is a finite-sized particle and inter-electron potential is given by

$\frac{e^2}{r_{12}}\to \frac{e^2}{r_{12}}\text{erf}\left(\sqrt{{\xi }_e}{r}_{12}\right)=\frac{2e^2}{\sqrt{\pi }{r}_e}{F}_0\left(\frac{r_{12}^2}{r_e^2}\right)=\frac{2e^2\sqrt{{\xi }_e}}{\sqrt{\pi }}{F}_0({\xi }_e{r}_{12}^2)$ (B1)

where each of ${\xi }_e$ and re is the size parameter and related to each other as $r_e=1/\sqrt{{\xi }_e}$ . We may use the classical radius of the electron as re = 2.81794092(−15) m = 0.532513619(−4) bohr. The value of the size parameter may be somewhat smaller than this value because the size parameter r0 of the GCDD model of the finite nucleus is somewhat smaller than the experimental radius of the finite nucleus as given by $r_0=\sqrt{2/3}RMS$ , where RMS is the experimental root-mean-square radius [21]. Using the potential, Eq. (B1), we derive the electron-repulsion integral (ERI) for the finite-sized electron as follows:

$ER{I}^{(F)}=e^2{N}_A{N}_B{N}_C{N}_D{V}_{ABCD}^{(F)}$ (B2)

where

$\text{exp}[-{\zeta }_A{r}_{1A}-{\zeta }_B{r}_{1B}-{\zeta }_C{r}_{2C}-{\zeta }_D{r}_{2D}]$ (B3)

Using Eq. (1.10). we have

$V_{ABCD}^{(F)}=\frac{1}{{(4\pi )}^2}\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}{\zeta }_C^{1+{\epsilon }_C}{\zeta }_D^{1+{\epsilon }_D}}{\Gamma \left(1+{\epsilon }_A\right)\Gamma \left(1+{\epsilon }_B\right)\Gamma \left(1+{\epsilon }_C\right)\Gamma \left(1+{\epsilon }_D\right)}{\displaystyle {\int }_0^{\infty }d{S}_1}{\displaystyle {\int }_0^{\infty }d{S}_2}{\displaystyle {\int }_0^{\infty }d{S}_3}{\displaystyle {\int }_0^{\infty }d{S}_4}{(S_1{S}_2{S}_3{S}_4)}^{-3/2}$ $\left[\frac{{\zeta }_A^2}{2S_1}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4S_1{t}_1^2}\right]\left[\frac{{\zeta }_B^2}{2S_2}{{\displaystyle {\int }_0^{1}dt}}_2\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{{\displaystyle {\int }_0^{1}dt}}_2\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4S_2{t}_2^2}\right]$ $\left[\frac{{\zeta }_C^2}{2S_3}{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{4+{\epsilon }_C}}-{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{2+{\epsilon }_C}}\right]exp\left[-\frac{{\zeta }_C^2}{4S_3{t}_3^2}\right]\left[\frac{{\zeta }_D^2}{2S_4}{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{4+{\epsilon }_D}}-{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{2+{\epsilon }_D}}\right]exp\left[-\frac{{\zeta }_D^2}{4S_4{t}_4^2}\right]I_1^{(F)}$ (B4)

where

$I_1^{(F)}=\frac{2\sqrt{{\xi }_e}}{\sqrt{\pi }}{{\displaystyle \int }}^{\text{}}d\overrightarrow{r_1} d\overrightarrow{r_2} F_0({\xi }_e{r}_{12}^2)\text{exp}[-S_1{r}_{1A}^2-S_2{r}_{1B}^2-S_3{r}_{2C}^2-S_4{r}_{2D}^2]$ (B5)

We use the Gaussian product rule and have

$I_1^{\left(F\right)}=exp\left[-\frac{S_1{S}_2}{S_{12}}{\overline{AB}}^2-\frac{S_3{S}_4}{S_{34}}C{D}^2\right]I_2^{(F)}$ (B6)

where

$I_2^{(F)}=\frac{2\sqrt{{\xi }_e}}{\sqrt{\pi }}{{\displaystyle \int }}^{\text{}}d\overrightarrow{r_1}d\overrightarrow{r_2}{F}_0({\xi }_e{r}_{12}^2)\text{exp}[-S_{12}{r}_{1P}^2-S_{34}{r}_{2Q}^2]$ (B7)

We know $F_0\left({\xi }_e{r}_{12}^2\right)={\displaystyle {\int }_0^{1}dt}\text{exp}[-{\xi }_e{r}_{12}^2{t}^2]$ (B8)

Then we have

$I_2^{(F)}=\frac{2\sqrt{{\xi }_e}}{\sqrt{\pi }}{\displaystyle {\int }_0^{1}dt}{{\displaystyle \int }}^{\text{}}d\overrightarrow{r_1}\exp \left[-S_{12}{r}_{1P}^2\right]I_3^{(F)}$ (B9)

where

$I_3^{(F)}={{\displaystyle \int }}^{\text{}}d\overrightarrow{r_2}\text{exp}[-S_{34}{r}_{2Q}^2-{\xi }_e{t}^2{r}_{12}^2]$ (B10)

We use the Gaussian product rule and have

$I_3^{(F)}=exp\left[-\frac{S_{34}{\xi }_e{t}^2}{S_{34}+{\xi }_e{t}^2}{r}_{1Q}^2\right]{{\displaystyle \int }}^{\text{}}d\overrightarrow{r_2}\text{exp}[-(S_{34}+{\xi }_e{t}^2)r_{2R}^2]=exp\left[-\frac{S_{34}{\xi }_e{t}^2}{S_{34}+{\xi }_e{t}^2}{r}_{1Q}^2\right]{\left(\frac{\pi }{S_{34}+{\xi }_e{t}^2}\right)}^{3/2}$ (B11)

where $\overrightarrow{R}=\frac{s_{34}}{S_{34}+{\xi }_e{t}^2}\overrightarrow{Q}+\frac{{\xi }_e{t}^2}{S_{34}+{\xi }_e{t}^2}\overrightarrow{r_1}$ . Substituting Eq. (B11) into Eq. (B9), we have

$I_2^{\left(F\right)}=\frac{2\sqrt{{\xi }_e}}{\sqrt{\pi }}{\displaystyle {\int }_0^{1}dt}{\left(\frac{\pi }{S_{34}+{\xi }_e{t}^2}\right)}^{3/2}{{\displaystyle \int }}^{\text{}}d\overrightarrow{r_1}exp\left[-S_{12}{r}_{1P}^2-\frac{S_{34}{\xi }_e{t}^2}{S_{34}+{\xi }_e{t}^2}{r}_{1Q}^2\right]$ (B12)

We use the Gaussian product rule and have

$I_2^{(F)}=\frac{2\sqrt{{\xi }_e}}{\sqrt{\pi }}{\displaystyle {\int }_0^{1}dt}{\left(\frac{\pi }{S_{34}+{\xi }_e{t}^2}\right)}^{3/2}exp\left[-\frac{S_{12}{S}_{34}{\xi }_e{t}^2}{\delta (S_{34}+{\xi }_e{t}^2)}{\stackrel{-}{PQ}}^2\right]{{\displaystyle \int }}^{\text{}}d\overrightarrow{r_1}\text{exp}[-\delta r_{1S}^2]$ $=\frac{2\sqrt{{\xi }_e}}{\sqrt{\pi }}{\displaystyle {\int }_0^{1}dt}{\left(\frac{\pi }{S_{34}+{\xi }_e{t}^2}\right)}^{3/2}{\left(\frac{\pi }{\delta }\right)}^{3/2}exp\left[-\frac{S_{12}{S}_{34}{\xi }_e{t}^2}{\delta (S_{34}+{\xi }_e{t}^2)}{\overline{PQ}}^2\right]=\frac{2\sqrt{{\xi }_e}{\pi }^3}{\sqrt{\pi }}{\displaystyle {\int }_0^{1}dt}\frac{1}{{[S_{12}{S}_{34}+S_{1234}{\xi }_e{t}^2]}^{3/2}}exp\left[-\frac{S_{12}{S}_{34}{\xi }_e{t}^2{\overline{PQ}}^2}{S_{12}{S}_{34}+S_{1234}{\xi }_e{t}^2}\right]$ (B13)

where $S_{1234}=S_{12}+S_{34}$ , $\delta =\frac{S_{12}{S}_{34}+S_{1234}{\xi }_e{t}^2}{S_{34}+{\xi }_e{t}^2}$ and $\overrightarrow{S}=\frac{S_{12}}{\delta }\overrightarrow{P}+\frac{S_{34}{\xi }_e{t}^2}{S_{12}{S}_{34}+S_{1234}{\xi }_e{t}^2}\overrightarrow{Q}$ . We use $c_1=S_{12}{S}_{34}$ and $c_2=S_{1234}{\xi }_e$ and have

$I_2^{(F)}=\frac{2\sqrt{{\xi }_e}{\pi }^3}{\sqrt{\pi }}{\displaystyle {\int }_0^{1}dt}\frac{1}{{(c_1+c_2{t}^2)}^{3/2}}exp\left[-\frac{c_1{\xi }_e{t}^2{\overline{PQ}}^2}{c_1+c_2{t}^2}\right]$ (B14)

We change the integral variable as $x^2=\frac{(c_1+c_2)t^2}{C_1+C_2{t}^2}$ as x → 0 for t → 0 and x → 1 for t → 1. We know $dt=\frac{c_1+c_2}{c_1}\frac{t^3}{x^3}dx$ and have

$I_2^{(F)}=\frac{2\sqrt{{\xi }_e}{\pi }^3}{\sqrt{\pi }}{\displaystyle {\int }_0^{1}dx}\frac{(c_1+c_2)t^3}{c_1{x}^3}{\left(\frac{x^2}{(c_1+c_2)t^2}\right)}^{3/2}exp\left[-\frac{c_1{\xi }_e}{c_1+c_2}{\overline{PQ}}^2{x}^2\right]=\frac{2\sqrt{{\xi }_e}{\pi }^3}{\sqrt{\pi }}\frac{1}{c_1\sqrt{c_1+c_2}}{\displaystyle {\int }_0^{1}dx}exp\left[-z_0{x}^2\right]=\frac{2\sqrt{{\xi }_e}{\pi }^3}{\sqrt{\pi }}\frac{1}{c_1\sqrt{c_1+c_2}}{F}_0(z_0)$ (B15)

where $z_0=\frac{c_1}{c_1+c_2}{\xi }_e{\overline{PQ}}^2$ . We change the integral variables as $z=S_{1234}$ , $w=\frac{S_{12}}{S_{1234}}$ , $u=\frac{S_1}{S_{12}}$ , and $v=\frac{S_3}{S_{34}}$ . The Jacobian is $\frac{\partial (S_1,S_2,S_3,S_4)}{\partial (z, w, u, v)}=w(1-w)z^3$ . Then, after substituting Eq. (B15) into Eq. (B6) and doing Eq. (B6) into Eq. (B4), we have

$V_{ABCD}^{(F)}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}{\zeta }_C^{1+{\epsilon }_C}{\zeta }_D^{1+{\epsilon }_D}}{\Gamma \left(1+{\epsilon }_A\right)\Gamma \left(1+{\epsilon }_B\right)\Gamma \left(1+{\epsilon }_C\right)\Gamma \left(1+{\epsilon }_D\right)}{\displaystyle {\int }_0^{1}dw}{\displaystyle {\int }_0^{1}du}{\displaystyle {\int }_0^{1}dv}{\displaystyle {\int }_0^{\infty }dz}\frac{\text{exp}[-u\left(1-u\right)wz{\overline{AB}}^2-v\left(1-v\right)\left(1-w\right)z{\overline{CD}}^2]}{{[u\left(1-u\right)v(1-v)]}^{3/2}{[w(1-w)]}^2{z}^3}$ $\frac{2\sqrt{{\xi }_e}}{\sqrt{\pi }}\frac{{\pi }^3}{{(4\pi )}^2}\frac{F_0(z_0)}{w(1-w)z^2\sqrt{w\left(1-w\right)z^2+{\xi }_{e}z}}\left[\frac{{\zeta }_A^2}{2uwz}{\displaystyle {\int }_0^{1}dt}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}dt}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4uwz{t}_1^2}\right]$ $\left[\frac{{\zeta }_B^2}{2\left(1-u\right)wz}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4\left(1-u\right)wz{t}_2^2}\right]\left[\frac{{\zeta }_C^2}{2v\left(1-w\right)z}{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{4+{\epsilon }_C}}-{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{2+{\epsilon }_C}}\right]exp\left[-\frac{{\zeta }_C^2}{4v\left(1-w\right)z{t}_3^2}\right]$ $\left[\frac{{\zeta }_D^2}{2\left(1-v\right)\left(1-w\right)z}{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{4+{\epsilon }_D}}-{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{2+{\epsilon }_D}}\right]exp\left[-\frac{{\zeta }_D^2}{4\left(1-v\right)\left(1-w\right)z{t}_4^2}\right]$ (B16)

We separate the integral over z as similarly to Eq. (2.8) as given by

${\displaystyle {\int }_0^{\infty }dz}=a^2{\displaystyle {\int }_0^{1}ds}+a^2{\displaystyle {\int }_0^{1}ds}\frac{1}{s^2}$ (B17)

Then we have the final formula of the electron-repulsion integral over Dirac wave functions for the finite-sized electron as given by

$V_{ABCD}^{(F)}=\frac{{\zeta }_A^{1+{\epsilon }_A}{\zeta }_B^{1+{\epsilon }_B}{\zeta }_C^{1+{\epsilon }_C}{\zeta }_D^{1+{\epsilon }_D}}{\Gamma \left(1+{\epsilon }_A\right)\Gamma \left(1+{\epsilon }_B\right)\Gamma \left(1+{\epsilon }_C\right)\Gamma \left(1+{\epsilon }_D\right)}\frac{\sqrt{\pi }}{8a^9}{\displaystyle {\int }_0^{1}dw}{\displaystyle {\int }_0^{1}du}{\displaystyle {\int }_0^{1}dv}{[u\left(1-u\right)v(1-v)]}^{-3/2}{[w(1-w)]}^{-3}$ $\{{\displaystyle {\int }_0^{1}ds}\frac{\exp \left[-u\left(1-u\right)ws{a}^2{\overline{AB}}^2-v\left(1-v\right)\left(1-w\right)s{a}^2{\overline{CD}}^2\right]}{s^{11/2}{[1+w\left(1-w\right)s{a}^2{r}_e^2]}^{1/2}}{F}_0({\delta }_1{\overline{PQ}}^2)\left[\frac{{\zeta }_A^2}{2uws{a}^2}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{{\zeta }_A^2}{4uws{a}^2{t}_1^2}\right]$ $\left[\frac{{\zeta }_B^2}{2\left(1-u\right)ws{a}^2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{{\zeta }_B^2}{4\left(1-u\right)ws{a}^2{t}_2^2}\right]\left[\frac{{\zeta }_C^2}{2v\left(1-w\right)s{a}^2}{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{4+{\epsilon }_C}}-{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{2+{\epsilon }_C}}\right]exp\left[-\frac{{\zeta }_C^2}{4v\left(1-w\right)s{a}^2{t}_3^2}\right]$ $\left[\frac{{\zeta }_D^2}{2\left(1-v\right)\left(1-w\right)s{a}^2}{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{4+{\epsilon }_D}}-{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{2+{\epsilon }_D}}\right]exp\left[-\frac{{\zeta }_D^2}{4\left(1-v\right)\left(1-w\right)s{a}^2{t}_4^2}\right]+{\displaystyle {\int }_0^{1}ds}\frac{\exp \left[-u\left(1-u\right)w{a}^2{\overline{AB}}^2/s-v\left(1-v\right)\left(1-w\right)a^2{\overline{CD}}^2/s\right]}{{[1+w\left(1-w\right)a^2{r}_e^2/s]}^{1/2}}{s}^{7/2}{F}_0({\delta }_1{\overline{PQ}}^2)$ $\left[\frac{s{\zeta }_A^2}{2uw{a}^2}{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{4+{\epsilon }_A}}-{\displaystyle {\int }_0^{1}d{t}_1}\frac{{(1-t_1)}^{{\epsilon }_A}}{t_1^{2+{\epsilon }_A}}\right]exp\left[-\frac{s{\zeta }_A^2}{4uw{a}^2{t}_1^2}\right]\left[\frac{s{\zeta }_B^2}{2\left(1-u\right)w{a}^2}{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{4+{\epsilon }_B}}-{\displaystyle {\int }_0^{1}d{t}_2}\frac{{(1-t_2)}^{{\epsilon }_B}}{t_2^{2+{\epsilon }_B}}\right]exp\left[-\frac{s{\zeta }_B^2}{4\left(1-u\right)w{a}^2{t}_2^2}\right]$ $\left[\frac{s{\zeta }_C^2}{2v\left(1-w\right)a^2}{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{4+{\epsilon }_C}}-{\displaystyle {\int }_0^{1}d{t}_3}\frac{{(1-t_3)}^{{\epsilon }_C}}{t_3^{2+{\epsilon }_C}}\right]exp\left[-\frac{s{\zeta }_C^2}{4v\left(1-w\right)a^2{t}_3^2}\right]\left[\frac{s{\zeta }_D^2}{2\left(1-v\right)\left(1-w\right)a^2}{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{4+{\epsilon }_D}}-{\displaystyle {\int }_0^{1}d{t}_4}\frac{{(1-t_4)}^{{\epsilon }_D}}{t_4^{2+{\epsilon }_D}}\right]exp\left[-\frac{s{\zeta }_D^2}{4\left(1-v\right)\left(1-w\right)a^2{t}_4^2}\right]\}$ (B18)

where ${\delta }_1=\frac{w\left(1-w\right)s{a}^2{\xi }_e}{w\left(1-w\right)s{a}^2+{\xi }_e}=\frac{w\left(1-w\right)s{a}^2}{1+w\left(1-w\right)s{a}^2{r}_e^2}$ and ${\delta }_2=\frac{w(1-w)a^2{\xi }_e}{w\left(1-w\right)a^2+s{\xi }_e}=\frac{w(1-w)a^2}{s+w(1-w)a^2{r}_e^2}$ . When re → 0 in Eq. (B18), we see that Eq. (B18) coincides with Eq. (5.18). Integrals over w, u, v, s, t¹, t², t3, and t4 can be evaluated by the Gauss-Legendre quadrature. The 64 point-quadrature can give a good precision of 7 SFs. The calculated value is as ERI(F) = 0.2343726 × e2 for the finite-sized electron with ${\zeta }_A={\zeta }_B={\zeta }_C={\zeta }_D=1$ , $\overrightarrow{A}=(-\frac{\sqrt{8}}{3}, -\sqrt{\frac{8}{3}}, \frac{2}{3})$ , $\overrightarrow{B}=(-\frac{\sqrt{8}}{3}, \sqrt{\frac{8}{3}}, \frac{2}{3})$ , $\overrightarrow{C}=(\frac{\sqrt{24}}{3}, 0, \frac{2}{3})$ , and $\overrightarrow{D}=(0, 0, -2)$ , which is the case of four hydrogen atoms at $\overrightarrow{A}$ , $\overrightarrow{B}$ , $\overrightarrow{C}$ , and $\overrightarrow{D}$ . Thus, the value is the same as that for the point-like electron within the 7 SF precision, because $r_e^2$ = 0.2835707544(−8)bohr2 is very small. For the case of four carbon +5 cations, the 128-point quadrature is necessary to give the 5 SF precision. The value is as ERI(F) = 0.61259(−13) × e2 with ${\zeta }_A={\zeta }_B={\zeta }_C={\zeta }_D=6$ , and $\overrightarrow{A}$ , $\overrightarrow{B}$ , $\overrightarrow{C}$ , and $\overrightarrow{D}$ are the same as in the four-hydrogen case. Thus, the value is also the same as that for the point-like electron within the 5 SF precision.

記事について