Skip to main content
Matter Modeling

Return to Answer

added 8 characters in body
Source Link
Nike Dattani - No Free Time
Nike Dattani - No Free Time
  • 36.7k
  • 4
  • 100
  • 252

  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: $$N=\bigg(\frac{2\zeta}{\pi}\bigg)^{3/4}$$$$\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. 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 $$N=\bigg(\frac{2\zeta}{\pi}\bigg)^{3/4}$$ The normalization will be different for higher angular momentum functions (p, d, f, etc).

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

Source Link
Tyberius ♦
Tyberius ♦
  • 15.9k
  • 4
  • 31
  • 122

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 $$N=\bigg(\frac{2\zeta}{\pi}\bigg)^{3/4}$$ The normalization will be different for higher angular momentum functions (p, d, f, etc).