% uses explicit expression for one- and two-electron integrals
% basis set slater exponents
% no schmidt orthogonalization
global Z
Z=2;
ze=[1.5;3];


res12=two_electron(ze);
res1=two_electron([ze(1) ze(1)]);
res2=two_electron([ze(2) ze(2)]);

syms r r1
psi=sqrt(4*ze.^3).*exp(-ze*r);
ppsi=psi;
%overlap matrix 
clear smat
smat=[1 res12(1);res12(1) 1];


% kinetic energy matrix
kemat=[res1(2) res12(2);res12(2) res2(2)];
% electron-nuclear attraction energy in schmidt basis
v1mat=[res1(3) res12(3);res12(3) res2(3)];
   



% determine two electron integrals
clear v2
clear v3
% two-electron integrals
% here for all indices the same
v3(1,1,1,1)=res1(4);
v3(2,2,2,2)=res2(4);
% here for i=j=k=1 and l=2
v3(1,2,1,1)=res12(7);
v3(2,1,1,1)=res12(7);
v3(1,1,1,2)=res12(7);
v3(1,1,2,1)=res12(7);
% here for i=j=k=2 and l=1
v3(2,2,1,2)=res12(6);
v3(2,2,2,1)=res12(6);
v3(2,1,2,2)=res12(6);
v3(1,2,2,2)=res12(6);
% here for i=j=2 and k=l=1
v3(2,2,1,1)=res12(4);
v3(1,1,2,2)=res12(4);
v3(1,2,1,2)=res12(5);
v3(2,1,1,2)=res12(5);
v3(2,1,2,1)=res12(5);
v3(1,2,2,1)=res12(5);



% hartree-fock iterations

c=[1 0];

clear v2
c=[1 0];
emin=0;
for iter=1:10
% determine 2x2 matrix of electron repulsion
for i=1:2;for j=1:2
      v2mat(i,j)=double(c*squeeze(v3(i,j,:,:))*c');
   end;end

% determine full matrix of hamiltonian
hmat=double(kemat+v1mat+v2mat);

[evec eval]=eig(hmat,smat);
c=evec(:,1)';
% eigenvectors are now ordered by row
ehf=eval(1)*2-c*v2mat*c';
if abs(ehf-emin)>1.e-4
    emin=ehf;
else
    [iter eval(1) c ehf]
%         [iter eval(1) c eval(1)*2-c*v2mat*c']

    break
end

end


 