## RHF with manually computed integrals over s functions

In this project, we will compute the RHF energy of the diatomic hydrogen molecule using a basis composed of purely s-type Gaussian orbitals. Instead of relying on Psi4 for integrals, we will compute them manually. As we learned in class, all higher angular momentum integrals can be constructed from a linear combination of integrals over s functions using recursion relations. After completing this assignment, you should have a good understanding of how electronic structure codes compute integrals; you already dealt with integrals over s-functions, so if you _were_ to implement a recursion function, you would have a fully functional integrals code.

We will be attempting to mimic the integral and RHF energy outputs of the following Psi4 input file:

```
import psi4  
import numpy as np  
np.set_printoptions(linewidth=300)  
  
basis {  
assign s_orb  
[ s_orb ]  
cartesian  
****  
H     0  
S   1   1.00  
      0.50000000000      1.0000000000  
S   1   1.00  
      0.40000000000      1.0000000000  
S   1   1.00  
      0.30000000000      1.0000000000  
S   1   1.00  
      0.20000000000      1.0000000000  
****  
}  
  
molecule h2 {  
symmetry c1  
units bohr  
0 1  
H 0.000000000000 0.000000000000 -0.849220457955  
H 0.000000000000 0.000000000000  0.849220457955  
}  
  
set scf_type pk          # no density fitting  
set puream false         # cartesian basis (not spherical)  
set e_convergence 10  
set d_convergence 10  
e, wfn = energy('hf', return_wfn=True)  
  
print("Hartree-Fock energy: ", e)  
  
mol = wfn.molecule()  
bs = wfn.basisset()  
basis = psi4.core.BasisSet.build(mol)  
mints = psi4.core.MintsHelper(basis)  

S = mints.ao_overlap().np  
T = mints.ao_kinetic().np  
V = mints.ao_potential().np  
G = mints.ao_eri().np  

#print(S)  
#print(T)  
#print(V)  
#print(G)  
```


At first define a custom basis composed of 4 s-functions. They have different orbital exponents $\alpha$, but they all are composed of a single primitive (there are no contractions here, the contraction coefficient of each is 1). We then just define the geometry in bohr, do a Hartree-Fock computation, and then load in the integrals. We do this so that we can print and visually check Psi4's integral arrays with our own later on (by uncommenting the print statements at the bottom and running `psi4 input.dat` and the command line).

We will start by importing scipy and numpy. This programming project will receive point deductions if you use anything other than scipy and numpy.


```python
import numpy as np
from scipy import special, linalg
np.set_printoptions(linewidth=200)
```

Here we will construct arrays to mimic the Cartesian geometry and basis set information. We also need the charge of each atom for potential energy integrals and the nuclear repulsion energy.


```python
# Create a geometry in bohr which matches our psi4 input file
geom = np.array([[0.000000000000,0.000000000000,-0.849220457955],
                 [0.000000000000,0.000000000000, 0.849220457955]])
# Define charge, number of basis functions per atom, organize basis function data
charge = np.array([1.0,1.0])
nbf = 8
nbf_per_atom = np.array([4,4])
centers = np.repeat(geom, nbf_per_atom, axis=0)
exponents = np.tile(np.array([0.5,0.4,0.3,0.2]), 2)

```


```python
print(centers)
print(exponents.reshape(-1,1))
```

    [[ 0.          0.         -0.84922046]
     [ 0.          0.         -0.84922046]
     [ 0.          0.         -0.84922046]
     [ 0.          0.         -0.84922046]
     [ 0.          0.          0.84922046]
     [ 0.          0.          0.84922046]
     [ 0.          0.          0.84922046]
     [ 0.          0.          0.84922046]]
    [[ 0.5]
     [ 0.4]
     [ 0.3]
     [ 0.2]
     [ 0.5]
     [ 0.4]
     [ 0.3]
     [ 0.2]]


We are now ready to make some functions for computing integrals. First of all, we need a function which computes normalization constants for each Gaussian.


```python
def normalize(aa):
    '''Normalization constant for s-type primitive Gaussian basis functions. Argument is orbital exponent coefficient'''
    N = ((2 * aa) / np.pi)**(3/4)
    return N

```

We will also need at some point the 0th Boys Function for the electron-nuclear attraction integrals and electron repulsion integrals. Recall from the lecture notes that the 0th Boys function can be found by:

\begin{equation}
F_0(x) = \mathrm{erf}(\sqrt x) \frac{\sqrt{\pi}} {2 \sqrt {x}} 
\end{equation}

However, this equation numerically blows up when $x$ is near 0, so we need to employ a Taylor expansion for small x:
\begin{equation}
1 - \frac{x}{3} + \frac{x^2}{10} - \frac{x^3}{42}
\end{equation}


```python
def boys0(x):
    '''0th Boys Function'''
    if x < 1e-8:
        return 1 - x/3 + x**2/10 - x**3/42
    else:
        return special.erf(np.sqrt(x)) * np.sqrt(np.pi) /  (2 * np.sqrt(x))
```

