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

where and . The parameter of Schrodinger will give equations as:

The algorithm to compute Helium Ground State Energy is:

- Compute all matrix S,T,A,H,Q
- Gives random value to Cp, enter the value c as zeros matrix 4×1
- Self-consistency is done until the previous value of c and the new value of c converges
- Normalize c by dividing it with
- Compute Fock matrix, and enter the value of previous c to buffer variable
- Find eigenvector from the smallest value of eigenvalue, and update the value of c
- Compute the difference of the new and previous value of c, and if it has converged, stop the iteration
- 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 .

## Leave a Reply