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