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.

Advertisements