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.


2 Answers 2


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)

  • 1
    $\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$ Sep 11, 2020 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.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .