%Fig2_8a.m %Essential Electron Transport for Device Physics %tight-binding model in 1D %random search for objective function clear all; clc; clf; %close all; %plotting parameters FS=12; %label fontsize LW=1; %curve linewidth MS=7.5; %marker size % Change default axes fonts. set(0,'DefaultAxesFontName', 'Times'); set(0,'DefaultAxesFontSize', FS); % Change default text fonts. set(0,'DefaultTextFontname', 'Times'); set(0,'DefaultTextFontSize', FS); %parameters A=3; %total number of atoms on lattice 6 9 t=-1; %hopping energy alpha=3; %long-range exponent 3 Gamma=0.3*abs(t);%characteristic energy broadening minE=-4*abs(t); %minimum energy 3 maxE=4*abs(t); %maximum energy -3 nE=1000; %size of energy array in density of states plot E=linspace(minE,maxE,nE);%energy range for density of states plot L=12; %length of one-dimensional domain 12 %setup emax= 2; %minimum energy in 1D band emin=-2; %maximum energy in 1D band Eobj=linspace(emin,emax,A)%evenly spaced eigenvalues as objective function %Eobj=-2*t*cos(linspace(0,pi,A))%nearest neighbor tight-binding % note: rely on fact eigenenergies are in assending order NEobj=zeros(1,nE); %initialize objective density of states array %************************************************************************** maxcount=100000; %maximum number of iterations Jth=10^-6; %acceptance threshold for cost function 0.1 0.03 0.01 0.001 H=zeros(A,A); %initialize H J=zeros(maxcount,1);%initialize J Gamma2=(Gamma/2).^2; Gamma2pi=2*Gamma/(2*pi);%note factor 2 for spin %************************************************************************* %calculate objective density of states over energy range of interest for i=1:A NEobj=NEobj+(Gamma2pi./((E-Eobj(i)).^2+Gamma2));%note factor 2 for spin end NEobj=NEobj*(maxE-minE)/(nE);%density of states per unit energy interval %initial atom positions x=L*rand(A,1); %************************************************************************** %******************* main loop **************************** %************************************************************************** for k=1:maxcount %maxcount is bound on maximum number of iterations for i=1:A for j=1:A if i==j H(i,j)=0;%do not count i=j contributions, set diagonal to zero else dx=abs(x(i)-x(j));% difference in position dx < L H(i,j)=-t/(dx^alpha)-t/((L-dx)^alpha);%periodic boundary conditions in 1D %H(i,j)=-t/(dx^alpha);%non-periodic pairs end end end Esim= eig(H);% note: rely on fact eigen energies are in order %************************************************************************* %cost function is sum of differences in ordered eigenenergies squared J(k)=sum((Esim'-Eobj).^2); %************************************************************************* %if JJ(i-1); J(i)=J(i-1); end end ttl1=['\rm Fig2.8a L=',num2str(L),', A=',num2str(A),... ', t=',num2str(t),', \alpha=',num2str(alpha),', J=',num2str(J(k)),... ', J_{th}=',num2str(Jth),', n_{count}=',num2str(k)]; ttl2=['\rm Fig2.8a obj(red), L=',num2str(L),', A=',num2str(A),... ', t=',num2str(t),', \Gamma=',num2str(Gamma),... ', \alpha=',num2str(alpha),', J=',num2str(J(k)),... ', J_{th}=',num2str(Jth),', n_{count}=',num2str(k)]; %************************************************************************* %plot figures figure(1) loglog(J);title(ttl1) grid on;hold on axis([1 maxcount 10^-7 10^4]);% %axis([1 maxcount 0 10000*Jth]); xlabel('Iteration','Fontsize',FS);ylabel('Convergence','Fontsize',FS); figure(2) plot(E,NEobj,'r');hold on;plot(E,NE);grid on;title(ttl2) xlabel('Energy, E/t','Fontsize',FS);ylabel('Density of states, N(E)','Fontsize',FS); hold off x=sort(x) %sorted values of atom positions x y=zeros(length(x),1); figure(3) plot(x,y,'or');grid on axis([0 L -1 1]); xlabel('Position, x','Fontsize',FS);title(ttl1); return %************************************************************************* %if k