runstructure.m

AQUILA is a MATLAB toolbox for the one- or twodimensional simulation of ％the electronic properties

   
   %the flag is set. This means, that the algorithm has been trapped in a
   %local minimum and will not converge further. If this happens during the
   %predictor-corrector loop, then after the return from 'newton' new eigenvalues
   %are computed which normally cures the problem.
   if check==1
      temp=max(abs(gradf)'.*max(abs(phi),ones(size(phi)) ))/max(rho,0.5*n);
      if temp<1e-6
         check=1;
      else
         check=0;
      end
      phi=reshape(phi,nx,ny)';
      disp('runstructure: local minimum reached, no further progress!')
      return
   end
   
   %stop, if the change in the potential has dropped below the desired limit
   if (max(abs(phi-phiold)./max(abs(phi),ones(size(phi))))<aquila_control.poisson.phitol)&(itr>1)
      phi=reshape(phi,nx,ny)';
      if aquila_control.verbose>1
         disp('runstructure: convergence reached by poisson.phitol');
      end
      return
   end
   
   dphi=max(abs(phi-phiold));%save the change in the potential
   
   %again some output for the user
   if aquila_control.verbose>1
      os=sprintf('Poisson-Iteration %d, errf=%g, deltaphimax=%g',itr,rho,dphi);
      disp(os)
   end
end
%if we get here, we have not reached convergence
disp('runstructure: no convergence reached at poisson.maxiter!');
phi=reshape(phi,nx,ny)';
return

%***********************************
% here follow some helper routines
%***********************************

function showprogress(phix,dphi,itr)

%show progress information on screen

global aquila_structure aquila_control aquila_material
constants
sp=size(aquila_structure.pbox);
sp=sp(1);
%if we have output at all
if sp>0
   %textual output
   os=sprintf('Newton step %d: max(deltaphi)=%g eV',itr,dphi);
   disp(os)
   
   %graphical output
   %compute charge density, conduction band and valence band
   ch=sumcharge(phix);
   cb=aquila_material.ec-phix;
   vb=aquila_material.ev-phix;
   %how many plots do we need
   plotnr=length(unique([find(aquila_structure.pbox(:,1)-aquila_structure.pbox(:,3)==0)' ...
         find(aquila_structure.pbox(:,2)-aquila_structure.pbox(:,4)==0)']));
   for count=1:sp %for all PBOXes
      
      %PBOX is a slice in y-direction
      if (aquila_structure.pbox(count,1)==aquila_structure.pbox(count,3))&...
            (aquila_control.mode==2)
         %find positions
         ix=find(aquila_structure.xpos<=aquila_structure.pbox(count,1));
         ix=ix(end);
         iy=intersect(find(aquila_structure.pbox(count,2)<=aquila_structure.ypos),...
            find(aquila_structure.pbox(count,4)>=aquila_structure.ypos));         
         %plot the band diagram
         subplot(2,plotnr,count);
         pl=ones(length(iy),1)*aquila_control.Efermi-aquila_material.bias(iy,ix); %the Fermi energy
         if bitand(aquila_structure.pbox(count,5),CB)>0
            pl=[cb(iy,ix) pl]; %add conduction band to plot if desired
         end
         if bitand(aquila_structure.pbox(count,5),VB)>0
            pl=[pl vb(iy,ix)]; %add valence band to plot if desired
         end         
         plot(aquila_structure.ypos(iy),pl); %draw and label the plot     
         os=sprintf('bands at x=%d A',aquila_structure.pbox(count,1));
         title(os);
         xlabel('y-position [Angstrom]');
         ylabel('energy [eV]');
         %plot the charge density
         subplot(2,plotnr,count+plotnr);
         plot(aquila_structure.ypos(iy),ch(iy,ix)*1e24); %1e24 scales A^-3 to cm^-3
         os=sprintf('charge density at x=%d A',aquila_structure.pbox(count,1));
         title(os);
         xlabel('y-position [Angstrom]');
         ylabel('density [cm^-3]');
         %textual output of charge sheet density along slice
         %get correct width of the slice line
         if ix==1
            width=aquila_structure.hx(1)/2;
         elseif ix==length(aquila_structure.xpos)
            width=aquila_structure.hx(end)/2;
         else
            width=(aquila_structure.hx(ix)+aquila_structure.hx(ix-1))/2;
         end
         %integrate charge along slice line and form sheet density
         %factor 1e16 scales A^-2 to cm^-2
         tcharge=sum(ch(iy,ix).*aquila_structure.boxvol(iy,ix))/width*1e16;
         os=sprintf('%d<=y<=%d A at x=%d A sheet density %g cm^-2',...
            aquila_structure.pbox(count,2),aquila_structure.pbox(count,4),...
            aquila_structure.pbox(count,1),tcharge);
         disp(os);
         
         %PBOX is slice in x-direction or we have 1D simulation only
      elseif (aquila_structure.pbox(count,2)==aquila_structure.pbox(count,4))|...
            (aquila_control.mode==1)
         %find correct y-index and x-index
         if aquila_control.mode==2
            iy=find(aquila_structure.ypos<=aquila_structure.pbox(count,2));
            iy=iy(end);
         else
            iy=1;
         end
         ix=intersect(find(aquila_structure.pbox(count,1)<=aquila_structure.xpos),...
            find(aquila_structure.pbox(count,3)>=aquila_structure.xpos));
         %plot the band diagram
         subplot(2,plotnr,count);
         pl=ones(length(ix),1)*aquila_control.Efermi-aquila_material.bias(iy,ix)'; %Fermi energy
         if bitand(aquila_structure.pbox(count,5),CB)>0
            pl=[cb(iy,ix)' pl]; %add conduction band to plot if desired
         end
         if bitand(aquila_structure.pbox(count,5),VB)>0
            pl=[pl vb(iy,ix)']; %add valence band to plot if desired
         end         
         plot(aquila_structure.xpos(ix),pl); %draw and label the plot
         os=sprintf('bands at y=%d A',aquila_structure.pbox(count,2));
         title(os);
         xlabel('x-position [Angstrom]');
         ylabel('energy [eV]');
         %plot the charge density
         subplot(2,plotnr,count+plotnr);
         plot(aquila_structure.xpos(ix),ch(iy,ix)*1e24); %1e24 scales A^-3 to cm^-3
         if aquila_control.mode==2
            os=sprintf('charge density at y=%d A',aquila_structure.pbox(count,2));
         else
            os='charge density'; %for 1D simulation
         end
         title(os);
         xlabel('x-position [Angstrom]');
         ylabel('density [cm^-3]');
         %textual output of charge sheet density along slice
         %compute correct width
         if iy==1
            if aquila_control.mode==2
               width=aquila_structure.hy(1)/2;
            else
               width=1;
            end
         elseif iy==length(aquila_structure.ypos)
            width=aquila_structure.hy(end)/2;
         else
            width=(aquila_structure.hy(iy)+aquila_structure.hy(iy-1))/2;
         end
         %integrate along slice line
         %1e16 scales A^-2 to cm^-2
         tcharge=sum(ch(iy,ix).*aquila_structure.boxvol(iy,ix))/width*1e16;
         if aquila_control.mode==2
            os=sprintf('%d<=x<=%d A at y=%d A sheet density %g cm^-2',...
               aquila_structure.pbox(count,1),aquila_structure.pbox(count,3),...               
               aquila_structure.pbox(count,2),tcharge);
         else
            os=sprintf('%d<=x<=%d A sheet density %g cm^-2',...
               aquila_structure.pbox(count,1),aquila_structure.pbox(count,3),tcharge);
         end   
         disp(os);
      else
         %textual output only, if PBOX is a real box and not a slice
         ix=intersect(find(aquila_structure.pbox(count,1)<=aquila_structure.xpos),...
            find(aquila_structure.pbox(count,3)>=aquila_structure.xpos));
         iy=intersect(find(aquila_structure.pbox(count,2)<=aquila_structure.ypos),...
            find(aquila_structure.pbox(count,4)>=aquila_structure.ypos));
         %integrate over PBOX, 1e8 scales A^-1 to cm^-1
         tcharge=sum(sum(ch(iy,ix).*aquila_structure.boxvol(iy,ix)))*1e8;
         os=sprintf('%d<=x<=%d A, %d<=y<=%d A line density %g cm^-1',...
            aquila_structure.pbox(count,1),aquila_structure.pbox(count,3),...
            aquila_structure.pbox(count,2),aquila_structure.pbox(count,4),tcharge);
         disp(os);
      end
   end
   drawnow %draw the plot
end

%*****************************

function [phi,rho,fvec,check]=linesearch(phiold,rhoold,gradf,deltaphi,stpmax,phi_last)

%linesearch algorithm for Newtons method
%parts of this routine have been taken from Numerical Recipes

global aquila_control
check=0;
no=norm(deltaphi);
if no>stpmax
   deltaphi=deltaphi*(stpmax/no);
end
slope=gradf*deltaphi;
lambdamin=1e-8./max(abs(deltaphi)./max(abs(phiold),ones(size(phiold))));
lambda=1;
count=0;
while 1
   count=count+1;
   phi=phiold+lambda*deltaphi;%try full Newton step
   [rho,fvec]=fmin(phi,phi_last);
   if lambda<lambdamin
      phi=phiold;
      check=1;
      disp('linesearch: lambdamin reached');
      return
   elseif rho<rhoold+1e-4*lambda*slope %1e-4 incoporates average rate of decrease
      if aquila_control.verbose>1
         disp('linesearch: sufficient function decrease');
      end
      return
   else
      if lambda==1 %first iteration
         tmplambda=-slope/(2*(rho-rhoold-slope));%quadratic approximation
      else
         rhs1=rho-rhoold-lambda*slope;
         rhs2=rho2-rhoold2-lambda2*slope;
         a=(rhs1/lambda^2-rhs2/lambda2^2)/(lambda-lambda2);
         b=(-lambda2*rhs1/lambda^2+lambda*rhs2/lambda2^2)/(lambda-lambda2);
         if a==0
            tmplambda=-slope/(2*b);
         else
            tmplambda=(-b+sqrt(b*b-3*a*slope))/(3*lambda);
         end
         if tmplambda>0.5*lambda %no too big iteration steps
            tmplambda=0.5*lambda;
         end
      end
   end
   lambda2=lambda;
   rho2=rho;
   rhoold2=rhoold;
   lambda=max(tmplambda,0.1*lambda); %no too small iteration steps
   if aquila_control.verbose>1
      os=sprintf('It2=%d, err=%g, lambda=%g',count,rho,lambda);
      disp(os)
   end
end

%*****************************

function [ferr,F]=fmin(phi,phi_last)

%compute error in Poisson equation and is square

F=poierr(phi,phi_last);
ferr=0.5*norm(F)^2;

%*****************************

function F=poierr(phi,phi_last)

%compute error in Poisson equation

global aquila_structure aquila_control poimatrix
nx=length(aquila_structure.xpos);
ny=length(aquila_structure.ypos);
n=nx*ny;
phi=reshape(phi,nx,ny)'; %make phi a matrix again
%if we have eigenvalues, use predictor density, otherwise use classical density
if bitget(aquila_control.progress_check,7)>0
   ch=sumcharge(phi,phi-phi_last);
else
   ch=sumcharge(phi);
end
rhs=genrhspoi(ch); %form Poissons equation right hand side
phi=phi';
phi=phi(:);
F=poimatrix*phi-rhs; %the error of Poisson equation