runstructure.m

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

function runstructure

%RUNSTRUCTURE performs computation
%
%runstructure
%
%performs the self-consistent solution of Schroedinger and Poisson equation
%for the previously defined structure. If no PBOX is defined, no output
%during the  computation is generated. The result of this procedure are the
%global variables 'phi' and 'aquila_subbands' that contain the computed electrical
%potential and the wave functions and energies of the electrons in the quantum
%regions. See file 'structures' for a desription of aquila_subbands.

%Copyright 1999 Martin Rother
%
%This file is part of AQUILA.
%
%AQUILA is free software; you can redistribute it and/or modify
%it under the terms of the GNU General Public License as published by
%the Free Software Foundation; either version 2 of the License, or
%(at your option) any later version.
%
%AQUILA is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with AQUILA; if not, write to the Free Software
%Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

global phi aquila_control aquila_structure poimatrix

%tell the user, who we are
disp('AQUILA 1.0');
disp('(c) 1999 Martin Rother');

%we have to generate the structure first
buildstructure;

%get the size of the problem
nx=length(aquila_structure.xpos);
ny=length(aquila_structure.ypos);
n=nx*ny;

%generate the Poisson matrix. It does not change during the program run, so
%generate it only once. It is stored in the global variable 'poimatrix'.
genmatrixpoi
   
%redefine Fermi level in midgap of GaAs
%this is necessary, because the user may have changed the temperature
%after calling 'initaquila', where 'Efermi' has been computed before
aquila_control.Efermi=gaasmaterial(0,'E_c')-0.5*gaasmaterial(0,'E_G6G8');


%check, whether the given startpotential has correct size
%the potential 'phi' is already a all-zero matrix of correct size, generated in
%'buildstructure', which is used to start the iteration if no startpotential was given
if exist('aquila_control.phistart')
   if size(aquila_control.phistart)==size(phi) %size is OK
      phi=aquila_control.phistart; %set it as the startpotential
      if ~isempty(aquila_structure.qbox) %if we have QBOXes, enable the Schroedinger solver
         aquila_control.progress_check=bitset(aquila_control.progress_check,7);
      end
   else %size is wrong
      disp('runstructure: external startpotential has wrong size, using default all zero potential');
   end
end

%give some output, that the main thing is starting now
if aquila_control.verbose>0
   disp('runstructure: self-consistent classical Poisson solution')               
end

%if no QBOXes are defined then redefine values for stopping the iteration
if isempty(aquila_structure.qbox) 
   aquila_control.poisson.phitol=aquila_control.phitol;
   aquila_control.poisson.ftol=0.1*aquila_control.phitol;
   aquila_control.poisson.maxiter=aquila_control.maxiter;
end

%now we call the main solver
%The idea is: If we have no user-defined start potential, then use classical
%(non-linear Poisson equation) first. If the iteration has converged enough,
%then additionally incorporate the Schroedinger equation via a predictor-
%corrector scheme. If the user has given a startpotential, we expect it to
%be good enough the start the Schroedinger-Poisson solution immediately
%via the predictor-corrector scheme.
if ~exist('aquila_control.phistart') %we have no user-defined startpotential
   %some output
   if (aquila_control.verbose>1)|(~isempty(aquila_structure.pbox))
      os=sprintf('Classical computation');
      disp(os)
   end
   newton(1,[]); %perform classical computation first
else
   schrsolve(phi,0); %solve Schroedingers equation for the user-defined startpotential
   %some output
   if (aquila_control.verbose>1)|(~isempty(aquila_structure.pbox))
      os=sprintf('Schroedinger step 0');
      disp(os)
   end
   newton(0,phi); %perform Newton scheme. Same as classical computation, but already with
                  %quantum mechanical charge density
end
%if we do classical computation only, we are finished now