If you're feeling sassy, you can instead use a more general function which calls "Kummer's confluent hypergeometric function" $M$ (`scipy.special.hyp1f1`) using the relation

\begin{equation}
F_n(x) = \frac{M(n + \frac{1}{2}, n + \frac{3}{2}, -x)}{2n + 1}
\end{equation}


```python
def boys_general(n, x):
    '''F_n(x) Boys Function of any order '''
    denom = 2 * n + 1
    num = special.hyp1f1(n+0.5,n+1.5,-x)
    return num / denom

print(boys_general(0, 0.5))
print(boys0(0.5))

print(boys_general(0, 0.0))
print(boys0(0.0))

print(boys_general(0, 25.0))
print(boys0(25.0))
```

    0.855624391892
    0.855624391892
    1.0
    1.0
    0.17724538509
    0.17724538509


You can also use the definition in terms of the Gamma functions in the notes. 

## Overlap Integrals

The overlap integral over two _s_ functions from the lecture notes is 

\begin{equation}
(s \mid s) = N_a N_b \left(\frac{\pi}{\alpha_a + \alpha_b}\right)^{3/2} \exp \left( \frac{-\alpha_a \alpha_b (\boldsymbol{A} - \boldsymbol{B})^2}{\alpha_a + \alpha_b} \right)
\end{equation}

Reading off this equation, we can write a Python function which computes this as a function of the two centers and two orbital exponents.


```python
def overlap(A, B, aa, bb):
    Na = normalize(aa)
    Nb = normalize(bb)
    return Na * Nb * (np.pi / (aa + bb))**(3/2) * np.exp((-aa * bb * np.dot(A-B,A-B)) / (aa + bb))

```

To construct the overlap matrix, we need to compute the overlap integral between every possible basis function with every other possible basis function. We can do this with two loops over the number of basis functions, using our `exponents` and `centers` arrays from before:


```python
def compute_overlap(exponents, centers, nbf):
    S = np.zeros((nbf,nbf))
    for i in range(nbf):
        for j in range(nbf):
            A = centers[i]
            B = centers[j]
            aa = exponents[i]
            bb = exponents[j]
            S[i,j] = overlap(A, B, aa, bb)
    return S

```


```python
S = compute_overlap(exponents, centers, nbf)
print(S)
```




    array([[ 1.        ,  0.99072638,  0.9527489 ,  0.85881166,  0.48618047,  0.521857  ,  0.55472353,  0.56875481],
           [ 0.99072638,  1.        ,  0.9846545 ,  0.91545206,  0.521857  ,  0.56161411,  0.60050475,  0.62315241],
           [ 0.9527489 ,  0.9846545 ,  1.        ,  0.96984744,  0.55472353,  0.60050475,  0.6487517 ,  0.68606652],
           [ 0.85881166,  0.91545206,  0.96984744,  1.        ,  0.56875481,  0.62315241,  0.68606652,  0.74940917],
           [ 0.48618047,  0.521857  ,  0.55472353,  0.56875481,  1.        ,  0.99072638,  0.9527489 ,  0.85881166],
           [ 0.521857  ,  0.56161411,  0.60050475,  0.62315241,  0.99072638,  1.        ,  0.9846545 ,  0.91545206],
           [ 0.55472353,  0.60050475,  0.6487517 ,  0.68606652,  0.9527489 ,  0.9846545 ,  1.        ,  0.96984744],
           [ 0.56875481,  0.62315241,  0.68606652,  0.74940917,  0.85881166,  0.91545206,  0.96984744,  1.        ]])



Running the above psi4 input file and printing the overlap matrix S with `print(S)` we find it is the same array as what we have just generated.


The rest will be left for you to figure out. The lecture notes contain the formulas you need. Once you code how to compute a single kinetic integral, a single potential integral, and a single electron repulsion integral, you can write functions for constructing the entire arrays T, V, and I that we would normally load in from Psi4. Once you have this done and they appear to match Psi4, just copy and paste your RHF code.  You will need to figure out how to compute the nuclear repulsion energy on your own. Once finished, your energy should match Psi4. If it does, you pass, if it doesn't, you fail. If your code imports Psi4 anywhere, you fail.


```python

```

ff96dde5d8459d05a64a278e458596133ece9807
ipynb
Jupyter Notebook
2020/programming/3_RHF_from_scratch.ipynb
GQCG-edu/chem-8950
a5f58a5feacbae16b02fddd2c74723da1486b8d2
Apache-2.0
2017-09-29T20:25:50.000Z
2022-03-31T17:40:32.000Z
2020/programming/3_RHF_from_scratch.ipynb
GQCG-edu/chem-8950
a5f58a5feacbae16b02fddd2c74723da1486b8d2
Apache-2.0
2020-07-13T12:11:58.000Z
2020-07-13T15:31:47.000Z
2020/programming/3_RHF_from_scratch.ipynb
GQCG-edu/chem-8950
a5f58a5feacbae16b02fddd2c74723da1486b8d2
Apache-2.0
2017-09-10T10:33:28.000Z
2021-10-05T07:39:50.000Z
Qwen/Qwen-72B
1. YES
2. YES