As per @user1271772 's suggestion I am asking this question here again.

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

Do an SCF calculation for the helium-atom ground state using basis set of two 1s STOs with exponent $\zeta_1$ = 1.45 and $\zeta_2$=2.91. [Ref: Roetti and Clementi J. Chem. Phys., 60, 4725 1974] $\chi_1 = 2\zeta_1^{3/2}e^{-\zeta_1 r}Y^0_0$, and $\chi_2 = 2\zeta_2^{3/2}e^{-\zeta_2 r}Y^0_0$

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

clear all;
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)));

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 Levine book as :

$$ (rs|tu) = \int \int \frac{\chi^*_r(1)\chi_s(1)\chi_t^*(2)\chi_u(2)}{r_{12}} dv_1dv_2 $$

user @TAR86 from the Chemistry SE suggested that

You replaced the 6-fold integration by one in spherical coordinates. Not sure if that can work as easily as you wrote it

But I was thinking as the function have no angular dependence, at least in above case, its integral should be really straight forward.


Your Octave code is trying to do the integral by quadrature, which makes very little sense since it will have a huge problems with the cusp.

Since this is a one-center problem, the best approach is to use the Legendre expansion for $|r_1-r_2|^{-1}$, which decomposes the interaction into a radial part and an angular part: $r_{12}^{-1} = \frac {4\pi} {r_>} \sum_{L=0}^\infty \frac 1 {2L+1} \left( \frac {r_<} {r_>} \right)^L \sum_{M=-L}^L Y_L^M (\Omega_1) (Y_L^M (\Omega_2))^*$.

You only have $s$ orbitals; this means that the angular parts are trivial and only a single term drops out; you're left with the radial integral $\int_0^\infty {\rm d}r r^2 \int_0^\infty {\rm d}r' {r'}^2 \chi_r(r) \chi_s(r) \chi_t(r') \chi_u(r') r_>^{-1} $ that you can solve by standard techniques i.e. dividing the integration in two parts for $r'\leq r$ and $r'>r$ and then evaluating these integrals separately.

This trick is also what makes fully numerical calculations possible on atoms, see e.g. my recent review in Int J Quantum Chem 119, 19, e25968 (arXiv:1902.01431) and the application to finite-element calculations on atoms in Int J Quantum Chem 119, 19, e25945 (arXiv:1810.11651)

| cite | improve this answer | |
  • $\begingroup$ I recommend using display mathematics whenever the formulae get as complicated as here. Invoke it with $$ ... $$ or \begin{align|equation} ... \end{...} or if really necessary inline with $\displaymath ... $. I'd have edited it, but these are too few characters. $\endgroup$ – Martin - マーチン Sep 11 at 17:00

This can be solved analytically, a complete solution can be found here

To refrain from rewriting the entire derivation I will only say that you need to integrate over all 3 dimensional degrees of freedom for both electrons, so TAR86 is correct.

In the derivation at the link, the distance between the electrons ($\mid r_1 - r_2 \mid \equiv r_{12}$) is better represented in polar coordinates (equation 1196).

Going through the algebra and calculus one can then end up with the correct solution for the Coulomb integral, -5/2 E0.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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