%otherwise we have a good potential now and are ready to switch on the Schroedinger solver
%if it has not already been switched on.
%This of course makes sense only, when QBOXes are defined
if ~isempty(aquila_structure.qbox)
   %tell the user, what we are doing
   if aquila_control.verbose>0
      disp('runstructure: starting self-consistent quantum Poisson solution')               
   end
   
   %prepare the predictor-corrector scheme
   
   phi_old=phi; %save the potential
   schrsolve(phi,0); %compute eigenvalues
   %tell the user
   if (aquila_control.verbose>1)|(~isempty(aquila_structure.pbox))
      os=sprintf('Schroedinger step %d',1);
      disp(os)
   end
   %set flag, indicating that we now have eigenvalues
   aquila_control.progress_check=bitset(aquila_control.progress_check,7); 
   newton(0,phi_old); %make some modified Newton steps with predictor charge density
   %we now have a new potential in 'phi'
   rho=max(max(abs(phi-phi_old))); %maximum change in potential
   
   %stop, if the iteration has converged enough
   %this happens at this early stage in the computation, if there is no charge in the
   %structure. The self-consistency is already reched at this point.
   if rho<aquila_control.phitol
      if aquila_control.verbose>1
         disp('runstructure: convergence reached by phitol');
      end
      return
   end

   %now comes the outer loop of the predictor-corrector scheme
   itr=1; %reset the iteration counter
   while itr<aquila_control.maxiter %too many iterations ?
      itr=itr+1;
      phi_old=phi;
      
      %now compute new eigenvalues
      %if the change in the potential is small, then use inverse vectoriteration
      %otherwise use MATLABs 'eigs' (Arnoldi algorithm ?)
      if rho>aquila_control.tracklimit
         schrsolve(phi,0);
      else
         schrsolve(phi,1);
      end
      %some output
      if (aquila_control.verbose>1)|(~isempty(aquila_structure.pbox))
         os=sprintf('Schroedinger step %d',itr);
         disp(os)
      end
      
      %we now have new eigenvalues (corrector step)
      %now lets do a few Newton steps with them (predictor step)
      newton(0,phi);
      
      %the maximum change in the potential
      rho=max(max(abs(phi-phi_old)));
            
      %stop, if the iteration has converged enough
      if rho<aquila_control.phitol
         if aquila_control.verbose>1
            disp('runstructure: convergence reached by phitol');
         end
         return
      end
   end
   
   %stop, we have no convergence
   disp('runstructure: no convergence reached at maxiter!');
end
return

%************************************
%some subroutines doing the main work
%************************************

function newton(qstop,phi_last)

%perform modified Newtons method
%in this iteration the error of non-linear Poissons equation is
%minimized. We do not directly minimize the error but the error squared
%thus giving global convergence.
%
%qstop is a flag. qstop=1 means, stop the iteration when the change in the
%potential becomes small enough to activate the Schroedinger solver
%
%phi_last is the potential from the last iteration. This is the potential
%for which the eigenvalues have been computed and is needed here to
%approximate the eigenvalues of the new potential during the iteration
%without a complete recomputation of the spectrum.
%
%the result of this routine is a new electric potential which
%is returned in the global variable phi

global phi aquila_control aquila_structure poimatrix

nx=length(aquila_structure.xpos);
ny=length(aquila_structure.ypos);
n=nx*ny;

%make the matrix of the nodes of the electric potential a vector
phi=phi';
phi=phi(:);
dphi=aquila_control.schrenable;

%compute the error of Poissons equation for potential phi
%fvec is the error vector itself, rho is 0.5*fvec'*fvec
%phi_last is the potential for wich the eigenvalues have been computed
%and is needed here to find an approximation of the eigenvalues of the
%actual potential phi
[rho,fvec]=fmin(phi,phi_last);
      
%reset the iteration counter and a flag
check=0;
itr=0;

stpmax=100*max(norm(phi),length(phi));   
%the iteration loop
while itr<aquila_control.poisson.maxiter %limit the number of iterations
   itr=itr+1;

   phix=reshape(phi,nx,ny)'; %store the potential in matrix form
   
   %display the progress of the iteration for the user
   if itr>1
      showprogress(phix,dphi,itr-1);
   end   
   
   %return if limit for Schrodinger solver is reached
   if (dphi<aquila_control.schrenable)&(bitget(aquila_control.progress_check,7)==0)&...
         (~isempty(aquila_structure.qbox))&(qstop==1)
      if aquila_control.verbose>1
         disp('runstructure: returning for Schroedinger solution');
      end
      phi=phix; %make phi a matrix again
      return
   end
   
   %compute derivative of the square of the error of the Poisson equation
   dF=poimatrix;
   
   %compute the derivative of the charge
   if bitget(aquila_control.progress_check,7)>0 %if already have wavefunctions
      ch=sumcharge(phix,phix-phi_last,1);%derivative of quantum charge
   else
      ch=sumcharge(phix,1);%derivative of charge
   end
   %generate the right hand side of Poissons equation and subtract it from the main diagonal
   dF=spdiags(spdiags(poimatrix,0)-genrhspoi(ch,1),0,dF);
   
   gradf=fvec'*dF; %the gradient
   
   deltaphi=-dF\fvec; %solve the system and compute the correction vector
         
   phiold=phi; %save phi
   %we adapt the length of the correction vector for the potential
   %by a line minimalization algorithm to emsure, that the error of
   %Poissons equation decreases.
   %The algorithm returns the new potential, the square of the error, the errorvector
   %and a flag
   if rho~=0
      [phi,rho,fvec,check]=linesearch(phi,rho,gradf,deltaphi,stpmax,phi_last);
   end
   
   %stop, if the error of Poissons equation has dropped below the desired limit
   if (max(abs(fvec))<aquila_control.poisson.ftol)&(itr>1)
      phi=reshape(phi,nx,ny)'; %make the potential a matrix again
      if aquila_control.verbose>1
         disp('runstructure: convergence reached by poisson.ftol');
      end
      return
   end