% huckel for butadiene
% number of C atoms
n=4;
% number of pi orbitals
norb=n/2;
hmat=zeros(n,n);
for i=1:n-1;hmat(i,i+1)=1;end
hmat=hmat+triu(hmat)';
hmat=-hmat;
[evec eval]=eig(hmat);
e=diag(eval);
e=e(1:norb);

% resonance energy
er=(2*sum(e)+n)/n;

% to determine bond order

bo=sum((evec(1:n-1,1:norb).*evec(2:n,1:norb))');
bo=2*bo(1:norb);
disp('butadiene');
disp(['energies:  ' num2str(e')])
disp(['resonance energy: ' num2str(er)])
disp(['bond orders: ' num2str(bo)])

% to plot the orbital

% make sure coefficient of first orbital is positive
if evec(1,1)<0
   % if not, change the sign of the eigenvector
   evec(1,:)=-evec(1,:);
end
if evec(2,1)<0
   evec(2,:)=-evec(2,:);
end
plot([1:4],evec(:,1:2),'-o');set(gca,'xtick',[1:4],'xlim',[0.5 4.5],'fontsize',14)
hold
plot([-.5 4.5],[0 0],'-k')
hold
xlabel('Atom index i');ylabel('C_{ki}')