0
$\begingroup$

I am trying to understand SCF cycle by trying to code up solved example from Quantum Chemistry by Levine (page 443, 5th edition). The problem is shown below:

scf helium problem statement

One electron integrals are straight forward and I was able to get the correct answer, however I still can get correct values for two election integral, lets say (11|11). Given below is my attempt in octave

clear all;
N=2000;
zeta1 = 1.45;
zeta2 = 2.91;

r = linspace(0.000001,10,N)';
dr = r(2)-r(1);
chi = @(zetad,x) (2*zetad.^(3/2))*exp(-zetad*x).*x;
chichi = 0;
for i =1:N
    chichi = chichi + dr*(chi(zeta1,r(i))*chi(zeta1,r(i))*chi(zeta1,r')*(chi(zeta1,r)./((r(i)-r) + 0.000001)));
end
chichi*dr

However my values are way off in this case. Can anyone please shed a light on it? Value of (11|11) = 5/8 zeta1 = 0.9062.

Two electron integrals are defined in the book as:

two electron integrals

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ You replaced the 6-fold integration by one in spherical coordinates. Not sure if that can work as easily as you wrote it and I'm not an expert on that transformation, but I think that some multiplicative constants from the angle integrations are missing. $\endgroup$
    – TAR86
    May 29 '20 at 17:43
  • 1
    $\begingroup$ I am closing this question as it has been asked and answered on Matter Modeling: One-center two-electron integrals between 1s STO $\endgroup$
    – Martin - マーチン ♦
    Sep 11 '20 at 17:03
0
$\begingroup$

@TAR86 was correct, I mistook the volume integral as simple 1d integral. Also integrating $1/|r_2 - r_1|$ would result in a singularity at $r_1 == r_2$ which would yield really wrong results.

As per @user1271772 suggestions, materials modelling stack exchange was quite useful (see question here). There as per Susi Lehota's suggestion I was able to get correct result by integrating Legendre Expansion.

My octave code for two electron integral (11|11) follows:

zeta1=1.45;
chi1111 = @(r1,r2) 16*zeta1^3*zeta1^3*exp(-2*zeta1.*(r1)).*exp(-2*zeta1.*(r2)).*r1.*r1.*r2.*r2./max(r1,r2);
dblquad(chi1111,0.000001,10,0.0000001,10,1E-3)
ans =  0.90623

Similarly other integrals can be encoded, and I get results exactly as Levine example.

Thank you everyone

Improve this answer
$\endgroup$

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