Taken from Phys5870 : Modern Computational Methods in Solids by Adrien E. Feiguin and was my homework assignment for Simulation of Material Science course. If she does it with Fortran, I’ll use Octave instead :p.

The Schrodinger equation for Helium ground state energy $(l = 0)$ using Hartree-Fock Method, which base is

$\Phi ( \vec r) = e^{- \alpha_p r^2}$

where $\alpha_1 = 0.297104, \alpha_2 = 1.236745, \alpha_3 = 5.749982$ and $\alpha_4 = 38.21667$. The parameter of Schrodinger will give equations as:

$\sum_{pq} (h_{pq}+ \sum_{rs} c_r c_s Q_{pqrs})c_q = E \sum_{pq} S_{pq} c_q$
$h_{pq} = \langle \Phi_p | - \frac{1}{2} \Delta^2-\frac{2}{r}|\rangle$
$Q_{pqrs} = \frac{2 \pi^5/2}{(\alpha_p+\alpha_q)(\alpha_r+\alpha_s)\sqrt{\alpha_p+\alpha_q+\alpha_r+\alpha_s}}$
$S_{pq} = \langle \Phi_p | \Phi_q \rangle$

The algorithm to compute Helium Ground State Energy is:

1. Compute all matrix S,T,A,H,Q
2. Gives random value to Cp, enter the value c as zeros matrix 4×1
3. Self-consistency is done until the previous value of c and the new value of c converges
4. Normalize c by dividing it with $\sqrt{c'*S*c}$
5. Compute Fock matrix, and enter the value of previous c to buffer variable
6. Find eigenvector from the smallest value of eigenvalue, and update the value of c
7. Compute the difference of the new and previous value of c, and if it has converged, stop the iteration
8. Compute the energy.

%%Find the parameter S, Q, dan H
alpha = [0.297104; 1.236745; 5.749982; 38.216677];
N = 4;

% Construct S T A H
for p=1:N
for q=1:N
S(p,q) = (pi/(alpha(p)+alpha(q)))^(3/2);
T(p,q) = 3*alpha(p)*alpha(q)*pi^(1.5)/(alpha(p)+alpha(q))^(5/2);
A(p,q) = -4*pi/(alpha(p)+alpha(q));
H(p,q) = T(p,q)+A(p,q);
end
end

% Construct Q
for p=1:N
for q=1:N
for r=1:N
for s=1:N
Q(p,q,r,s) = 2*pi^(2.5)/(((alpha(p)+alpha(q))*(alpha(r)+alpha(s))*(alpha(p)+alpha(q)+alpha(r)+alpha(s))^(0.5)));
end
end
end
end

%% Self-Consistency Hartree-Fock
c = zeros(4,1); c(1) = 1; % random cp, the value of q r s is set to null
c1 = zeros(4,1) % dummy for updating the value of c
F = zeros(4,4); v = zeros(4,4);
testing = 0;
const = 1; NormS= 0;

while(const >=1e-10)
NormS = c’*S*c;
% Normalized
c = c/((NormS)^(0.5));

testing = testing+1 % just testing…

% Create Fock Matrix F
for p=1:N
for q=1:N
F(p,q) = H(p,q);
for s=1:N
for r=1:N
F(p,q) = F(p,q) + c(r)*c(s)*Q(p,q,r,s);
end
end
end
end

% find the eigenvector from the smallest eigenvalue
c1 = c;    %update eigenvector
% F*V = S*V*D % V as cv, D as E
[V, D] = eig(F,S);
c = V(:,1)

% converged? y: stop; n: loop
const = (c-c1)’*(c-c1);
end

E = 0;
% Find the Energy
for p=1:N
for q=1:N
E = E + 2*c(p)*c(q)*H(p,q)
for r=1:N
for s=1:N
E = E + Q(p,q,r,s)*c(p)*c(q)*c(r)*c(s)
end
end
end
end

If you compute it right, you will get the Helium Energy as $E = -2.8552 Hartree$.