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

Question

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.

Answer 1

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.