Ribbon_Leads.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - Ribbon
%    Copyright (C) 2019 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 unit_tests Unit Tests
%> @{
%> @file Ribbon_Leads.m
%> @brief A class for calculations on a structure made of two connected leds.
%> @Available
%> EQuUs v5.0.144 or later
%> @}
%
%> @file Ribbon_Leads.m
%> @brief A class for calculations on a structure made of two connected leds.
classdef Ribbon_Leads < Ribbon

   
    properties ( Access = public )  
    end
    
    
methods ( Access = public )
%% Contructor of the class
%> @brief Constructor of the class.
%> @param Opt An instance of the structure #Opt.
%> @param varargin Cell array of optional parameters. For details see #InputParsing.
%> @return An instance of the class
    function obj = Ribbon_Leads( varargin )             
        obj = obj@Ribbon();         
        
        
    if strcmpi( class(obj), 'Ribbon_Leads') 
        % processing the optional parameters
        obj.InputParsing(varargin{:})                 

        
        obj.display(['EQuUs:Utils:', class(obj), ':Ribbon_Leads: Creating a Ribbon object'])
        
        % create the shape of the scattering region
        obj.createShape();
        
        % set the Fermi level
        obj.setFermiEnergy();
        
        % exporting calculation parameters into an XML format
        createOutput( obj.filenameOut, obj.Opt, obj.param );
        
        % create class intances and initializing class attributes
        obj.CreateHandles();
        
        % create Hamiltonians and coordinates of the unit cells 
        obj.CreateRibbon();
    end        
        

        
    end 



