```function [ H,S] = calcHamiltonian( orb,Z,xAtom,yAtom,zAtom,universallimit)
%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
%

constants

if nargin<6
universallimit = 5.2917721E-11*50; %a0*50
end

x_lim_l = -universallimit;
x_lim_u = universallimit;

y_lim_l = -universallimit;
y_lim_u = universallimit;

z_lim_l = -universallimit;
z_lim_u = universallimit;

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

num_orbs_per_atom = cellfun('size',orb,2);
nopa_cumsum = cumsum([0; 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): %i-%s, %i-%s\n', k, l, atom_k, orblist{k}, atom_l, orblist{l});
end
end

S = triu(S) + triu(S,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
end
end

end

```