Density.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - Density
%    Copyright (C) 2009-2017 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 Density.m
%> @brief A class to calculate the onsite desnity useful for self-consistent calculations.
%> @} 
%> @brief A class to calculate the onsite desnity useful for self-consistent calculations.
%> @Available
%> EQuUs v4.9 or later
%%
classdef Density < DOS %NEW superclass

    properties ( Access = public )
        
    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 = Density( Opt, varargin )        

        obj = obj@DOS( Opt, varargin{:} );  
        
        if strcmpi(class(obj), 'Density') 
            obj.InputParsing( obj.varargin{:} );
        end
        
        
    end



%% DensityCalc
%> @brief Calculates the onsite desnity using the method in Nanotechnology 25 (2014), 465201
%> @param varargin Cell array of optional parameters:
%> @param 'Edb' The number of the energy points over the contour integral.
%> @param 'DeltaPhi' A parameter to control the incident angle of the integration contour near the real axis. (Default value is $$\Delta\Phi=0.5\pi$$.
%> @param 'Evec' The complex energy points in case of custom integration path. (Overrides all other optional parameters.)
%> @param 'Emin' The minimum of the energy array.
%> @param 'JustScatter' Logical value. True if only an isolated scattering center should be considered in the self-consistent calculations, false otherwise.
%> @return An array of the calculated density and a structure describing the geometry of the sites.
    function [density_vec, junction_sites] = DensityCalc( obj, varargin )
        
        p = inputParser;
          p.addParameter('Edb', 511 ); % number of energy points on the contour 
        p.addParameter('DeltaPhi', 0.5*pi );
        p.addParameter('Evec', [] );
        p.addParameter('Emin', [] );
        p.addParameter('JustScatter', false ); %logical value. True if only an isolated scattering center should be considered in the self-consistent calculations, false otherwise.
        
        p.parse(varargin{:});      
        Edb      = p.Results.Edb;
        Evec     = p.Results.Evec;
        DeltaPhi = p.Results.DeltaPhi;   
        Emin     = p.Results.Emin;
        JustScatter = p.Results.JustScatter;

        obj.create_Hamiltonians( 'junction', obj.junction )
        
        phivec = (0.5 + 0.5*sin( -0.5*pi:pi/(Edb+1):0.5*pi )) * DeltaPhi + (pi-DeltaPhi)/2;
        
        obj.display(['EQuUs:Utils:', class(obj), ':DensityCalc: Calculating the on-site density.'], 1)
        
        if JustScatter && isempty(Emin)
            obj.create_scatter();
            Hscatter = obj.junction.CreateH.Read('Hscatter'); 
            Emin = -abs(eigs( Hscatter, 1, 'lm'))*1.05;
        elseif  isempty(Emin)            
            obj.getBandWidth();        
            Emin = min([-obj.BandWidth.Lead.Emax,-obj.BandWidth.Scatter.Emax])*1.2;
        end
        
        
        if isempty(Evec)
            Emax = obj.junction.getFermiEnergy();
        
            if Emin >= Emax
                density_vec = 0;
                return
            end                   
        
            radius = abs( Emax-Emin) /2;  
            R = radius/sin(DeltaPhi/2);
        
            Evec = R*(cos(phivec) + 1i*sin(phivec)) - radius + Emax - 1i*R*sin((pi-DeltaPhi)/2);
        end
        deltaEvec = [0, diff(Evec)];     
        
        obj.display(['EQuUs:Utils:', class(obj), ':DensityCalc: Calculating the density between Emin = ', num2str(min(Evec)),' and Emax = ',num2str(max(Evec)),],1);        
       
             
        
        %> determine the number of sites in the calculations to preallocate memory
        if JustScatter
            Hscatter = obj.junction.CreateH.Read('Hscatter');        
            number_of_sites = size(Hscatter,1);
        else            
            Hscatter = obj.junction.CreateH.Read('Hscatter');
            Leads = obj.junction.FL_handles.Read('Leads');
            size_K0_leads = zeros(length(Leads), 1);
            size_K0_interface = zeros(length(Leads), 1);
            for kdx = 1:length(Leads)
                Leads{kdx}.Calc_Effective_Hamiltonians( Evec(1) );
                K0_lead = Leads{kdx}.Get_Effective_Hamiltonians();
                size_K0_leads(kdx) = size(K0_lead,1);
            
                obj.junction.CreateInterface( kdx )
                K0_interface = obj.junction.Interface_Regions{kdx}.Get_Effective_Hamiltonians();
                size_K0_interface(kdx) = size(K0_interface,1);
            end
            
            if obj.useSelfEnergy
                number_of_sites = size(Hscatter,1) + sum(size_K0_interface); %2*interface region
            else
                number_of_sites = size(Hscatter,1) + sum(size_K0_interface) + sum(size_K0_leads); %2*lead + 2*interface region
            end
        end        
        
        %> preallocate memory for Densitysurf
        Densitysurf = complex(zeros( length(Evec), length(obj.T), number_of_sites ));
        
        % 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;
        hcomplexDensity = @obj.complexDensity;
        hcreate_scatter = @obj.create_scatter;
        
        junction_loc = obj.junction;
        
        junction_sites = cell(length(Evec),1);
        
        %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, junction_sites_tmp] = feval(hcomplexDensity, junction_loc, JustScatter);
            catch errCause
               feval(hdisplay, ['EQuUs:Utils:', class(obj), ':DensityCalc: Error at energy point E = ', num2str(Evec(iidx)), ' caused by:'], 1);
               feval(hdisplay, errCause.identifier, 1 );
                for jjjdx = 1:length(errCause.stack)
                    feval(hdisplay, ['file: ',errCause.stack(jjjdx).file], 1 );
                    feval(hdisplay, ['name: ',errCause.stack(jjjdx).name], 1 );
                    feval(hdisplay, ['line: ', num2str(errCause.stack(jjjdx).line)], 1 );
                end
               feval(hdisplay, ['EQuUs:Utils:', class(obj), ':DensityCalc: Giving NaN for the calculated current at the given energy point and phase difference.'], 1); 
                densityvec2int = NaN;
            end
                
            Densitysurf( iidx, :, : ) = densityvec2int;
            
            if iidx == 1
               junction_sites{iidx} =  junction_sites_tmp;
            end
            
        end

        %> closing the parallel pool
        poolmanager.closePool();
          
        
          %> calculating the contour integrals
        density_vec = cell( size(obj.T) );
          for kkdx = 1:length( obj.T )
               density_vec{kkdx} = zeros(number_of_sites, 1);         
            for jjdx = 1:number_of_sites
                    densityvec2integrate = Densitysurf(:,kkdx,jjdx);
                    % omitting NaNs by substituting interpolated values
               indexes = isnan(densityvec2integrate); 
                    indexesall = 1:length(densityvec2integrate);
                       
               densityvec2integrate(indexes) = interp1(indexesall(~indexes), densityvec2integrate(~indexes), indexesall(indexes), 'spline' );
                    density = imag(obj.SimphsonInt( densityvec2integrate.*transpose(deltaEvec) ))/(pi);
                    density_vec{kkdx}( jjdx ) = density_vec{kkdx}( jjdx ) + density;
            end
        end    
          if length(obj.T) == 1
               density_vec = density_vec{1};
          end
        
        obj.display(['EQuUs:Utils:', class(obj), ':DensityCalc: Process finished.'], 1)
        
        %reset the system  for a new calculation
        obj.InputParsing( obj.varargin{:} );
        
        % selecting the junction_sites from the cell array
        junction_sites = junction_sites{1};
              
        

  
    end
    

    
