%============================================================
% Tony Hyun Kim
% 8.06: Hartree-Fock implementation on two hydrogen molecules
%       using four 1s-STO functions as the basis set.
%
% All files can be obtained at: 
%   http://web.mit.edu/kimt/www/8.06/paper/
%
% General notes:
%   - This is the second HF implementation I worked on for 
%       the 8.06 term paper. This has been written in a way that
%       corresponds exactly to the discussion given in the paper.
%   
% Dependencies:
%   sto1s.m: Defines the 1s Slater-type wavefunction
%============================================================

%============================================================
% Program parameters
%============================================================

% Grid parameters
x_lower = -1; x_upper = 3.5; x_steps = 10;
y_lower = -1; y_upper = 1;   y_steps = 10;
z_lower = -1; z_upper = 1;   z_steps = 10;

% To make it easy to pass the above nine parameters...
grid_params = [x_lower x_upper x_steps;...
               y_lower y_upper y_steps;...
               z_lower z_upper z_steps];
clear x_lower x_upper x_steps;
clear y_lower y_upper y_steps;
clear z_lower z_upper z_steps;
           
%============================================================
% Step 1: Enter the atomic coordinates and the basis functions
%============================================================
R = [0   0 0;... % Hydrogen 1
     0.5 0 0;... % Hydrogen 2
     2   0 0;... % Hydrogen 3
     2.5 0 0];   % Hydrogen 4
 
Z = [1 1 1 1]; % The corresponding charges

V_nucl_fun = @(r)nucl_potential(r,Z,R); % Corresponding nuclear
                                        % potential function
                                    
% The basis functions are 1s-STO functions centered at 
% each of the nuclei
phi1 = @(r)sto1s(r,R(1,:),1.24);
phi2 = @(r)sto1s(r,R(2,:),1.24);
phi3 = @(r)sto1s(r,R(3,:),1.24);
phi4 = @(r)sto1s(r,R(4,:),1.24);
 
phi = {phi1 phi2 phi3 phi4};
clear phi1 phi2 phi3 phi4;

%============================================================
% STEP 2:
%============================================================
% The following matrices need to be computed only ONCE per a
% set of iterations.
%
% Compute S:  the overlap matrix,
%         T:  the kinetic energy matrix,
%         V_nucl: Nuclear potential
%         H_core:  core-Hamiltonian, H = T + V_nucl
%
%         I: The two-electron integrals
%============================================================

%------------------------------------------------------------
% Compute the overlap matrix, S
%------------------------------------------------------------
S = zeros(4,4);
disp('Computing the overlap matrix S...');
tic; % For benchmarking purposes. We can measure the amount
     % of time taken for Matlab to perform the computations
for i = 1:4
    for j = i:4
        S(i,j) = inner_product(phi{i},phi{j},grid_params);
                        
        % Use the hermiticity of the overlap matrix
        S(j,i) = S(i,j);
    end
end
t_overlap = toc;
msg = sprintf('  Done! (Time elapsed: %-.4f seconds)', t_overlap);
disp(msg);

%------------------------------------------------------------
% Compute the kinetic energy matrix, T
%------------------------------------------------------------
T = zeros(4,4);
disp('Computing the kinetic energy matrix T...');
tic;
for i = 1:4
    for j = i:4
        T(i,j) = -1/2*inner_product(phi{i},...
                                    laplacian(phi{j},grid_params),...
                                    grid_params);
        
        T(j,i) = T(i,j);
    end
end
t_kinetic = toc;
msg = sprintf('  Done! (Time elapsed: %-.4f seconds)', t_kinetic);
disp(msg);

%------------------------------------------------------------
% Compute the nuclear potential matrix, V_nucl
%------------------------------------------------------------
V_nucl = zeros(4,4);
disp('Computing the nuclear potential matrix V_nucl...');
tic;
for i = 1:4
    for j = i:4
        V_nucl(i,j) = inner_product(phi{i},...
                                    @(r)(V_nucl_fun(r)*phi{j}(r)),...
                                    grid_params);
        V_nucl(j,i) = V_nucl(i,j);
    end
end
t_nucl_potential = toc;
msg = sprintf('  Done! (Time elapsed: %-.4f seconds)', t_nucl_potential);
disp(msg);

%------------------------------------------------------------
% Compute the core-Hamiltonian
%------------------------------------------------------------
disp('Core-Hamiltonian (H_core = T + V_nucl) obtained!');
H_core = T + V_nucl;

%------------------------------------------------------------
% Compute the two-electron integrals, I
%------------------------------------------------------------
I = zeros(4,4,4,4);
disp('Calculating the two-electron integrals I... (Sit tight!)');
tic;
for i = 1:4
    for j = i:4
        for k = j:4
            for l = k:4
                msg = sprintf('I(%d,%d,%d,%d)',i,j,k,l);
                disp(msg);
                I(i,j,k,l) = two_elec_integral(phi{i},phi{j},...
                                               phi{k},phi{l},...
                                               grid_params);
                                           
                % We apply the symmetry of I
                I(i,l,k,j) = I(i,j,k,l);
                I(j,i,l,k) = I(i,j,k,l);
                I(j,k,l,i) = I(i,j,k,l);
                I(k,j,i,l) = I(i,j,k,l);
                I(k,l,i,j) = I(i,j,k,l);
                I(l,i,j,k) = I(i,j,k,l);
                I(l,k,j,i) = I(i,j,k,l);
            end
        end
    end
end
toc;
t_two_elec_int = toc;
msg = sprintf('  Done! (Time elapsed: %-.4f seconds)', t_two_elec_int);
disp(msg);