%% SSP Projekt- 1s
%Lena Gasser, Michael Gobald                    20.06.2018
%
% In this script the roothan equations for a 1s atomic orbital will be
% solved.
% Therefor some function are used:
% see also: ao_1s.m     Heavysidestep.m     Laplace1s.m
%
% for more information about the used formulars check the page: 
% http://lampx.tugraz.at/~hadley/ss1/molecules/atoms/1s.php
%

close all
%% Set Values
%a0         Bohr radius in m
%hbar       reduced plack wirkungsquantum in Js
%ec         electron charge in J
%eps0       permetivity constant in As/Vm
%me         electron mass in kg
%Z_eff      effective nuclear charge
%delta      small value needed for the Heavysidestep function
%L & L_e    often used value

a0 = 5.2917721E-11; 
hbar = 1.05457266E-34;
ec = 1.60217733E-19; 
eps0 = 8.854187817E-12; 
me = 9.10938356*10^-31;
Z_eff = 1;
Za = 1;

delta = a0/(10*Z_eff);


L_e = 2*a0/Z_eff;
L = 8*a0/Z_eff;

%% Calculation of E1
% 
% The formular used to calculate the Energy comes from the Schrödinger
% equation
%
% n         numbers of iterations
% x,y,z     cartesian coordinates 

n = 10000;

E1 = zeros(n,1);

for i = 1:n
    x = L_e*(0.5-rand);
    y = L_e*(0.5-rand);
    z = L_e*(0.5-rand);
    xi = L_e*(0.5-rand);
    yi = L_e*(0.5-rand);
    zi = L_e*(0.5-rand);
    
    E1(i) = -hbar^2./(2*me)*Laplace1s(x,y,z,xi,yi,zi,Z_eff,a0)./ao_1s( x,y,z,xi,yi,zi,Z_eff,a0 )-Z_eff*ec^2./(4*pi*eps0*sqrt((x-xi)^2+(y-yi)^2+(z-zi)^2));
end


E1 = sum(E1)/length(E1) % Energie in Joule

%% Calculation of the S Elements
% The S-Elements of the Roothan equations built the overlap matrix
% Because of the normed wavefunction the Elements S_11 and S_22 are equal
% to 1
% The Element S_21 is the same as S_12

%% Define xi, yi, zi
% For a Hydogen molecule r_1 = -0.38 e_x Angstrom and r_2 = 0.38 e_x
% Angstrom
% e_x describes the unit vector in x direction, this causes, that only the
% x-elements have values unequal zero

x1 = -0.38*10^-10;      
x2 = 0.38*10^-10;
y1 = 0;
y2 = 0;
z1 = 0;
z2 = 0;

%Integration borders
xmin = -L/2;
xmax = -xmin;
ymin = xmin;
ymax = xmax;
zmin = xmin;
zmax = xmax;

%% Integration
% for the numeric intergration the matlab function 'integral3' is used.

integrand_S12 = @(x,y,z) ao_1s(x,y,z,x1,y1,z1,Z_eff,a0).*ao_1s(x,y,z,x2,y2,z2,Z_eff,a0);

S_12 = integral3(integrand_S12,xmin, xmax, ymin, ymax,zmin,zmax)


%% Calculation of the H Elements
% The H elements are builing the Hamiltonian matrix of the Roothan 
% equations 
%
% To claculate these values, some approximations were used (see the link
% above)

%% H11 Element
% 
% C_11      Constant part of the sum (comes from the approximations)

C_H11 = Z_eff^4*ec^2*delta^2/(2*pi*eps0*a0^3)*exp(-2*Z_eff*sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2)/a0);

%Integration
integrand_H11 = @(x,y,z) Heavysidestep( x, y, z, x2, y2, z2, delta).*ao_1s(x,y,z,x1,y1,z1,Z_eff,a0) ...
    .*ao_1s(x,y,z,x1,y1,z1,Z_eff,a0).*(Z_eff*ec.^2)./(4*pi*eps0*sqrt(((x-x2).^2)+((y-y2).^2)+((z-z2).^2)));

Integral_H11 = integral3(integrand_H11,xmin, xmax, ymin, ymax,zmin,zmax);


H11 = (E1-C_H11-Integral_H11)/ec

%% Monte Carlo H11
% The Monte Carlo Method is a used on the website to calculate the
% Integral.
%
% m         Number of times the Monte Carlo Methode is calculated
% N         Number fo iterations for the monte carlo