%% CustomDysonFunc
%> @brief Custom Dyson function for a general two terminal arrangement.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'gfininv' The inverse of the Greens function of the scattering region. For default the inverse of the attribute #G is used.
%> @param 'constant_channels' Logical value. Set true (default) to keep constant the number of the open channels in the leads for each energy value, or false otherwise.
%> @param 'onlyGinverz' Logical value. Set true to calculate only the inverse of the total Green operator, or false (default) to calculate #G as well.
%> @param 'recalculateSurface' A vector of the identification numbers of the lead surfaces to be recalculated.
%> @param 'decimate' Logical value. Set true (default) to eliminate all inner sites in the Greens function and keep only the selected sites. Set false to omit the decimation procedure.
%> @param 'kulso_szabfokok' Array of sites to be kept after the decimation procedure. (Use parameter 'keep_sites' instead)
%> @param 'selfEnergy' Logical value. Set true to use the self energies of the leads in the Dyson equation, or false (default) to use the surface Green function instead.
%> @param 'keep_sites' Name of sites to be kept in the resulted Green function (POssible values are: 'scatter', 'interface', 'lead').
%> @return [1] The calculated Greens function. 
%> @return [2] The inverse of the Green operator. 
%> @return [3] An instance of structure #junction_sites describing the sites in the calculated Green operator.
    function [Gret, Ginverz, junction_sites] = CustomDysonFunc( obj, varargin )
        
    p = inputParser;
    p.addParameter('gfininv', []);
    p.addParameter('constant_channels', true);
    p.addParameter('onlyGinverz', false );
    p.addParameter('recalculateSurface', 1:length(obj.param.Leads) );
    p.addParameter('decimate', true );
    p.addParameter('kulso_szabfokok', []); %The list of sites to be left after the decimation procedure
    p.addParameter('SelfEnergy', false); %set true to calculate the Dyson equation with the self energy
    p.addParameter('keep_sites', 'lead'); %identification of sites to be kept (scatter, interface, lead)
    p.parse(varargin{:});
    gfininv     = p.Results.gfininv;
    constant_channels = p.Results.constant_channels;
    onlyGinverz        = p.Results.onlyGinverz;
    recalculateSurface = p.Results.recalculateSurface;  
    decimate           = p.Results.decimate; 
    kulso_szabfokok = p.Results.kulso_szabfokok;
    useSelfEnergy         = p.Results.SelfEnergy;
    keep_sites        = p.Results.keep_sites;
    
    if isempty(gfininv)
        if ~isempty( obj.Ginv )
            gfininv = obj.Ginv;
        else
            rcond_G = rcond(obj.G);
            if isnan(rcond_G) || abs( rcond_G ) < 1e-15
                gfininv = obj.inv_SVD( obj.G );
            else
                gfininv = inv(obj.G);
            end
        end       
    end

     if ~isempty(recalculateSurface)
    
        % creating interfaces for the Leads
        if constant_channels
            shiftLeads = ones(length(obj.param.Leads),1)*obj.E;
        else
            shiftLeads = ones(length(obj.param.Leads),1)*0;
        end

        % creating objects describing the Leads
        Leads = obj.FL_handles.LeadCalc('shiftLeads', shiftLeads, ...
            'SelfEnergy', useSelfEnergy, 'SurfaceGreensFunction', ~useSelfEnergy, 'gauge_field', obj.gauge_field, 'leads', recalculateSurface, 'q', obj.q, ...
            'leadmodel', obj.leadmodel, 'CustomHamiltonian', obj.cCustom_Hamiltonians);
     else
          Leads = obj.FL_handles.Read( 'Leads' );
     end
    

    Neffs = obj.FL_handles.Get_Neff();   

    
    if useSelfEnergy
        Ginverz = CalcGinverzwithSelfEnergy();
    else
        Ginverz = CalcGinverz();
    end
    
    junction_sites = GetJunctionSites();
    
    if decimate
        GetSitesForDecimation();
        Ginverz = obj.DecimationFunction( kulso_szabfokok, Ginverz);
    end
    
    
    
    if onlyGinverz
        Gret = [];
        return
    end
    
    if issparse( Ginverz )
        err = MException(['EQuUs:Utils:', class(obj), ':CustomDysonFunc', 'Ginverz is sparse. Calculation of the inverse of a sparse matrix is not supported at this point']);
        save('Error_Ribbon_CustomDysonFunc.mat');
        throw(err);
    end
     
    rcond_Ginverz = rcond(Ginverz);
    if isnan( rcond_Ginverz )
          err = MException(['EQuUs:Utils:', class(obj), ':CustomDysonFunc'], 'NaN is the condition number for Ginverz');
        save('Error_Ribbon_CustomDysonFunc.mat');
        throw(err);
    end
    
    if abs( rcond_Ginverz ) < 1e-15
        obj.display( ['EQuUs:Utils:', class(obj), ':CustomDysonFunc: Regularizing Ginverz by SVD'],1);
        Gret = obj.inv_SVD( Ginverz );
    else
        Gret = inv( Ginverz );
    end  
    
    % nested functions
    
        %-------------------------------
        %> calculate the inverse Greens function
        function Ginverz = CalcGinverz()

            Gsurfinverz = cell(length(Neffs),1);
            for idx = 1:length(Neffs)
                Gsurfinverz{idx} = Leads{idx}.Read('gsurfinv');
            end
           
            K_interface = cell(length(Neffs),1);
            K1_sc = cell(length(Neffs),1);
            K1adj_sc = cell(length(Neffs),1);
            K1  = cell(length(Neffs),1);
            K1adj = cell(length(Neffs),1);
            for idx = 1:length(recalculateSurface)
                 obj.CreateInterface( recalculateSurface(idx) );
            end
            for idx = 1:length( obj.Interface_Regions )
                [K_interface{idx}, K1{idx}, K1adj{idx}, K1_sc{idx}, K1adj_sc{idx}] = obj.Interface_Regions{idx}.Get_Effective_Hamiltonians();                      
            end           
            
            % only leads
            %Ginverz = [Gsurfinverz{1}, -K1{1}; -K1adj{1}, Gsurfinverz{2}];
            
            
            % leads and the interface region
            
            % leads
            Ginverz_leads = [Gsurfinverz{1}, zeros(size(Gsurfinverz{1})); zeros(size(Gsurfinverz{1})), Gsurfinverz{2}];
            
            % interface region
            Ginverz_interface = [-K_interface{1}, -K1adj{2};  -K1{2}, -K_interface{2}];
                
            % coupling
            Ginverz_coupling     = [  K1_sc{1}; K1_sc{2}];
            Ginverz_coupling_adj = [  K1adj_sc{1}, K1adj_sc{2}];
            
            Ginverz = [Ginverz_leads, Ginverz_coupling; Ginverz_coupling_adj, Ginverz_interface];
            
        end
        
        
        %-------------------------------
        % calculate the inverse Greens function with the self energy
        function Ginverz = CalcGinverzwithSelfEnergy()
            
            error('Not developed yet')
            
        end  
        
        
        % creating site indexes corresponding to the elements of the calculated Green operator
        function junction_sites = GetJunctionSites()
            
            junction_sites = structures('junction_sites');
            
            % determine the matrix sizes
            size_K0_leads = zeros(length(Leads), 1);
            size_K0_interface = zeros(length(Leads), 1);
            for kdx = 1:length(Leads)
                K0_lead = Leads{kdx}.Get_Effective_Hamiltonians();
                size_K0_leads(kdx) = size(K0_lead,1);
            
                K0_interface = obj.Interface_Regions{kdx}.Get_Effective_Hamiltonians();
                size_K0_interface(kdx) = size(K0_interface,1);
            end
            
            % determine junction sites corresponding to the scattering region
            junction_sites.Scatter = structures('sites');
            junction_sites.Scatter.coordinates = obj.CreateH.Read('coordinates');    
            junction_sites.Scatter.kulso_szabfokok = obj.CreateH.Read('kulso_szabfokok');   
                
            if useSelfEnergy
                junction_sites.Scatter.site_indexes = sum(size_K0_interface) + (1:size(gfininv,1));  %including 2*interface regions            
            else
                junction_sites.Scatter.site_indexes = sum(size_K0_leads) + sum(size_K0_interface) + (1:size(gfininv,1));  %including 2*interface regions and 2 leads             
            end       
            
            % determine junction sites corresponding to the interface
            % regions and leads
            junction_sites.Leads = cell(length(Leads),1);
            junction_sites.Interface = cell(length(Leads),1);
            [~, coordinates_interface] = obj.getCoordinates();
            for kdx = 1:length(Leads)
                junction_sites.Interface{kdx} = structures('sites');
                junction_sites.Leads{kdx} = structures('sites');
            
                if useSelfEnergy                
                    junction_sites.Interface{kdx}.site_indexes = sum(size_K0_interface(1:kdx-1)) + (1:size_K0_interface(kdx));
                    junction_sites.Interface{kdx}.coordinates = coordinates_interface{kdx};
                    
                    junction_sites.Leads{kdx}.site_indexes =  sum(size_K0_interface(1:kdx-1)) + (1:size_K0_leads(kdx));
                    junction_sites.Leads{kdx}.coordinates = Leads{kdx}.Read('coordinates');
                    junction_sites.Leads{kdx}.kulso_szabfokok = Leads{kdx}.Read('kulso_szabfokok');
                else
                    junction_sites.Leads{kdx}.site_indexes =  sum(size_K0_leads(1:kdx-1)) + (1:size_K0_leads(kdx));
                    junction_sites.Leads{kdx}.coordinates = Leads{kdx}.Read('coordinates');
                    junction_sites.Leads{kdx}.kulso_szabfokok = Leads{kdx}.Read('kulso_szabfokok');
                    
                    junction_sites.Interface{kdx}.site_indexes = sum(size_K0_leads) + sum(size_K0_interface(1:kdx-1)) + (1:size_K0_interface(kdx));
                    junction_sites.Interface{kdx}.coordinates = coordinates_interface{kdx};
                end
            
                
               
            end
            
        end
        
        
        % Obtain the site indexes to be decimated out and removes the unnecesarry sites from the structure junction_sites according to the
        % performed decimation
        function GetSitesForDecimation()
                        
            if isempty(kulso_szabfokok)
                if strcmpi(keep_sites, 'scatter')
                    kulso_szabfokok = get_scatter_sites();
                             
                    junction_sites.Scatter.site_indexes = 1:length(kulso_szabfokok);
                    junction_sites.Interface = [];
                    junction_sites.Leads = [];
                    
                elseif strcmpi(keep_sites, 'interface')
                    kulso_szabfokok = get_interface_sites();
                    
                    junction_sites.Scatter = [];
                    interface_size = 0;
                    for idx = 1:length( junction_sites.Interface )
                        junction_sites.Interface{idx}.site_indexes = interface_size + (1: length(junction_sites.Interface{idx}.site_indexes));
                        interface_size = interface_size + length(junction_sites.Interface{idx}.site_indexes);
                    end
                    junction_sites.Leads = [];
                    
                elseif strcmpi(keep_sites, 'lead')
                    kulso_szabfokok = get_lead_sites();
                    
                    junction_sites.Scatter = [];
                    junction_sites.Interface = [];
                    leads_size = 0;
                    for idx = 1:length( junction_sites.Leads )
                        junction_sites.Leads{idx}.site_indexes = leads_size + (1: length(junction_sites.Leads{idx}.site_indexes));
                        leads_size = leads_size + length(junction_sites.Leads{idx}.site_indexes);
                    end
                else
                    error(['EQuUs:Utils:', class(obj), ':CustomDysonFunc'], ['Undifined Sites: ', keep_sites]);
                end                
                
            else
                % for custom given kulso_szabfokok
                
                % determine sites in the scattering center
                junction_sites.Scatter = SelectSites( junction_sites.Scatter );
                
                
                % determine sites in the interface regions
                for idx = 1:length(junction_sites.Interface)
                    junction_sites.Interface{idx} = SelectSites( junction_sites.Interface{idx} );
                end
                
                
                % determine sites in the interface regions
                for idx = 1:length(junction_sites.Leads)
                    junction_sites.Leads{idx} = SelectSites( junction_sites.Leads{idx} );
                end
                            
                
            end
            
            %-------------------------------------------------------------
            % an auxilary function to select the sites to be kept
            function sites_ret = SelectSites( sites )
                site_indexes_tmp = sites.site_indexes;
                start_index = find( min(site_indexes_tmp) <= kulso_szabfokok, 1, 'first');
                end_index = find( max(site_indexes_tmp) >= kulso_szabfokok, 1, 'last');
                
                if isempty(start_index)
                    sites_ret = [];
                    return
                end
                
                sites_ret = structures('sites');
                sites_ret.coordinates = sites.coordinates;
                
                % selecting site indexes
                if isempty(end_index)
                    sites_ret.site_indexes = 1:length(kulso_szabfokok(start_index:end));
                else
                    sites_ret.site_indexes = 1:length(kulso_szabfokok(start_index:end_index));
                end
                        
                % selecting coordinate sof the sites
                names = fieldnames( sites.coordinates );
                for iidx = 1:length(names)
                    fieldname = names(iidx);
                    if strcmpi( fieldname, 'a') || strcmpi( fieldname, 'b')
                        continue
                    end
                        
                    field_tmp = sites_ret.coordinates.(fieldname);
                        
                    if isempty(field_tmp)
                        continue
                    end
                        
                    if isempty(end_index)
                        field_tmp = field_tmp( kulso_szabfokok(start_index:end) - min(site_indexes_tmp) + 1 );
                    else
                        field_tmp = field_tmp( kulso_szabfokok(start_index:end_index) - min(site_indexes_tmp) + 1 );
                    end
                    sites_ret.coordinates.(fieldname) = field_tmp;
                end         
                
                
            end
            
        end
        
        
        %-------------------------------
        % determine the site indexes of the scattering region within Ginverz
        function kulso_szabfokok = get_scatter_sites()
            kulso_szabfokok = junction_sites.Scatter.site_indexes;
        end
        
        %-------------------------------
        % determine the site indexes of the interface region within Ginverz
        function kulso_szabfokok = get_interface_sites()
            size_interface = 0;
            for idx = 1:length( junction_sites.Interface )
                    size_interface = size_interface + length(junction_sites.Interface{idx}.site_indexes);
            end
            
            kulso_szabfokok = zeros( size_interface, 1);
            size_interface = 0;
            for idx = 1:length( junction_sites.Interface )
                indexes = size_interface + ( 1:length(junction_sites.Interface{idx}.site_indexes) );
                kulso_szabfokok(indexes) = junction_sites.Interface{idx}.site_indexes;
                size_interface = size_interface + length(junction_sites.Interface{idx}.site_indexes);
            end     
                
        end        
        
        %-------------------------------
        % determine the site indexes of the leads within Ginverz
        function kulso_szabfokok = get_lead_sites()
                        
                size_leads = 0;
                for idx = 1:length( junction_sites.Leads )
                    size_leads = size_leads + length(junction_sites.Leads{idx}.site_indexes);
                end
            
                kulso_szabfokok = zeros( size_leads, 1);
                size_leads = 0;
                for idx = 1:length( junction_sites.Leads )
                    indexes = size_leads + ( 1:length(junction_sites.Leads{idx}.site_indexes) );
                    kulso_szabfokok(indexes) = junction_sites.Leads{idx}.site_indexes;
                    size_leads = size_leads + length(junction_sites.Leads{idx}.site_indexes);
                end                                 
                 

        end
        
        % end nested functions
        
    end



