# Implementation of one-elecron integrals in a Hartree-Fock program

Absolute newbie in quantum chemistry here. I was recently looking through some Fortran libraries from my university to see if I understand all of the mathematics. In a Hartree-Fock library, I found the following code snippet for calculation of the overlap-integral between two contracted s-type GTOs.

 rab = sum( (r_a-r_b)**2 )

do i=1,npa
do j=1,npb
eab = alp(i)+bet(j)
oab = 1.0_wp/eab
cab = ci(i)*cj(j)
xab = alp(i)*bet(j)*oab
est = rab*xab
ab = exp(-est)
s00 = cab*ab*sqrt(pi*oab)**3

!overlap
sab = sab+s00
enddo
enddo


with r_a and r_b being the positions of the two gaussians in cartesian coordinates and alp(i) and bet(i) containing the orbital coefficients of the primitives. From this code snippet I read that the overlap-integral is calculated via: $$K_{AB}=e^{-r_{AB}(\frac{\alpha \beta}{\alpha+\beta})}\tag{1}$$ $$S= K_{AB} \cdot \sqrt{\frac{\pi}{\alpha+\beta}}^3\tag{2}$$

But if I try to get to this formula myself, I get a slightly different result:

For the overlap integral of two Gaussians $$\phi_{A}$$ and $$\phi_{B}$$

$$\phi_{A}(\alpha,r-R_A)=e^{-\alpha(r-R_A)^2}\tag{3}$$

$$\phi_{B}(\beta,r-R_B)=e^{-\beta(r-R_B)^2}\tag{4}$$

we can use the Gaussian-product-theorem to combine both Gaussians to one function $$\phi_{AB}$$: $$\phi_{A}(\alpha,r-R_A) \cdot \phi_{B}(\beta,r-R_B) = K_{AB} \cdot \phi_{AB}(\alpha+\beta,r-\frac{\alpha R_A + \beta R_B}{\alpha + \beta})\tag{5}$$ $$\phi_{AB} = K_{AB} \cdot \left (\frac{2\alpha\beta}{(\alpha+\beta)\pi} \right )^{\frac{3}{4}} \cdot e^{-\alpha\beta(r-\frac{\alpha R_A+\beta R_B}{\alpha+\beta})^2}\tag{6}$$

Now integrating this Gaussian should yield the overlap integral between our starting functions. The integral of an arbitrary gaussian can be solved via: $$\int \phi_A(\alpha ,r-R_A) = \sqrt{\frac{\pi}{\alpha}}\tag{7}$$ $$\tag{8} \int\phi_{AB}=K_{AB} \cdot \left (\frac{2\alpha\beta}{(\alpha+\beta)\pi}\right )^{\frac{3}{4}} \cdot \sqrt{\frac{\pi}{\alpha\beta}}$$

Which, as far as I can see, is already a different equation to the one I read from the code. Where is my mistake? I know this is some pretty basic stuff, any help is appreciated.

There is a bunch of errors in your question.

First of all you seem to ignore the contraction coefficients in your code snippet. If however you ignore them (set them to 1 and pick only one pair $$i,j$$), still, your code calculates:

$$S_{00} = e^{-\left (|\vec R_{AB}|^2 \frac{\alpha \beta}{\alpha+\beta}\right ) } \sqrt{\frac{\pi}{\alpha+\beta}}^3\tag{1}$$

which means you missed the square in what you call $$K_{AB}$$

Then you try to compare to a "by hand" calculation that you do yourself, for which you use the Gaussian product rule in principle correctly stating that:

$$\phi_{A}(\alpha,r-R_A) \cdot \phi_{B}(\beta,r-R_B) = K_{AB} \cdot \phi_{AB}\left(\alpha+\beta,r-\frac{\alpha R_A + \beta R_B}{\alpha + \beta}\right)\tag{2}$$

However, in the next line, you write

$$\phi_{AB} = K_{AB} \cdot \left (\frac{2\alpha\beta}{(\alpha+\beta)\pi} \right )^{\frac{3}{4}} \cdot e^{-\alpha\beta\left(r-\frac{\alpha R_A+\beta R_B}{\alpha+\beta}\right)^2}\tag{3}$$

which suffers from three issues, for one your definition of $$K_{AB}$$ is incorrect since you missed the square, additionally your pre-exponential factor is suddenly $$\alpha \beta$$ instead of $$\alpha + \beta$$. Moreover, it is unclear where

$$\left (\frac{2\alpha\beta}{(\alpha+\beta)\pi} \right )^{\frac{3}{4}}\tag{4}$$

is coming from.

To the best of my knowledge the gaussian product rule should have resulted in:

$$\phi_{AB} = K_{AB} e^{-(\alpha + \beta)r_P}\tag{5}$$

where $$\vec r_P$$ is $$\vec r - \vec P$$ and $$\vec P$$ is the place: $$\begin{eqnarray} \vec P = \frac{\alpha \vec R_A - \beta \vec R_B}{\alpha + \beta}\tag{6} \end{eqnarray}$$ and $$K_{AB} = e^{-\frac{\alpha \; \beta }{\alpha + \beta} \vec R_{AB}^2}\tag{7}$$

You than carry these errors into your result.

I guess all that was needed was some double checking on your formulas.