%Fig2_7.m %Essential Electron Transport for Device Physics %tight-binding model in 2D %guided random walk towards objective function clear all; clf; %plotting parameters FS=12; %label fontsize LW=1; %curve linewidth MS=7.5; %marker size npo=20; %number of printouts for iteration loop % Change default axes fonts. set(0,'DefaultAxesFontName', 'Times'); set(0,'DefaultAxesFontSize', FS); % Change default text fonts. set(0,'DefaultTextFontname', 'Times'); set(0,'DefaultTextFontSize', FS); tic;%initiate tic-toc timer %switches plotting=1; %option to generate figures during optimization contplot=1; %option to generate contour plot at end x_initial=1; %'1'=use random atom positions, '2'=use estimated atom positions from figure %using x_initial=2 requires that A1=12,A2=4,Lx=Ly=L=80,dx=0.1 x_obj=2; %'1'=use randomly generated atomic positions as objective; %'2'=use custom atomic configuration %'3'=use test configuration (to determine interatomic distance cutoff) randtype=3; %type of randomization: %'1'=completely new random arrangement of all atoms %'2'=random walk, one atom selected at random %'3'=random walk, all atoms maxcount=1e5; %maximum number of iterations 1e4 Jth=0.0002; %acceptance threshold for cost function 0.005 %miscellaneous options periodic=1; %implement periodic boundary conditions ('1'=yes,'0'=no) allints=0; %option to include sum of interactions between given atom and %all other atoms in central cell and in every surrounding cell. options=[periodic,allints]; levels=60; %number of contour lines (levels) for contour plot 20 conthresh=0.50; %threshold for contour plot 0.2 %parameters A=8; %total number of atoms 8 t=1; %magnitude of hopping energy alpha=2; %long-range decay exponent 3 2 1 Gamma=0.2828*t;%characteristic energy broadening EX=4; %hopping strength multiplier (energy limits) 10 minE=-EX*t; %minimum energy 3 maxE= EX*t; %maximum energy -3 nE=1000; %size of energy array in density of states plot dminmaxE=(maxE-minE)/(nE); %prefactor density of states per unit energy interval E=linspace(minE,maxE,nE);%energy range for density of states plot L=40; %length of two-dimensional square domain 40 Lx=L; %length of 2D domain (along x-direction) Ly=L; %length of 2D domain (along y-direction) dx=0.1; %size of spatial increment (between lattice sites) 0.1 if x_initial==2 && isequal([A1,A2,Lx,Ly,dx],[12,4,80,80,0.1])==0 error(['If x_initial==2, the following must be true: A1=12,A2=4,Lx=Ly=80,dx=0.1']); else %store relevant parameters for parallel loop params=[randtype,A,Lx,Ly,t,alpha,dx]; end %objective function setup config_obj=sparse(Ly,Lx); if x_obj==1 %random starting configuration a=0; while a=minE && max(Esim)<=maxE Ecsw=0; end end J(1)=sum((Esim-Eobj).^2); %compute error % Esim = -Esim; %for whatever reason, the hopping energy seems to be %flipped in textbook figure NE=sum(Gamma2pi./((E-Esim).^2+Gamma2));%note factor 2 for spin NE=NE*dminmaxE; %density of states per unit energy interval %************************************************************************** %******************* main loop **************************** %************************************************************************** for k=1:maxcount %maxcount is bound on maximum number of iterations if mod(k,round(maxcount/npo))==0 percentdone=100*k/maxcount; disp([' Total iterations are ',num2str(percentdone),'% complete']); end if J(k)=minE && max(Esim)<=maxE Ecsw=0; end end %************************************************************************* %cost function is sum of differences in ordered eigenenergies squared J(k+1)=sum((Esim-Eobj).^2); %compute cost for new configuration %************************************************************************* %decide whether to keep new configuration if J(k+1)