%Fig2_8b.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 3 6 9 t=-1; %hopping energy alpha=3; %long-range decay exponent 3 Gamma=0.3*abs(t); %characteristic energy broadening L=12; %length of one-dimensional domain 4*A 10 20 %setup emax=2*abs(t); %minimum energy in 1D band emin=-2*abs(t); %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 %************************************************************************** maxcount=100000; %maximum number of iterations x=L*rand(A,1); %initialize atom positions x H=zeros(A,A); %initialize Hamiltonian H J=zeros(maxcount,1); %iniitialize cost function J Jav=zeros(maxcount,1); %initialize average cost function Jav %************************************************************************** %******************* main loop **************************** %************************************************************************** kavmax=100; %average over 100 runs for kav=1:kavmax for k=1:maxcount %maxcount is bound on maximum number of iterations %************************************************************************** % call EvalEig to find simulated eigenvalues, Esim %************************************************************************** Esim = EvalEig(x,A,t,alpha,L); %************************************************************************* %cost function is sum of differences in ordered eigenenergies squared J(k) = sum((Esim'-Eobj).^2);%distance is abs(Esim'-Eobj) or (Esim'-Eobj).^2 %************************************************************************* %************************************************************************* x=L*rand(A,1); end %prepare cost function for convergence plot for i=2:k if J(i)>J(i-1); J(i)=J(i-1); end end Jav=(Jav+J); end ttl=['\rm Fig2.8b L=',num2str(L),', A=',num2str(A),', t=',num2str(t),... ', \alpha=',num2str(alpha),', J=',num2str(J(k)),... ', n_{count}=',num2str(k),', n_{av}=',num2str(kavmax)]; %************************************************************************* %plot figure Jav from 10th interation to maxcount iteration %************************************************************************* figure(1) loglog(Jav/kavmax);title(ttl) grid on;hold on; axis([1 maxcount 10^-7 10^4]);% %axis([1 maxcount 0 Jav(10)]);% xlabel('Iteration','Fontsize',FS);ylabel('Convergence','Fontsize',FS); %************************************************************************* % function EvalEig to find simulated eigenvalues, Esim %************************************************************************* function Esim = EvalEig(x,A,t,alpha,L) Hltrig=zeros(A,A);%initialize lower-triangle and ensure H diagonal is zero for i=2:A for j=1:i-1 dx=abs(x(i)-x(j));% difference in position dx < L Hltrig(i,j)=-t/(dx^alpha)-t/((L-dx)^alpha);%periodic boundary conditions in 1D end end %Hutrig=Hltrig', note upper-triangle is transpose of lower-triangle H=Hltrig+Hltrig';%H with diagonals zero Esim = eig(H); %eigenergies note: rely on fact eigenenergies are in order end