%% CreateClone
%> @brief Creates a clone of the present object.
%> @return Returns with the cloned object.
    function ret = CreateClone( obj )
        
        ret = Ribbon_Leads( 'width', obj.width, ...
            'height', obj.height, ...
            'filenameIn', obj.filenameIn, ...
            'filenameOut', obj.filenameOut, ...
            'E', obj.E, ...
            'EF', 0, ...
            'phi', obj.phi, ...
            'silent', obj.silent, ...
            'transversepotential', obj.transversepotential, ...
               'Opt', obj.Opt, ...
               'param', obj.param, ...
            'q', obj.q, ...
            'leadmodel', obj.leadmodel, ...
            'interfacemodel', obj.interfacemodel);
        
        ret.EF = obj.EF;
        ret.CreateH = obj.CreateH.CreateClone();
        ret.FL_handles = obj.FL_handles.CreateClone();
        ret.Scatter_UC = obj.Scatter_UC.CreateClone();
        ret.Interface_Regions = cell(size(obj.Interface_Regions));
        for idx = 1:length(obj.Interface_Regions)
            ret.Interface_Regions{idx} = obj.Interface_Regions{idx}.CreateClone();
        end
        if ~isempty( obj.PeierlsTransform_Scatter )
            ret.PeierlsTransform_Scatter = obj.PeierlsTransform_Scatter.CreateClone();
        end
        
        if ~isempty( obj.PeierlsTransform_Leads )
            ret.PeierlsTransform_Leads   = obj.PeierlsTransform_Leads.CreateClone();    
        end
        ret.gauge_field              = obj.gauge_field;  % function handle for the scalar field to transform the vector potential from Landauy to Landaux
        
    end  