m = 100;
H11_m = zeros(1,m);
N = 1000000;
for j = 1:m
    for i = 1:N
        x_m = L*(0.5-rand());
        y_m = L*(0.5-rand());
        z_m = L*(0.5-rand());
    
        H11_m(j) = H11_m(j) - Heavysidestep(x_m,y_m,z_m,x2,0,0,delta)*ao_1s(x_m,y_m,z_m,x1,0,0,1,a0)*ao_1s(x_m,y_m,z_m,x1,0,0,1,a0)./sqrt(((x_m-x2).^2)+((y_m).^2)+((z_m).^2));
    end
    
    H11_m(j) = L^3*H11_m(j)*ec/(4*pi*eps0*N);
    H11_m(j) = H11_m(j)-C_H11/ec+E1/ec;
end

%% Differnce between Monte Carlo and exact Method
% The Monte Carlo Method is quite usefule, but it is not that exact. 
% To show that, the Methode is doneseveral times and each result is plotted. 
% The with the numerical integral calculated value and the mean value of 
% all Monte Carlo Calculations are also shown in the figure.
%
% H11_m_average         Average of the H11_m Values
% H11_plot              H11 Integral value put in a vector to plot them as
%                       a line

H11_m_average = zeros(1,length(H11_m))+sum(H11_m)/length(H11_m);

H11_plot = zeros(1,m)+H11;

fig_H11 = figure;
plot(1:m,H11_plot,'-k')
hold on
plot(1:m,H11_m,'xr')
plot(1:m,H11_m_average,'--b')
hold off
xlim ([1 m])
ylim auto
title({'Calculation of the H11-Element:';'Camparison of the exact integral and the Monte Carlo method'})
xlabel('Number of iteration')
ylabel('H11 / eV')
legend('Exact Integral','Monte Carlo Method','Average of the Monte Carlo values','Location','best')

%% H12 Element
%
% C_12      Constant part of the sum (comes from the approximations)

C_H12 = Z_eff^4*ec^2*delta^2/(2*pi*eps0*a0^3)*exp(-Z_eff*abs(x2-x1)/a0);


%Integral
integrand_H12 = @(x,y,z) Heavysidestep( x, y, z, x1, y1, z1, delta).*ao_1s(x,y,z,x1,y1,z1,Z_eff,a0) ...
    .*ao_1s(x,y,z,x2,y2,z2,Z_eff,a0).*(Z_eff*ec.^2)./(4*pi*eps0*sqrt(((x-x1).^2)+((y-y1).^2)+((z-z1).^2)));

Integral_H12 = integral3(integrand_H12,xmin, xmax, ymin, ymax,zmin,zmax);


H12 = (E1*S_12-C_H12-Integral_H12)/ec

%% Monte Carlo H12
% The Monte Carlo Method is a used on the website to calculate the
% Integral.
%
% m2         Number of times the Monte Carlo Methode is calculated
% N         Number fo iterations for the monte carlo

m2 = 100;
H12_m = zeros(1,m2);
N = 1000000;

for j = 1:m
    for i = 1:N
        x_m2 = L*(0.5-rand);
        y_m2 = L*(0.5-rand);
        z_m2 = L*(0.5-rand);
        H12_m(j) = H12_m(j) + Heavysidestep( x_m2, y_m2, z_m2, x1, y1, z1, delta).*ao_1s(x_m2,y_m2,z_m2,x1,y1,z1,Z_eff,a0) ...
        .*ao_1s(x_m2,y_m2,z_m2,x2,y2,z2,Z_eff,a0)./sqrt((x_m2-x1).^2+(y_m2-y1).^2+(z_m2-z1).^2);
    end

H12_m(j) = H12_m(j)*L^3*ec/(4*pi*eps0*N);
H12_m(j) = (E1*S_12-C_H12)/ec-H12_m(j);
end

%% Differnce between Monte Carlo and exact Method
% The Monte Carlo Method is quite usefule, but it is not that exact. 
% To show that, the Methode is doneseveral times and each result is plotted. 
% The with the numerical integral calculated value and the mean value of 
% all Monte Carlo Calculations are also shown in the figure.
%
% H12_m_average         Average of the H12_m Values
% H12_plot              H12 Integral value put in a vector to plot them as
%                       a line

H12_m_average = zeros(1,length(H12_m))+sum(H12_m)/length(H12_m);

H12_plot = zeros(1,m2)+H12;

fig_H12 = figure;
plot(1:m2,H12_plot,'-k')
hold on
plot(1:m2,H12_m,'xr')
plot(1:m2,H12_m_average,'--b')
hold off
xlim ([1 m2])
ylim auto
title({'Calculation of the H12-Element:';'Camparison of the exact integral and the Monte Carlo method'})
xlabel('Number of iteration')
ylabel('H12 / eV')
legend('Exact Integral','Monte Carlo Method','Average of the Monte Carlo values','Location','best')

%% Savings
% save the claculated values in a .mat file to use them later if needed 
% save the figures as well

save('H_1s_values','H11','H12','H11_m','H12_m')

saveas(fig_H11,'H11_figure.png')
saveas(fig_H12,'H12_figure.png')
