# How to compute the overlap matrix in Python

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:
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

• It was getting a bit long here, I moved the comments to a chat room
– Tyberius
Commented Dec 19, 2021 at 16:23
• Ok, I write there . Commented Dec 19, 2021 at 16:25
• I gave my +1 here a long time back, but a user has strongly criticized questions like yours here: mattermodeling.meta.stackexchange.com/q/312/5 You may want to take a look, especially since they selected this exact question as an example. Commented Jan 6, 2022 at 2:49

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).

• +1 and thanks for updating the user's code, then being so thorough in your explanations! Commented Dec 20, 2021 at 6:31
• First, thank you for your help. 1- I realized that my code run over zeta, but instead of moving from s1 basis set to s2 basis, my code repeat the s1 basis so this cause a problem. Now, the track you did, the code worked well 2- About the normalization, in general, yes I need to normalize, but I have the basis function normalized, so I didn't do it in my code. Commented Dec 20, 2021 at 9:56