%% DensityCalcSmall
%> @brief Calculates the onsite desnity using the method in Nanotechnology 25 (2014), 465201. This  method is usable for small isolated scattering centers for which an exact diagonalization can be performed
%> @param num_occupied_states Number of occupied states
%> @return An array of the calculated density and a structure describing the geometry of the sites.
    function [density_vec, junction_sites, num_occupied_states] = DensityCalcSmall( obj, num_occupied_states )
        
        Max_size = 10000;
        obj.create_scatter();
        Hscatter = obj.junction.CreateH.Read('Hscatter'); 
        coordinates = obj.junction.CreateH.Read('coordinates'); 
        
        if size(Hscatter,1) > Max_size
            error('EQuUs:Utils:Density:DensityCalcSmall', 'Matrix size is too large');          
            density_vec = [];
            junction_sites = [];
            return
        end
        
        
        [eigenvectors, eigenvalues] = eig(full(Hscatter),'nobalance');
        
        eigenvalues = diag(eigenvalues);
        
        [eigenvalues,IX] = sort(eigenvalues, 'ascend');
        
        %eigenvalues = eigenvalues( 1:num_occupied_states );
        
        if ~exist( 'num_occupied_states', 'var') || isempty(num_occupied_states)
            num_occupied_states = find(eigenvalues < obj.junction.EF, 1, 'last');
        end
        
        
        eigenvectors = eigenvectors( :, IX(1:num_occupied_states) );            
        density_vec = sum(transpose(eigenvectors.^2));
        density_vec = reshape( density_vec, numel(density_vec), 1);            
        
        
        %>creating site indexes corresponding to the elements of the density vector
        junction_sites = structures('junction_sites');
        junction_sites.Scatter = structures('sites');
        junction_sites.Scatter.coordinates = coordinates; 
        junction_sites.Scatter.site_indexes = 1:size(Hscatter,1); 
        
        
    end



