2
$\begingroup$

I'd want to implement SCF method using Modern Quantum Chemistry (Szabo). But when I want to calculate overlap matrix for H2 (STO-3G basis) my results aren't equal to true overlap matrix that there's in the book. However when I compare my results in other basis (STO-1G) I get exact overlap matrix. Below my code for calculation overlap values.

class gauss():

  def __init__(self, alpha, Ra, k = 1):

    self.alpha = alpha
    self.Ra = np.array(Ra)
    self.k = k

  def __mul__(self, other):
    'calculate integral'
    lenght = (self.Ra - other.Ra)
    scalar = np.dot(lenght, lenght)

    return np.exp(-self.alpha*other.alpha*scalar/(self.alpha + other.alpha))*np.pi/(self.alpha + other.alpha)**(3/2)



def overlap(first, second):

  s11 = 0
  for i in first:
    for j in first:

        s11 += i * j * i.k * j.k

  s12 = 0
  for i in first:
    for j in second:

        s12 += i * j * i.k * j.k

  return s11, s12

## gaussian primitives:
s1 = gauss(0.168855404, (0.0, 0.0, 0.0), 0.444634542)
s2 = gauss(0.623913730, (0.0, 0.0, 0.0), 0.535328142) 
s3 = gauss(3.42525091, (0.0, 0.0, 0.0), 0.154328967)   
s1b = gauss(0.168855404, (1.4, 0.0, 0.0), 0.444634542)
s2b = gauss(0.623913730, (1.4, 0.0, 0.0), 0.535328142)
s3b = gauss(3.42525091, (1.4, 0.0, 0.0), 0.154328967)



first = [s1, s2, s3]
second = [s1b, s2b, s3b]

s11, s12 = overlap(first, second)

print(s12 / s11) # print normalized s12

0.781180936263708

but exact value from book = 0.6593 (1.4 Bohr)

and if I choose basis as:

s = gauss(0.4166127, (0.0, 0.0, 0.0))
sb = gauss(0.4166127, (1.4, 0.0, 0.0))

I get exact value:

0.6647924143035879

from book: 0.6648

$\endgroup$

closed as off-topic by A.K., Todd Minehardt, Nuclear Chemist, Soumik Das, tschoppi Feb 19 at 9:10

  • This question does not appear to be about chemistry within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 1
    $\begingroup$ I have down-voted this question, because it has little to nothing to do with chemistry. I believe this is a programming issue, since all that is given here is code. It is not even explained what this code actually is supposed to do; it is not even stated in which language it is written. I am surprised that no-one has voted to close it as off-topic. $\endgroup$ – Martin - マーチン Jan 21 at 14:52
  • $\begingroup$ @Martin-マーチン after looking the question over more, I agree. While the issue the OP had wound up being chemistry related (misunderstanding about how to normalize overlap integrals), the question as phrased is more about debugging this particular code than addressing the underlying math. The OP does say they are trying to develop an SCF code, but not specifying the language doesn't help. I think a version of this question exists that would be on-topic, but as is I agree that it should be closed for now. $\endgroup$ – Tyberius Jan 21 at 19:43
  • 6
    $\begingroup$ I'm voting to close this question as off-topic because it is about code debugging and not an underlying chemical principle. $\endgroup$ – Tyberius Jan 21 at 19:44
  • $\begingroup$ I'm voting to close this question as off-topic because stack exchange's chemistry area should be devoted to the discussion of chemistry and not how to write code for a bit of computational chemistry software. To my mind it is like going on a diesel engine design forum and asking how to change the oil filter on a specific car when the forum is devoted to the theroy of how to design an engine. $\endgroup$ – Nuclear Chemist Feb 10 at 15:24
5
$\begingroup$

You have two errors in your code, one minor and one major.

The minor error is in your multiplication function. The value you return should have $\pi^{3/2}$, but the way it is arranged now the power is only being applied to the denominator.

The major error is how you are normalizing your overlap integrals. You can't simply divide $\langle1|2\rangle$ by $\langle1|1\rangle$ in order to normalize it. You have to explicitly normalize each of the Gaussian basis functions. To do this you can include a factor of ((4*self.alpha*other.alpha)/np.pi**2)**(3/4) in the value you are returning from your multiplication (in terms readability of your code, it may be wiser to define variables within the function like proportionality_constant, gauss_integral, and normalization and then return the product of these). Or you could simplify and just multiply the exponential by ((4*self.alpha*other.alpha)/(self.alpha+other.alpha)**2)**(3/4) This also means you only need the second part of your overlap function, since the multiplication will take care of the normalization.

The reason your code works for STO-1G is that dividing by $\langle1|1\rangle$ accidentally gives the correct normalization since there is only one normalization constant associated with the one Gaussian function of STO-1G.

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.