DOS.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - DOS
%    Copyright (C) 2016 Peter Rakyta, Ph.D.
%
%    This program 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 3 of the License, or
%    (at your option) any later version.
%
%    This program 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 this program.  If not, see http://www.gnu.org/licenses/.
%
%> @addtogroup utilities Utilities
%> @{
%> @file DOS.m
%> @brief A class to calculate the density of states along a one-dimensional energy array.
%> @} 
%> @brief A class to calculate the density of states along a one-dimensional energy array.
%> @Available
%> EQuUs v4.9 or later
%%
classdef DOS < UtilsBase
  
    
    properties (Access = public)
        % A structure containing the calculated DOS.
        DOSdata
        
    end
    
    
methods (Access=public)
    
%% Contructor of the class
%> @brief Constructor of the class.
%> @param Opt An instance of structure #Opt.
%> @param varargin Cell array of optional parameters. For details see #InputParsing.
%> @return An instance of the class
    function obj = DOS( Opt, varargin )     
  
        obj = obj@UtilsBase( Opt, varargin{:} );  
        
        obj.DOSdata = structures('DOS');
        
        if strcmpi(class(obj), 'DOS') 
            obj.InputParsing( obj.varargin{:} );
        end
        
        
    end    

    
    
    
%% DOSCalc
%> @brief Calculates the onsite desnity of states.
%> @param Evec One-dimensional array of the energy points.
%> @param varargin Cell array of optional parameters:
%> @param 'JustScatter' Logical value. Set true to omit the contacts from the system
    function DOSCalc( obj, Evec, varargin ) 
        
        p = inputParser;
        p.addParameter('JustScatter', false ); %logical value. True if only an isolated scattering center should be considered in the calculations, false otherwise.
        
        p.parse(varargin{:});      
        JustScatter = p.Results.JustScatter;

        obj.create_Hamiltonians( 'junction', obj.junction )
        Hscatter = obj.junction.CreateH.Read('Hscatter'); 
        
        
        obj.display(['EQuUs:Utils:', class(obj), ':DOSCalc: Calculating the DOS between Emin = ', num2str(min(Evec)),' and Emax = ',num2str(max(Evec)),],1);        
       
              
        %> preallocate memory for Densitysurf
        DOSvec = complex(zeros( size(Evec) ));
        
        % starting parallel pool
        poolmanager = Parallel( obj.junction.FL_handles.Read('Opt') );
        poolmanager.openPool(); 

        %> creating function handles for parallel for
        hdisplay = @obj.display;
        hSetEnergy = @obj.SetEnergy;
        hcomplexSpectral = @obj.complexSpectral;
        hcreate_scatter = @obj.create_scatter;
        
        junction_loc = obj.junction;
        
        %for iidx = 1:length(Evec)
        parfor (iidx = 1:length(Evec), poolmanager.NumWorkers)
            
            %> setting the current energy
            feval(hSetEnergy, Evec(iidx), 'junction', junction_loc );  
            
            %> creating the Hamiltonian of the scattering region
            feval(hcreate_scatter, 'junction', junction_loc );
           
            try 
                %> evaluating the energy resolved density
                densityvec2int = feval(hcomplexSpectral, junction_loc, JustScatter);
            catch errCause
               feval(hdisplay, ['EQuUs:Utils:', class(obj), ':DOSCalc: Error at energy point E = ', num2str(Evec(iidx)), ' caused by:'], 1);
               feval(hdisplay, errCause.identifier, 1 );
                for jjjdx = 1:length(errCause.stack)
                    feval(hdisplay, errCause.stack(jjjdx), 1 );
                end
               feval(hdisplay, ['EQuUs:Utils:', class(obj), ':DOSCalc: Giving NaN for the calculated current at the given energy point.'], 1); 
                densityvec2int = NaN;
            end
                
            DOSvec( iidx ) = -imag( sum( densityvec2int ) )/pi;
            
        end

        %> closing the parallel pool
        poolmanager.closePool();
        
        % setting the attribute values
        obj.DOSdata.DOS = DOSvec;
        obj.DOSdata.Evec = Evec;
          
        
        obj.display(['EQuUs:Utils:', class(obj), ':DOSCalc: Process finished.'], 1)
        
        %reset the system  for a new calculation
        obj.InputParsing( obj.varargin{:} );
        
        





         
    end
    
    
