How to compute the overlap matrix in Python

Question

I want to calculate the overlap integral (S), I made the code, but it is only worked when I use a basis function that describe 1s and when I include 2s orbital, I start to get error. The First Problem is solved as Tyberius said in the comments, But I still got wrong overlap integral. I noticed that the error is because that part of my equation is not like what I should get. In details, Here what I have for example in zeta.txt.

#it is the zeta from the basis function description
    18.7311370000
    2.8253937000
    0.6401217000
    0.1612778000
    18.7311370000
    2.8253937000
    0.6401217000
    0.1612778000

My code below when I print(float(zetas[k])), I got most of the number correct but some of them are wrong, so what should be the error ?

   nb = 4
    S = np.zeros([nb,nb]) # nb is the number of basis set and I have it in the code 
    print(distancess) # where distancess is an array that contains the coordinate 
    # here is the shape of distancess =  [[0.0, 0.0, 0.0], [0.0, 0.0, 1.3999091131307069]]
    print(nprim) #is an array that contains the number of primitive for each basis function 
    # nprim = [3, 1, 3, 1]
    with open("zeta.txt","r") as z_inp:
        with open("cijk.txt","r") as cijk_inp:
            zetas = z_inp.readlines()
            cijk = cijk_inp.readlines()
            for i in range(nb):
                for j in range(nb):
    
                    n_prim_i = nprim[i]
                    n_prim_j = nprim[j]
    
                    for k in range(n_prim_i):
                        for l in range(n_prim_j):
                           Q = distancess[atom_for_basis[i]] - 
                              distancess[atom_for_basis[j]]
                           Q2 = np.dot(Q,Q)

                           zetazeta = float(zetas[k]) + float(zetas[l])
                           print(float(zetas[k]))
                
                           zi = float(zetas[k]) * 
                                float(zetas[l]) / zetazeta
                           c1c2 = float(cijk[k]) * 
                           float(cijk[l])
                           gaussint = ((math.pi/zetazeta))**(3/2)

                           prefac = math.exp(-zi*Q2)
                           S[i,j] +=  c1c2 * prefac * gaussint
    
    print(S)

The output from my code for print(float(zetas[k])) is,

18.731137
18.731137
18.731137
2.8253937
2.8253937
2.8253937
0.6401217
0.6401217
0.6401217
18.731137
2.8253937
0.6401217
18.731137
18.731137
18.731137
2.8253937
2.8253937
2.8253937
0.6401217
0.6401217
0.6401217
18.731137
2.8253937
0.6401217
18.731137
18.731137
18.731137
18.731137
18.731137
18.731137
18.731137
18.731137
18.731137
18.731137
18.731137
2.8253937
2.8253937
2.8253937
0.6401217
0.6401217
0.6401217
18.731137
2.8253937
0.6401217
18.731137
18.731137
18.731137
2.8253937
2.8253937
2.8253937
0.6401217
0.6401217
0.6401217
18.731137
2.8253937
0.6401217
18.731137
18.731137
18.731137
18.731137
18.731137
18.731137
18.731137
18.731137

The output that should I get for print(float(zetas[k])).

18.731137
18.731137
18.731137
2.8253937
2.8253937
2.8253937
0.6401217
0.6401217
0.6401217
18.731137
2.8253937
0.6401217
18.731137
18.731137
18.731137
2.8253937
2.8253937
2.8253937
0.6401217
0.6401217
0.6401217
18.731137
2.8253937
0.6401217
0.1612778
0.1612778
0.1612778
0.1612778
0.1612778
0.1612778
0.1612778
0.1612778
18.731137
18.731137
18.731137
2.8253937
2.8253937
2.8253937
0.6401217
0.6401217
0.6401217
18.731137
2.8253937
0.6401217
18.731137
18.731137
18.731137
2.8253937
2.8253937
2.8253937
0.6401217
0.6401217
0.6401217
18.731137
2.8253937
0.6401217
0.1612778
0.1612778
0.1612778
0.1612778
0.1612778
0.1612778
0.1612778
0.1612778

Answer 1

Here is an updated version of your code that I believe now works. I had to add in a few of the values that were computed elsewhere in your code, but otherwise it is mostly the same.

import numpy as np

nb = 4
S = np.zeros([nb,nb]) # nb is the number of basis set and I have it in the code
distancess=np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.4]])
#is an array that contains the number of primitive for each basis function
nprim = [3, 1, 3, 1]
atom_for_basis=[0,0,1,1]

zetas = 2*[18.7311370000,2.8253937000,0.6401217000,0.1612778000]
cijk = 2*[0.0334946434,.2347269535,0.8137573261,1.00000000]
norms= 2*[6.4170171 , 1.55317145, 0.51004325,0.18138065]

current_i=0
for i in range(nb):
    n_prim_i = nprim[i]
    current_j=0
    for j in range(nb):
        n_prim_j = nprim[j]

        for k in range(current_i,current_i+n_prim_i):
            for l in range(current_j,current_j+n_prim_j):
               Q = distancess[atom_for_basis[i]] - distancess[atom_for_basis[j]]
               Q2 = np.dot(Q,Q)

               zetazeta = float(zetas[k]) + float(zetas[l])

               zi = float(zetas[k]) * float(zetas[l]) / zetazeta
               c1c2 = float(cijk[k]) * float(cijk[l])
               gaussint = ((np.pi/zetazeta))**(3/2)


               prefac = np.exp(-zi*Q2)
               S[i,j] +=  c1c2 * prefac * gaussint*norms[k]*norms[l]
        current_j+=n_prim_j
    current_i+=n_prim_i
print(S)

So what still needed fixing?

1. Your loops over k and l always started at zero for every iteration, meaning you never got past the first three basis functions when accessing zeta and cijk. I added current_i and current_j to keep track of how far you had already travelled through these arrays. Having the loops include these variables and incrementing at the end of each i and j loop, it will now go through all of zeta and cijk.

2. You need to normalize the integrals or each of the individual basis functions (I discuss this for a similar question on Chem SE). Here, I explicitly included values to normalize each function, but you will need to add code to do this in general. For an s-type function, the normalization formula is: $$\tag{1}N=\bigg(\frac{2\zeta}{\pi}\bigg)^{3/4}$$ The normalization will be different for higher angular momentum functions (p, d, f, etc).