%% getBandWidth
%> @brief Determines the band width of the leads and of the scattering region.
%> @return ret An instance of structure BandWidth containing the bandwidths.
    function ret = getBandWidth( obj )
     obj.display('Getting Bandwidths', 1);
        
        if ~isempty(obj.BandWidth)
            return
        end
        
        poolmanager = Parallel( obj.junction.FL_handles.Read('Opt') );
        poolmanager.openPool();
        
        Surface_1 = obj.junction.FL_handles.SurfaceGreenFunctionCalculator(1, 'Just_Create_Hamiltonians', true, 'q', obj.junction.q, 'leadmodel', obj.junction.leadmodel, ...
            'CustomHamiltonian', obj.junction.cCustom_Hamiltonians);
        W = BandWidthCalculator( Surface_1 );
        
        param = obj.junction.FL_handles.Read( 'param' );         
        pair_potential = min([abs(param.Leads{1}.pair_potential), abs(param.Leads{2}.pair_potential)]);
        
        obj.BandWidth.Lead = structures('BandWidthLimits');
        obj.BandWidth.Lead.Emin = abs(pair_potential);
        obj.BandWidth.Lead.Emax = W;
        
        
        if strcmpi( class(obj.junction), 'Ribbon' )
            obj.junction.CreateRibbon('justHamiltonians', true);        
            W = BandWidthCalculator( obj.junction.Scatter_UC );
        elseif strcmpi( class(obj.junction), 'NTerminal' )
            obj.junction.CreateScatter();
            W = 1.3*BandWidthCalculatorMolecule( obj.junction.CreateH ); % multiplied by 1.3 just for sure to include normal bound states
        else
            W = 0;
        end
        obj.BandWidth.Scatter = structures('BandWidthLimits');
        obj.BandWidth.Scatter.Emin = 0;
        obj.BandWidth.Scatter.Emax = W;
        
        poolmanager.closePool();

        ret = obj.BandWidth;
         
         %--------------------------------------------
         % nested functions
         function W = BandWidthCalculator( Scatter_UC )

          H0 = Scatter_UC.Read('H0');
            if ~isempty( obj.junction.q )
                H1_transverse = Scatter_UC.Read('H1_transverse');
                H0 = H0 + H1_transverse*exp(1i*obj.junction.q) + H1_transverse'*exp(-1i*obj.junction.q); 
            end
          H1 = Scatter_UC.Read('H1'); 
          u = rand(size(H0,1)); 
          u = u/norm(u);     
            
            hgetMaxEigenvalue = @getMaxEigenvalue;
            hsecular_H = @secular_H;
            k = -pi:2*pi/50:pi;
            arnoldi = zeros( length(k),1);
            parfor (idx = 1:length(k), poolmanager.NumWorkers)
                Heff = feval(hsecular_H,H0,H1,k(idx));
                arnoldi(idx) = abs(feval(hgetMaxEigenvalue, Heff ));                
            end
            W = max( arnoldi );
            

          % --------------------------------------
          % Hamiltonian for the secular equation
          function H = secular_H(H0,H1,k)     
               H = H0 + H1*exp(1i*k) + H1'*exp(-1i*k);    
          end
         
            %--------------------------------------------
          % getting the maximal eigenvalue with Arnoldi iteration
          % Numerical Linear Algebra" by Trefethen and Bau. (p250) or
          % http://en.wikipedia.org/wiki/Arnoldi_iteration
          function ret = getMaxEigenvalue( Heff )
               %u = rand(size(H0,1)); 
               %u = u/norm(u);
               maxIter = 1e2;
               tolerance = 1e-5;
              ret = 0;
               lambda = 1;
               for iidx=1:maxIter
                    u_new = Heff*u;
                    norm_new = norm(u_new);
                    lambda_new = norm_new/norm(u);
               
                    if abs((lambda - lambda_new)/lambda) < tolerance
                         ret = lambda;
                         return
                    end                 
                    u = u_new/norm_new;
                    lambda = lambda_new;          
               end
          end
         end
         
         
         %--------------------------------------------
         function W = BandWidthCalculatorMolecule( CreateH )
            Hscatter  = CreateH.Read('Hscatter');
            W = abs(eigs(Hscatter,1));
         end       
         
         % end nested functions
        
    end
    
    
    end % public methods
    
    
    methods (Access=protected)    
        