%% LocalDOSCalc
%> @brief Calculates the local desnity of states for a given energy.
%> @param Energy The energy value in the units of the Hamiltonians
%> @param varargin Cell array of optional parameters:
%> @param 'JustScatter' Logical value. Set true to omit the contacts from the system
%> @param 'ChoseSites' Function handle ret = ChoseSites( #Coordinates ) returning an array of logical values. The local DOS is calculated for those sites in the scattering region, which are labeled with true in the array ret. (ret should be the same size as the x,y,z arrays in the input class #Coordinates.)
%> @return Returns with the calculated local DOS. (An instance of structure #LocalDOS)
    function cLocalDOS = LocalDOSCalc( obj, Energy, varargin ) 
        
        p = inputParser;
        p.addParameter('JustScatter', false ); %logical value. True if only an isolated scattering center should be considered in the calculations, false otherwise.
        p.addParameter('ChoseSites', [] );
        
        p.parse(varargin{:});      
        JustScatter = p.Results.JustScatter;
        ChoseSites  = p.Results.ChoseSites;

        obj.create_Hamiltonians( 'junction', obj.junction )
        Hscatter = obj.junction.CreateH.Read('Hscatter'); 
        cCoordinates_scatter = obj.junction.CreateH.Read('coordinates'); 
        
        
        obj.display(['EQuUs:Utils:', class(obj), ':LocalDOSCalc: Calculating the local DOS at energy E = ', num2str(Energy)],1);        
       
                      
        junction_loc = obj.junction;
        
        % setting the current energy
        obj.SetEnergy( Energy, 'junction', junction_loc );  
            
        % creating the Hamiltonian of the scattering region
        obj.create_scatter( 'junction', junction_loc );
        
        % setting the function handle for displaying the message
        hdisplay = @obj.display;
           
        try 
            % evaluating the energy resolved density
            densityvec2int = obj.complexSpectral( junction_loc, JustScatter, 'ChoseSites', ChoseSites);
        catch errCause
          obj.display( ['EQuUs:Utils:', class(obj), ':LocalDOSCalc: Error at energy point E = ', num2str(Energy), ' caused by:'], 1);
               feval(hdisplay, errCause.identifier, 1 );
            for jjjdx = 1:length(errCause.stack)
                obj.display( errCause.stack(jjjdx), 1 );
            end
               obj.display( ['EQuUs:Utils:', class(obj), ':LocalDOSCalc: Giving NaN for the calculated current at the given energy point.'], 1); 
            densityvec2int = NaN;
        end
        
        DOSvec = -imag( densityvec2int )/pi;
        
        
        % creating the data structure of the local density of states
        cLocalDOS = structures( 'LocalDOS');
        
        cLocalDOS.Energy = Energy;
        cLocalDOS.DOSvec = DOSvec;
        
        if ~isempty( ChoseSites )
            cCoordinates_scatter = cCoordinates_scatter.KeepSites( ChoseSites( cCoordinates_scatter ) );
        end
          cLocalDOS.coordinates = cCoordinates_scatter;
        
        obj.display(['EQuUs:Utils:', class(obj), ':DOSCalc: Process finished.'], 1)
        
        %reset the system  for a new calculation
        obj.InputParsing( obj.varargin{:} );
        
        





         
    end    
  
    
    
    end % public methods
    
    
    methods (Access=protected)  
        
        
%% complexSpectral
%> @brief Calculates the spectral function at a complex energy
%> @param junction_loc An instance of class NTerminal or its subclass representing the junction.
%> @param JustScatter logical value. True if only an isolated scattering center should be considered in the calculations, false otherwise.
%> @param varargin Cell array of optional parameters:
%> @param 'ChoseSites' Function handle ret = ChoseSites( #Coordinates ) returning an array of logical values. The local DOS is calculated for those sites in the scattering region, which are labeled with true in the array ret. (ret should be the same size as the x,y,z arrays in the input class #Coordinates.)
%> @return spectral_function The calculted spectral function and a data structure containing geometry data of the junction.
     function [spectral_function, junction_sites] = complexSpectral( obj, junction_loc, JustScatter, varargin )
        
        p = inputParser;
        p.addParameter('ChoseSites', [] );
        
        p.parse(varargin{:});      
        ChoseSites  = p.Results.ChoseSites;
        
        % Obtaining the Hamiltonian of the scattering center
        Hscatter = junction_loc.CreateH.Read('Hscatter');
        Hscatter_transverse = junction_loc.CreateH.Read('Hscatter_transverse');
        
        q = junction_loc.CreateH.Read('q');
        if ~isempty( q ) && ~junction_loc.CreateH.Read('HamiltoniansDecimated')
            Hscatter = Hscatter + Hscatter_transverse*diag(exp(1i*q)) + Hscatter_transverse'*diag(exp(-1i*q)); 
        end
        
        % creating the inverse of the Green opertor related to the scattering center
        gfininv = speye(size(Hscatter))*junction_loc.getEnergy()-Hscatter;
        
        
            
        if JustScatter
            % in case when just the scattering region is involved in the calculations
            Ginverz = gfininv;
        else                            
            % Use Dyson equation to attach contacts to the scattering center
            [~, Ginverz, junction_sites] = junction_loc.CustomDysonFunc( 'UseHamiltonian', true, 'onlyGinverz', true, ...
                'decimate', false, 'gfininv', gfininv, 'constant_channels', false, ...
                'SelfEnergy', obj.useSelfEnergy, 'keep_sites', 'scatter');            
        end
        
        if isempty( ChoseSites )
            % if no geometrical constriction on the sites of the scatter
            % region, then all the sites are kept in the calculations
            spectral_function = obj.diagInv( Ginverz );   
            spectral_function = spectral_function(end-size(Hscatter,1)+1:end);
        else
            
            % deterirmine sites in the scattering region to be calculated
            cCoordinates_scatter = junction_loc.CreateH.Read('coordinates'); 
            ChosenSites_scatter = ChoseSites( cCoordinates_scatter );
            
            % determine sites to be calculated in the full system
            indexes2calc = false( size(Ginverz,1), 1 );
            indexes2calc( end-size(Hscatter,1)+1:end) = ChosenSites_scatter;
            
            Ginverz = [ Ginverz( ~indexes2calc, ~indexes2calc) , Ginverz(~indexes2calc, indexes2calc);
                     Ginverz(indexes2calc, ~indexes2calc), Ginverz(indexes2calc, indexes2calc)];
            
            try
                spectral_function = obj.partialInv( Ginverz, length(cCoordinates_scatter.x(ChosenSites_scatter)) );
                spectral_function = diag(spectral_function);
            catch
                spectral_function = NaN(length(cCoordinates_scatter.x(ChosenSites_scatter)), 1);
            end
        end
             

    end        


%% create_Hamiltonians
%> @brief Creates the Hamiltonians of the system
%> @param varargin Cell array of optional parameters:.
    function create_Hamiltonians(  obj, varargin )
        
        p = inputParser;
        p.addParameter('junction', obj.junction); %for parallel computations
        p.parse(varargin{:});       
        junction_loc     = p.Results.junction;     
        
        
        % cerating Hamiltonian of the scattering region
        obj.create_scatter( 'junction', obj.junction )
        
        % creating the Lead Hamiltonians
        supClasses = superclasses(obj.junction);
        if sum( strcmp( supClasses, 'Ribbon' ) ) > 0
            coordinates_shift = [-2, junction_loc.height+1] + junction_loc.shift; 
            junction_loc.FL_handles.LeadCalc('coordinates_shift', coordinates_shift, 'transversepotential', junction_loc.transversepotential, ...
                'gauge_field', junction_loc.gauge_field, 'q', junction_loc.getTransverseMomentum(), 'leadmodel', junction_loc.leadmodel);
        else
            junction_loc.FL_handles.LeadCalc( 'CustomHamiltonian', junction_loc.cCustom_Hamiltonians , 'gauge_field', junction_loc.gauge_field, 'q', junction_loc.getTransverseMomentum(), 'leadmodel', junction_loc.leadmodel);
        end
              
        
  
        
    end        


%% create_scatter
%> @brief Creates the Hamiltonian of the scattering center
%> @param varargin Cell array of optional parameters:.
    function create_scatter(  obj, varargin )
        
        p = inputParser;
        p.addParameter('gauge_trans', true); % logical: true if want to perform gauge transformation on the Green's function and Hamiltonians
        p.addParameter('junction', obj.junction); %for parallel computations
        p.parse(varargin{:});       
        gauge_trans     = p.Results.gauge_trans;
        junction_loc     = p.Results.junction;     
        
        % create the unit cell of the scattering center
        junction_loc.CreateScatter();
        
        % apply custom potential in the scattering center
        if ~isempty(obj.scatterPotential)          
            junction_loc.ApplyPotentialInScatter( junction_loc.CreateH, obj.scatterPotential)  
        end             
                   
                
        % gauge transformation of the calculated effective Hamiltonians
        if gauge_trans && ~isempty(  junction_loc.PeierlsTransform_Scatter )
            Hscatter = junction_loc.CreateH.Read('Hscatter');
            coordinates_scatter = junction_loc.CreateH.Read('coordinates');
            Hscatter = junction_loc.PeierlsTransform_Scatter.gaugeTransformation( Hscatter, coordinates_scatter, junction_loc.gauge_field );
            junction_loc.CreateH.Write('Hscatter', Hscatter);
            %ribbon_loc.PeierlsTransform_Scatter.gaugeTransformationOnLead( ribbon_loc.Scatter_UC, ribbon_loc.gauge_field );  
        end
  
        
    end




%% InputParsing
%> @brief Parses the optional parameters for the class constructor.
%> @param varargin Cell array of optional parameters:
%> @param 'T' The absolute temperature in Kelvins.
%> @param 'scatterPotential' A function handle y=f( #coordinates coords ) or y=f( #CreateHamiltonians CreateH, E ) of the potential across the junction.
%> @param 'useSelfEnergy' Logical value. Set true (default) to solve the Dyson equation with the self energies of the leads, or false to use the surface Green operators.
%> @param 'gfininvfromHamiltonian' Logical value. Set true calculate the surface Greens function of the scattering region from the Hamiltonaian of the scattering region, or false (default) to calculate it by the fast way (see Phys. Rev. B 90, 125428 (2014) for details).
%> @param 'junction' An instance of a class #NTerminal or its subclass describing the junction.
    function InputParsing( obj, varargin )
        
        p = inputParser;
        p.addParameter('T', obj.T);
        p.addParameter('scatterPotential', []);
        p.addParameter('gfininvfromHamiltonian', false);
        p.addParameter('useSelfEnergy', true);
        p.addParameter('junction', []);     
        
        p.parse(varargin{:});
        
        InputParsing@UtilsBase( obj, 'T', p.Results.T, ...
                            'junction', p.Results.junction, ...
                            'scatterPotential', p.Results.scatterPotential, ...
                            'gfininvfromHamiltonian', p.Results.gfininvfromHamiltonian, ...
                            'useSelfEnergy', p.Results.useSelfEnergy);

               

    end % 
    
    
end % protected methods
end