```function [ H,S] = calcHamiltonian( orb,Z,xAtom,yAtom,zAtom )
%CALCHAMILTONIAN calculates Hamiltonian and overlap matrix
% Uses molecular orbitals and doesn't take nuclei-nuclei interaction into
% account.
%
% Z ... a [m x 1] matrix containing the atomic numbers of all atoms in the
%       molecule. Atoms with n orbitals have to mentioned n times
% xAtom, yAtom and zAtom are [m x 1] matrizes containing relative positions
%       of all molecules (in Meter).
% orb ... {m x 1} cell containing {1 x n} cells containing strings
%       describing the n orbitals which
%       should be used for calculation (1s, 2s, 2px, 2py, 2pz ...).
%
%H is a matrix containing energy expectation values, S is the overal matrix

a0 = 5.2917721E-11; % Bohr radius in meters
%
%  x_lim = a0*[-300 300];
%  y_lim = a0*[-300 300];
%  z_lim = a0*[-300 300];

universallimit = 50*a0;

x_lim_l = -universallimit;
x_lim_u = universallimit;

y_lim_l = -universallimit;
y_lim_u = universallimit;

z_lim_l = -universallimit;
z_lim_u = universallimit;

hbar = 1.05457266E-34;
me = 9.1093897E-31;
eps0 = 8.854187817E-12;
ec = 1.60217733E-19;
k1 = 1/2*hbar^2/me;
k2 = ec^2/(4*pi*eps0);

%---------------------Set tolerance of integrals---------------------------
ATol = ec/1E2;
RTol = 1E-5;

num_orbs_per_atom = cellfun('size',orb,2);
nopa_cumsum = cumsum(num_orbs_per_atom);
N_atom = length(Z);
N = sum(num_orbs_per_atom);
orblist = {};
for k = 1:length(orb)
orblist = [orblist orb{k}];
end

H = NaN(N);
S = H;

for k = 1:N
for l = k:N
atom_k = find(k >= nopa_cumsum,1,'last');
atom_l = find(l >= nopa_cumsum,1,'last');

if(k==l) S(k,l) = 1;
else
S(k,l) = integral3(@(x,y,z) atorbit(x-xAtom(atom_k),y- yAtom(atom_k), z-zAtom(atom_k),Z(atom_k),false,orblist{k}).* ...
atorbit(x-xAtom(atom_l),y- yAtom(atom_l),z- zAtom(atom_l),Z(atom_l),false,orblist{l}), ...
x_lim_l,x_lim_u,y_lim_l,y_lim_u,z_lim_l,z_lim_u,'AbsTol',ATol,'RelTol',RTol);
end
fprintf('Done S(%i,%i)\n',k,l);
end
end
S(tril(true(N),-1)) = S(triu(true(N),1));

for k = 1:N
for l = 1:N
atom_k = find(k >= nopa_cumsum,1,'last');
atom_l = find(l >= nopa_cumsum,1,'last');

orb1 = orblist{k};
orb2 = orblist{l};

pos1 = [xAtom(atom_k) yAtom(atom_k) zAtom(atom_k)];
pos2 =  [xAtom(atom_l) yAtom(atom_l) zAtom(atom_l)];

%n1 = str2double(orb1(1));
n2 = str2double(orb2(1));

%eigenE1 = -me*ec^4/(2*(4*pi*eps0)^2*hbar^2)/n1^2*Z1^2;
eigenE2 = -me*ec^4/(2*(4*pi*eps0)^2*hbar^2)/n2^2*Z(atom_l)^2;

H(k,l) = eigenE2*S(k,l);

% Exclude atom corresponding to right eigenfunction (atom_l)
for j = [1:atom_l-1 atom_l+1:N_atom]
H(k,l) = H(k,l) - k2*Z(j)*...
integral3(@(x,y,z) atorbit(x-pos1(1),y-pos1(2),z-pos1(3),Z(atom_k),false,orb1)...
./sqrt((x-xAtom(j)).^2+(y-yAtom(j)).^2+(z-zAtom(j)).^2).* ...
atorbit(x-pos2(1),y-pos2(2),z-pos2(3),Z(atom_l),false,orb2),...
x_lim_l,x_lim_u,y_lim_l,y_lim_u,z_lim_l,z_lim_u,'AbsTol',ATol,'RelTol',RTol);
fprintf('Done H(%i,%i) - Atom %i\n',k,l,j);
end

%
% mask = 1:N_atom ~= atom_l;
% H(k,l) = H(k,l) - k2*...
%         integral3(@(x,y,z) expecval(x,y,z,xAtom(mask),yAtom(mask),zAtom(mask),Z(mask),orb1,orb2,pos1,pos2,Z(atom_k),Z(atom_l)),...
%         x_lim_l,x_lim_u,y_lim_l,y_lim_u,z_lim_l,z_lim_u,'AbsTol',1e-20,'RelTol',1e-5);
end
end

% integral3(@(x,y,z) atorbit(x,y,z,0,0,0,1,false,'s2'),-1,1,-1,1,-1,1);

end

function H = expecval(x,y,z,xAtom,yAtom,zAtom,Z,orb1,orb2,pos1,pos2,Z1,Z2)
%atom at pos2 has to be excluded from atommatrix xAtom,yAtom etc.
%
%   hbar = 1.05457266E-34;
%  me = 9.1093897E-31;
%  eps0 = 8.854187817E-12;
%  ec = 1.60217733E-19;
%  %k1 = 1/2*hbar^2/me;
%  k2 = ec^2/(4*pi*eps0);
%
%  %n1 = str2double(orb1(1));
%  n2 = str2double(orb2(1));
%
%  %eigenE1 = -me*ec^4/(2*(4*pi*eps0)^2*hbar^2)/n1^2*Z1^2;
%  eigenE2 = -me*ec^4/(2*(4*pi*eps0)^2*hbar^2)/n2^2*Z2^2;
%
% H = eigenE2*atorbit(x-pos1(1),y-pos1(2),z-pos1(3),Z1,false,orb1).*atorbit(x-pos2(1),y-pos2(2),z-pos2(3),Z2,false,orb2);

H=0;
for j = 1:length(Z)
H = H - atorbit(x-pos1(1),y-pos1(2),z-pos1(3),Z1,false,orb1).* ...
Z(j)./sqrt((x-xAtom(j)).^2+(y-yAtom(j)).^2+(z-zAtom(j)).^2).* ...
atorbit(x-pos2(1),y-pos2(2),z-pos2(3),Z2,false,orb2);
end

end
```