%% complexDensity        
%> @brief Calculates the density at a complex energy
%> @param junction_loc An instance of class #NTerminal or its subclass describing the junction.
%> @param JustScatter Logical value. True if only an isolated scattering center should be considered in the self-consistent calculations, false otherwise.
     function [density, junction_sites] = complexDensity( obj, junction_loc, JustScatter )            
        
            
        [spectral_function, junction_sites] = obj.complexSpectral( junction_loc, JustScatter );
             

        fE = obj.Fermi(junction_loc.getEnergy());            
            
          density = zeros(length( fE ), length(spectral_function) );

        for idx = 1:length( fE )
            spectral_function_tmp = spectral_function;
            if junction_loc.Opt.BdG
                coordinates = junction_loc.CreateH.Read('coordinates');
                BdG_indexes = coordinates.BdG_u;
                spectral_function_tmp( BdG_indexes) = spectral_function_tmp( BdG_indexes)*fE(idx);
                spectral_function_tmp( ~BdG_indexes) = spectral_function_tmp( ~BdG_indexes)*(1-fE(idx));
            else
                spectral_function_tmp = spectral_function_tmp*fE(idx);
            end

            density(idx,:) = spectral_function_tmp;

        end
         
    end        



    

%% InputParsing
%> @brief Parses the optional parameters for the class constructor.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'T' The absolute temperature in Kelvins.
%> @param 'silent' Set true to suppress output messages.
%> @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('useSelfEnergy', true);
        p.addParameter('gfininvfromHamiltonian', false);
        p.addParameter('junction', []);     
        
        p.parse( varargin{:});
        
        InputParsing@DOS( 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