I have a set of basis functions defined as:
def HO_wavefunction(x, n, x0, omega, m=1):
N = 1.0 / math.sqrt(2**n * math.factorial(n)) * ((m * omega)/(math.pi))**(0.25) # Normaliziation constant
y = (np.sqrt(m * omega)) * (x - x0)
return N * np.exp(-y * y / 2.0) * sp.hermite(n)(y)
#Define the basis
def enol_basis(x, n):
return HO_wavefunction(x, n, x0=Enolminx, omega=wenol)
I now want to compute the overlap integrals Sii = integral((SiSi)dx), Sjj = integral((SjSj)dx) and Sij = integral((Si*Sj)dx) of my basis functions and store them in some type of array. I tried the following:
G = 10
S = np.empty([G,G])
for n in range (G-1):
for m in range (G-1):
S[n][m]= np.trapz(enol_basis(x,n)*enol_basis(x,m),x)
print (S[n][m])
This only returns a single value instead of all the results stored in an array. If anyone could help me compute the overlap integrals as I defined them above and store the results in an array I would really appreciate it!