end % methods public

    
methods (Access=protected)        


%% 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 'width' Integer. The number of the nonsingular atomic sites in the cross section of the ribbon.
%> @param 'height' Integer. The height of the ribbon in units of the lattice vector.
%> @param 'filenameIn' The input filename containing the computational parameters. (Use parameters 'Op' and 'param' instead)
%> @param 'filenameOut' The output filename to export the computational parameters.
%> @param 'WorkingDir' The absolute path to the working directoy.
%> @param 'E' The energy value used in the calculations (in the same units as the Hamiltonian).
%> @param 'EF' The Fermi energy in the same units as the Hamiltonian. Attribute #E is measured from this value. (Use for equilibrium calculations in the zero temperature limit. Overrides the one comming from the external source)
%> @param 'silent' Set true to suppress output messages.
%> @param 'transversepotential' A function handle pot = f( #Coordinates ) or pot=f( #CreateLeadHamiltonians, Energy) of the transverse potential applied in the lead. (Instead of #CreateLeadHamiltonians can be used its derived class)
%> @param 'leadmodel' A function handle #Lead=f( idx, E, varargin ) of the alternative lead model with equivalent inputs and return values as #Transport_Interface.SurfaceGreenFunctionCalculator and with E standing for the energy.
%> @param 'interfacemodel' A function handle f( #InterfaceRegion ) to manually adjus the interface regions. (Usefull when 'leadmodel' is also given. For example see @InterfaceModel)
%> @param 'Opt' An instance of the structure #Opt.
%> @param 'param' An instance of the structure #param.
%> @param 'q' The transverse momentum quantum number.
    function InputParsing( obj, varargin )
        % calling the input parsing of the superclass
        InputParsing@Ribbon( obj, varargin{:});  
        
    end
    
end % methdos protected



end