% huckel for butadiene
% number of C atoms
n=4;
% number of pi orbitals
norb=n/2;

%determine the Hamiltonian matrix
hmat=zeros(n,n);
% set the upper triangle of the matrix
for i=1:n-1
   hmat(i,i+1)=1;
end
% symmetrize the matrix
hmat=hmat+triu(hmat)';
% change the sign
hmat=-hmat;

% diagonalize 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);

% print out the results

disp('butadiene');
disp(['energies:  ' num2str(e')])
disp(['resonance energy: ' num2str(er)])
disp(['bond orders: ' num2str(bo)])

% plot the orbital

% ensure that the 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 the expansion coefficients
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}')

butadiene
energies:  -1.618    -0.61803
resonance energy: -0.11803
bond orders: 0.89443     0.44721
Current plot held
Current plot released

Published with MATLAB® 7.14