Computing molecular integrales in STO-3G basis with the Obara Saika recursion formulas

I'm trying to implement the Hartree-Fock algorithm on a water molecule using the STO-3G basis. This far I have only computed the overlap integrals with the Obara-Saika formulas but I'm running into problems with the ones involving p-orbitals. My question is,on formulas like this one,


$$
(p\_i||s) = (P\_i-A\_i)(s||s)
$$


What function is exactly the $(s|$ on the right hind side? My guess would be that it is the same as $(p\_i|$ (same exponent and coefficient) but with $l\_i = 0$ instead of 1. In case this is indeed the correct interpretation I'll leave the code I have written, and hopefully my mistake is really obvious.



```
import sys
import math
import numpy as np
from scipy import special
from scipy import linalg

def d_fact(n):
    product = 1.0
    for i in range(3,n+1,2):
        product *= i
    return product

class primitive_gaussian():

    def __init__(self, alpha, coeff, coordinates, l1, l2, l3):
        
        self.alpha = alpha
        self.coeff = coeff
        self.coordinates = np.array(coordinates)
        self.l = [l1,l2,l3]
        self.N = (2.0*alpha/math.pi)**(3/4)*(4*alpha)**((l1+l2+l3)/2)*(d_fact(2*l1-1)*d_fact(2*l2-1)*d_fact(2*l3-1))**(-1/2)

def overlap(i,j):

    N = i.N*j.N
    cicj = i.coeff*j.coeff
                    
    zeta = i.alpha+j.alpha
    xi = i.alpha*j.alpha/zeta
    AB = i.coordinates-j.coordinates
    AB2 = np.dot(AB,AB)

    P = (i.alpha*i.coordinates+j.alpha*j.coordinates)/zeta
    
    if((np.count_nonzero(i.l)==0)and(np.count_nonzero(j.l)==0)):
        return N*cicj*math.exp(-xi*AB2)*(math.pi/zeta)**(3/2)
    elif((np.count_nonzero(i.l)==1)and(np.count_nonzero(j.l)==0)):
        i_m = np.nonzero(i.l)
        si = primitive_gaussian(i.alpha,i.coeff,i.coordinates,0,0,0)
        return (P[i_m]-i.coordinates[i_m])*overlap(si,j)
    elif((np.count_nonzero(i.l)==0)and(np.count_nonzero(j.l)==1)):
        j_m = np.nonzero(j.l)
        sj = primitive_gaussian(j.alpha,j.coeff,j.coordinates,0,0,0)
        return (P[j_m]-j.coordinates[j_m])*overlap(i,sj)
    elif((np.count_nonzero(i.l)==1)and(np.count_nonzero(j.l)==1)):
        i_m = np.nonzero(i.l)
        j_m = np.nonzero(j.l)
        si = primitive_gaussian(i.alpha,i.coeff,i.coordinates,0,0,0)
        sj = primitive_gaussian(j.alpha,j.coeff,j.coordinates,0,0,0)
        if(i_m==j_m):
            return (P[j_m]-j.coordinates[j_m])*overlap(i,sj)+overlap(si,sj)/(2*zeta)
        else:
            return (P[j_m]-j.coordinates[j_m])*overlap(i,sj)

def kinetic(i,j):

    zeta = i.alpha+j.alpha
    xi = i.alpha*j.alpha/zeta
    AB = i.coordinates-j.coordinates
    AB2 = np.dot(AB,AB)
    P = (i.alpha*i.coordinates+j.alpha*j.coordinates)/zeta

    if((np.count_nonzero(i.l)==0)and(np.count_nonzero(j.l)==0)):
        return xi*(3-2*xi*AB2)*overlap(i,j)
    elif((np.count_nonzero(i.l)==1)and(np.count_nonzero(j.l)==0)):
        i_m = np.nonzero(i.l)
        si = primitive_gaussian(i.alpha,i.coeff,i.coordinates,0,0,0)
        return (P[i_m]-i.coordinates[i_m])*kinetic(si,j)+2*xi*overlap(i,j)
    elif((np.count_nonzero(i.l)==0)and(np.count_nonzero(j.l)==1)):
        j_m = np.nonzero(j.l)
        sj = primitive_gaussian(j.alpha,j.coeff,j.coordinates,0,0,0)
        return (P[j_m]-j.coordinates[j_m])*kinetic(i,sj)+2*xi*overlap(i,j)
    elif((np.count_nonzero(i.l)==1)and(np.count_nonzero(j.l)==1)):
        i_m = np.nonzero(i.l)
        j_m = np.nonzero(j.l)
        si = primitive_gaussian(i.alpha,i.coeff,i.coordinates,0,0,0)
        sj = primitive_gaussian(j.alpha,j.coeff,j.coordinates,0,0,0)
        if(i_m==j_m):
            return (P[j_m]-j.coordinates[j_m])*kinetic(i,sj)+kinetic(si,sj)/(2*zeta)+2*xi*overlap(i,j)
        else:
            return (P[j_m]-j.coordinates[j_m])*kinetic(i,sj)+2*xi*overlap(i,j)

def make_matrixes(molecule):
    
    nbasis = len(molecule)
    
    S = np.zeros([nbasis, nbasis])
    T = np.zeros([nbasis, nbasis])
    
    for i in range(nbasis):
        for j in range(nbasis):
            
            nprimitives_i = len(molecule[i])
            nprimitives_j = len(molecule[j])
            
            for k in range(nprimitives_i):
                for l in range(nprimitives_j):

                    S[i,j] += overlap(molecule[i][k],molecule[j][l])
                    T[i,j] += kinetic(molecule[i][k],molecule[j][l])

    return S , T

molecule_coordinate = [[0.0,0.0,0.116321],[0.0,0.751155,-0.465285],[0.0,-0.751155,-0.465285]]

# create H2O molecule - sto-3g basis
O_pg1sa  = primitive_gaussian(0.1307093214E+03, 0.1543289673E+00, molecule_coordinate[0], 0, 0, 0)
O_pg1sb  = primitive_gaussian(0.2380886605E+02, 0.5353281423E+00, molecule_coordinate[0], 0, 0, 0)
O_pg1sc  = primitive_gaussian(0.6443608313E+01, 0.4446345422E+00, molecule_coordinate[0], 0, 0, 0)

O_pg2sa  = primitive_gaussian(0.5033151319E+01,-0.9996722919E-01, molecule_coordinate[0], 0, 0, 0)
O_pg2sb  = primitive_gaussian(0.1169596125E+01, 0.3995128261E+00, molecule_coordinate[0], 0, 0, 0)
O_pg2sc  = primitive_gaussian(0.3803889600E+00, 0.7001154689E+00, molecule_coordinate[0], 0, 0, 0)

O_pg2pxa = primitive_gaussian(0.5033151319E+01, 0.1559162750E+00, molecule_coordinate[0], 1, 0, 0)
O_pg2pxb = primitive_gaussian(0.1169596125E+01, 0.6076837186E+00, molecule_coordinate[0], 1, 0, 0)
O_pg2pxc = primitive_gaussian(0.3803889600E+00, 0.3919573931E+00, molecule_coordinate[0], 1, 0, 0)

O_pg2pya = primitive_gaussian(0.5033151319E+01, 0.1559162750E+00, molecule_coordinate[0], 0, 1, 0)
O_pg2pyb = primitive_gaussian(0.1169596125E+01, 0.6076837186E+00, molecule_coordinate[0], 0, 1, 0)
O_pg2pyc = primitive_gaussian(0.3803889600E+00, 0.3919573931E+00, molecule_coordinate[0], 0, 1, 0)

O_pg2pza = primitive_gaussian(0.5033151319E+01, 0.1559162750E+00, molecule_coordinate[0], 0, 0, 1)
O_pg2pzb = primitive_gaussian(0.1169596125E+01, 0.6076837186E+00, molecule_coordinate[0], 0, 0, 1)
O_pg2pzc = primitive_gaussian(0.3803889600E+00, 0.3919573931E+00, molecule_coordinate[0], 0, 0, 1)

H1_pg1sa = primitive_gaussian(0.3425250914E+01, 0.1543289673E+00, molecule_coordinate[1], 0, 0, 0)
H1_pg1sb = primitive_gaussian(0.6239137298E+00, 0.5353281423E+00, molecule_coordinate[1], 0, 0, 0)
H1_pg1sc = primitive_gaussian(0.1688554040E+00, 0.4446345422E+00, molecule_coordinate[1], 0, 0, 0)

H2_pg1sa = primitive_gaussian(0.3425250914E+01, 0.1543289673E+00, molecule_coordinate[2], 0, 0, 0)
H2_pg1sb = primitive_gaussian(0.6239137298E+00, 0.5353281423E+00, molecule_coordinate[2], 0, 0, 0)
H2_pg1sc = primitive_gaussian(0.1688554040E+00, 0.4446345422E+00, molecule_coordinate[2], 0, 0, 0)

np.array(molecule_coordinate[2])]
O_1s  = [O_pg1sa, O_pg1sb, O_pg1sc]
O_2s  = [O_pg2sa, O_pg2sb, O_pg2sc]
O_2px = [O_pg2pxa, O_pg2pxb, O_pg2pxc]
O_2py = [O_pg2pya, O_pg2pyb, O_pg2pyc]
O_2pz = [O_pg2pza, O_pg2pzb, O_pg2pzc]
H1_1s = [H1_pg1sa, H1_pg1sb, H1_pg1sc]
H2_1s = [H2_pg1sa, H2_pg1sb, H2_pg1sc]
molecule = [O_1s, O_2s, O_2px, O_2py, O_2pz, H1_1s, H2_1s]

S , T = make_matrixes(molecule)
for i in range(7):
    for j in range(i+1):
        print("{:.6f}".format(S[i][j]),end="\t")
    print("\n")

for i in range(7):
    for j in range(i+1):
        print("{:.6f}".format(T[i][j]),end="\t")
    print("\n")

```

