EigenProblemLead.m

Go to the documentation of this file.

%%   Eotvos Quantum Transport Utilities - EigenProblemLead
%    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 basic Basic Functionalities
%> @{
%> @file EigenProblemLead.m
%> @brief Class to solve the eigenproblem of a translational invariant leads and calculate the group velocities.
%> @} 
%> @brief Class to solve the eigenproblem of a translational invariant leads and calculate the group velocities.
%> @Available
%> EQuUs v4.8 or later
%%
classdef EigenProblemLead < handle & SVDregularizationLead
    
    
properties ( Access = protected )
%> true for calculating the retarded Green function, or false for the advanced Green function.
    retarted
%> The unsorted right-sided eigenstates.
    modusmtx
%> The unsorted left-sided eigenstates.
    modusmtx_left
%> The unsorted wave numbers in form exp(1i*k).
    expk
%> The unsorted group velocities.
    csoportseb
%> The wave numbers of the eigenstates, that propagates or decays in the positive direction.
    expk_p
%> The wave numbers of the eigenstates in form exp(1i*k), that propagates or decays in the negative direction.
    expk_m
%> The group velocities of the eigenstates in form exp(1i*k), that propagates or decays in the positive direction.
    vcsop_p
%> The group velocities of the eigenstates, that propagates or decays in the negative direction.
    vcsop_m
%> The right sided wave functions of the eigenstates, that propagates or decays in the positive direction.
    modusmtx_p
%> The right-sided wave functions of the eigenstates, that propagates or decays in the negative direction.
    modusmtx_m    
%> The left-sided wave numbers of the eigenstates, that propagates or decays in the positive direction.
    modusmtx_p_left
%> The left-sided wave functions of the eigenstates, that propagates or decays in the negative direction.
    modusmtx_m_left    
%> The dual basis of the right-sided wave functions of the eigenstates, that propagates or decays in the positive direction.
    d_modusmtx_p    
%> The dual basis of the right-sided wave functions of the eigenstates, that propagates or decays in the negative direction.
    d_modusmtx_m
%> A real number corresponding to the tolerance used to sort the left and right moving (decaying) modes.
    sort_tolerance
%> A structure #open_channels containing info on the open channels
    open_channels   
%> Logical matrix containing the degenerate k subspaces
    degenerate_k_subspaces
%> The energy value for which the #TrukkosSajatertekek eigenvalue problem was solved.
    E
end




methods ( Access = public )
%% constructorof the class
%> @brief Constructor of the class.
%> @param Opt An instance of the structure #Opt.
%> @param param An instance of structure #param.
%> @param varargin Cell array of optional parameters identical to #SVDregularizationLead.SVDregularizationLead.
%> @return An instance of the class
    function obj = EigenProblemLead(Opt, param, varargin)    
        obj = obj@SVDregularizationLead( Opt, param, varargin{:} ); 
    end

%% Determine_Open_Channels
%> @brief Determines the open channels in the lead. The data are storen within the attribute #open_channels.
    function Determine_Open_Channels( obj )
        if isempty(obj.expk_p)
            obj.display('Right and left going states were not determined yet. Use method EigenProblemLead::szetvalaszto first.', 1);
            return;
        end
        propagating_indexes_p = logical( abs( abs(obj.expk_p) - 1 ) < obj.sort_tolerance );
        propagating_indexes_m = logical( abs( abs(obj.expk_m) - 1 ) < obj.sort_tolerance );
        
        obj.open_channels = structures('open_channels');
        obj.open_channels.num = length(find(propagating_indexes_p));  
        obj.open_channels.open_channels_p = propagating_indexes_p;
        obj.open_channels.open_channels_m = propagating_indexes_m;
        
        coordinates = obj.Get_Effective_Coordinates();
        
        if obj.Opt.BdG
            if obj.isSuperconducting()
                obj.open_channels.isSuperconducting = true;
            else
                obj.open_channels.BdG_u_incoming = logical(diag(obj.modusmtx_p(coordinates.BdG_u,:)'*obj.modusmtx_p(coordinates.BdG_u,:)));
                obj.open_channels.BdG_u_outgoing = logical(diag(obj.modusmtx_m(coordinates.BdG_u,:)'*obj.modusmtx_m(coordinates.BdG_u,:)));
                obj.open_channels.isSuperconducting = false;
            end
        else
            obj.open_channels.isSuperconducting = false;
        end
        
    end
              
%% szetvalaszto
%> @brief Sorts the left and right propagating (decaying) modes.
%> @param tolerance The tolerance to be used during the sort.
    function szetvalaszto( obj, tolerance )
        
        if isempty(obj.csoportseb)
            return
        end
        
            
        if ~exist('tolerance', 'var')
            tolerance  = 1e-6;
        elseif tolerance < 1e-16
            err = MException(['EQuUs:', class(obj), ':szetvalaszto'], 'Failed to separate the left and right-going states');
            save(['EQuUs_', class(obj), '_szetvalaszto.mat']);
            throw(err); 
        end
            
            if obj.retarted
                indexes = (real(obj.csoportseb*obj.Lead_Orientation) > 0 & abs(abs(obj.expk)-1) < tolerance) | (1-abs(obj.expk))*obj.Lead_Orientation > tolerance;
            else
                indexes = (real(obj.csoportseb*obj.Lead_Orientation) < 0 & abs(abs(obj.expk)-1) < tolerance) | (1-abs(obj.expk))*obj.Lead_Orientation > tolerance;
            end
            obj.expk_p          = obj.expk(indexes);
            obj.vcsop_p         = obj.csoportseb(indexes);
            obj.modusmtx_p      = obj.modusmtx(:,indexes);
            if ~obj.Opt.usingDualModes
                obj.modusmtx_p_left = obj.modusmtx_left(indexes,:);
            end
            
            obj.expk_m          = obj.expk(~indexes);
            obj.vcsop_m         = obj.csoportseb(~indexes); 
            obj.modusmtx_m      = obj.modusmtx(:,~indexes);
            if ~obj.Opt.usingDualModes
                obj.modusmtx_m_left = obj.modusmtx_left(~indexes,:);
            end
            
            obj.sort_tolerance  = tolerance;
            if length(obj.expk_p) ~= length(obj.expk_m)
                obj.szetvalaszto( obj.sort_tolerance/10 );
            else
                obj.modusmtx   = [];
                obj.csoportseb = [];
                obj.expk          = []; 
            end
                          
        
        
    end

%% Group_Velocity
%> @brief Calculates the group velocities corresponding to the propagating states. The calculated group velocities are stored within the class.
    function Group_Velocity(obj)
        
        usingDualModes = obj.Opt.usingDualModes;        
        
        csoportseb = obj.csoportseb;
        
        needs_to_be_diagonalized = false;
        
        GroupVelocityCalc();
        
        if needs_to_be_diagonalized
            Diagonalize_Current_Operator();
            GroupVelocityCalc();
        end
        
        obj.csoportseb = diag(csoportseb);
        
        
        %------------------------------------------------------------------
        % nested functions
        function Diagonalize_Current_Operator()
            
            [U, eigvals] = eig( full(csoportseb));
            
            % rotate the wavefunctions
            obj.modusmtx = obj.modusmtx*U;
            if ~usingDualModes  
                obj.modusmtx_left = U'*obj.modusmtx_left;
            end            
            obj.NormalizeModes()
            
            % reordering the wavenumbers
            obj.expk = diag(U'*diag(obj.expk)*U);
            
            obj.degenerate_k_subspaces = logical(speye(length(obj.expk)));
            
        end
        %------------------------------------------------------------------
        
        %------------------------------------------------------------------
        function GroupVelocityCalc()
            
            [~, K1_eff, K1adj_eff] = obj.Get_Effective_Hamiltonians();
            Neff = obj.Get_Neff();
            csoportseb = sparse([],[],[], 2*Neff, 2*Neff);
            if usingDualModes   
                adj_modusmtx = obj.modusmtx';
                for idx=1:2*Neff
                    jdx_indexes = obj.degenerate_k_subspaces(idx,:);
                    if ~needs_to_be_diagonalized && length(csoportseb(idx, jdx_indexes)) >1
                        needs_to_be_diagonalized = true;
                    end
                    %Eq (18) in GNew J. Phys 093029 (2014), (modes are already normalized with the overlap integrals)
                    csoportseb(idx, jdx_indexes) = 1i*adj_modusmtx(idx,:)*(K1_eff*obj.expk(idx)-K1adj_eff*(1/obj.expk(idx)) )*obj.modusmtx(:,jdx_indexes);
                end
                
            else
                %equation (A4) in PRB 78, 035407 (difference in the evanescent group velocities)
                 for idx=1:2*Neff
                    jdx_indexes = obj.degenerate_k_subspaces(idx,:);
                    if ~needs_to_be_diagonalized && length(csoportseb(idx, jdx_indexes)) >1
                        needs_to_be_diagonalized = true;
                    end                    
                    csoportseb(idx, jdx_indexes) = 1i*obj.modusmtx_left(idx,:)*(K1_eff*obj.expk(idx)-K1adj_eff*(1/obj.expk(idx)) )*obj.modusmtx(:,jdx_indexes);
                end               
            end            
        end
        % end nested functions
        %------------------------------------------------------------------
        
        
       
    end
    

%% TrukkosSajatertekek
%> @brief Calculates the wave numbers corresponding to the propagating states at given energy. The calculated wave numbers are stored by attribute #expk.
%> @param E The energy value used in the calculations.
    function TrukkosSajatertekek(obj, E)
    
        if imag(E) < 0
            obj.retarted = false;
        else
            obj.retarted = true;
        end 
        
        obj.E = E;
        
        % create effective Hamiltonians with SVD regularization if necessary        
        obj.Calc_Effective_Hamiltonians( E );
        [K0_eff, K1_eff, K1adj_eff] = obj.Get_Effective_Hamiltonians();
        Neff = obj.Get_Neff();
        
        usingDualModes = obj.Opt.usingDualModes;
        
        if obj.isSuperconducting() || ~obj.Opt.BdG
            %  if normal system, or if BdG Hamiltonian with nonzero pair potential
            [obj.modusmtx,obj.expk, obj.modusmtx_left] = calcTrickyEigenvalues( K0_eff, K1_eff, K1adj_eff );
        else
            % if BdG Hamiltonian with nonzero pair potential
            M0 = Neff/2;
            K0_e = K0_eff(1:M0, 1:M0);
            K0_h = K0_eff(M0+1:Neff, M0+1:Neff);
            K1_e = K1_eff(1:M0, 1:M0);
            K1_h = K1_eff(M0+1:Neff, M0+1:Neff);
            K1adj_e = K1adj_eff(1:M0, 1:M0);
            K1adj_h = K1adj_eff(M0+1:Neff, M0+1:Neff);
                
            % eigenvalues for electron-like part
            [modusmtx_e, expk_e, modusmtx_left_e] = calcTrickyEigenvalues( K0_e, K1_e, K1adj_e );
                
            %eigenvalues for hole-like part
            [modusmtx_h, expk_h, modusmtx_left_h] = calcTrickyEigenvalues( K0_h, K1_h, K1adj_h );
            
            obj.modusmtx = [ modusmtx_e, zeros( size(modusmtx_e,1), size(modusmtx_e,2)); ...
                        zeros(size(modusmtx_h,1), size(modusmtx_e,2)), modusmtx_h ];
                    
            obj.expk = [expk_e;expk_h];            
            
            if ~usingDualModes
                obj.modusmtx_left = [ modusmtx_left_e, zeros( size(modusmtx_left_e,1), size(modusmtx_left_e,2)); ...
                            zeros(size(modusmtx_left_h,1), size(modusmtx_left_e,2)), modusmtx_left_h ];
            end
                
        end 
        
        obj.NormalizeModes();   
        obj.Determine_Degenerate_k_subspaces();

        obj.modusmtx_p     = [];
        obj.modusmtx_m     = [];
        obj.d_modusmtx_p     = [];
        obj.d_modusmtx_m     = [];
        
        
      
        
        
        %------------------------------------------------------------------
        % nested functions
        % calculates the tricky eigenvalue problem
        function [modusmtx,expk, modusmtx_left] = calcTrickyEigenvalues( H0, H1, H1adj )
            M = size(H0, 1);
            [H_eff,B_eff] = Effective_H(H0,H1,H1adj);
            
            if usingDualModes                
                [modusmtx,expk]  = eig(H_eff,B_eff);
                expk             = diag(expk); 
                modusmtx_left = []; 
            else
                [modusmtx, modusmtx_left, expk] = obj.eig(H_eff,B_eff);
            end
            
            if ~usingDualModes   
                % Eq (A3) in PRB 78, 035407
                %normfactors = diag(modusmtx_left*B_eff*modusmtx);
                %modusmtx_left  = diag(1./normfactors)*modusmtx_left;
                modusmtx_left = modusmtx_left(:,1:M);
                modusmtx_left  = diag(sqrt(expk))*modusmtx_left;
            end 
                 
        
            % Eq (7) in PRB 78, 035407
            modusmtx       = modusmtx(1:M,:);   
            modusmtx       = modusmtx*diag(1./sqrt(expk));
        end
        %------------------------------------------------------------------
                
                
        %------------------------------------------------------------------
        % Eq (6) PRB 78 035407
        function [H_eff,B_eff] = Effective_H(H0,H1,H1adj)
            M = size(H0,1);
            %H1inv = H1\eye(size(H1));%inv(H1);
            H_eff    = full([-H0,full(-H1adj);eye(M),zeros(M)]);  
            B_eff    = full([ H1, zeros(M); zeros(M), eye(M)]);
            
        end
        %------------------------------------------------------------------
        % end nested functions
    end


%% AddPotential
%> @brief Adds on-site potential to the Hamiltonian H0.
%> @param V The potential calculated on the sites.
    function AddPotential( obj, V )        

        AddPotential@CreateLeadHamiltonians( obj, V );
    end
    
%% Unitary_Transform
%> @brief Transforms the effective Hamiltonians by a unitary transformation
%> @param Umtx The matrix of the unitary transformation.
    function Unitary_Transform(obj, Umtx)
        
        Unitary_Transform@SVDregularizationLead( obj, Umtx );
        
        if ~isempty(obj.modusmtx_p)
            obj.modusmtx_p = Umtx*obj.modusmtx_p;
        end
        
        if ~isempty(obj.modusmtx_m)
            obj.modusmtx_m = Umtx*obj.modusmtx_m;
        end
        
        if ~isempty(obj.modusmtx_p_left)
            obj.modusmtx_p_left = Umtx*obj.modusmtx_p_left;
        end
        
        if ~isempty(obj.modusmtx_m_left)
            obj.modusmtx_m_left = Umtx*obj.modusmtx_m_left;
        end
        
        if ~isempty(obj.d_modusmtx_p)
            obj.d_modusmtx_p = obj.modusmtx_p*Umtx';
        end        
        
        if ~isempty(obj.d_modusmtx_m)
            obj.d_modusmtx_m = obj.modusmtx_m*Umtx';
        end
                 
        
    end        

%% CreateClone
%> @brief Creates a clone of the present object.
%> @return Returns with the cloned class.
%> @param varargin Cell array of optional parameters:
%> @param 'empty' Set true to create an empty class, or false (default) to copy all the attributes.
%> @return Returns with an instance of the cloned class.
    function ret = CreateClone( obj )
        
        p = inputParser;
        p.addParameter('empty', false );
        
        p.parse(varargin{:});       
        empty    = p.Results.empty;
        
        ret = EigenProblemLead(obj.Opt, obj.param,...
                                            'Hanyadik_Lead', obj.Hanyadik_Lead,...
                                            'Lead_Orientation', obj.Lead_Orientation, ...
                                            'q', obj.q);   
        if empty
            return
        end                                        
                                        
        meta_data = metaclass(obj);
        
        for idx = 1:length(meta_data.PropertyList)
            prop_name = meta_data.PropertyList(idx).Name;
            if strcmpi( prop_name, 'next_SVD_cycle' )
              if ~isempty( obj.next_SVD_cycle )
               ret.next_SVD_cycle = obj.next_SVD_cycle.CreateClone();
              end
            else
                 ret.Write( prop_name, obj.(prop_name));  
            end
        end   
        
    end   

%% Reset
%> @brief Resets all elements in the object.
    function Reset( obj )
        
        if strcmpi( class(obj), 'EigenProblemLead' )
            meta_data = metaclass(obj);
        
            for idx = 1:length(meta_data.PropertyList)
                prop_name = meta_data.PropertyList(idx).Name;
                if strcmp(prop_name, 'Opt') || strcmp( prop_name, 'param') || strcmp(prop_name, 'varargin')
                    continue
                end                
                obj.Clear( prop_name ); 
            end   
        end    
        
        obj.Initialize(); 
             
    end

%% Write
%> @brief Sets the value of an attribute in the class.
%> @param MemberName The name of the attribute to be set.
%> @param input The value to be set.
    function Write(obj, MemberName, input)       
        
        if strcmp(MemberName, 'param') || strcmp(MemberName, 'params')
            Write@CreateLeadHamiltonians( obj, MemberName, input );
            return        
        else
            try
               obj.(MemberName) = input;
               catch
                    error(['EQuUs:', class(obj), ':Read'], ['No property to write with name: ',  MemberName]);
               end
        end       
        
    end
%% Read
%> @brief Query for the value of an attribute in the class.
%> @param MemberName The name of the attribute to be set.
%> @return Returns with the value of the attribute.
    function ret = Read(obj, MemberName)
        
        if strcmp(MemberName, 'Hamiltonian2Dec')
            ret = Read@CreateLeadHamiltonians( obj, MemberName );
            return            
        else
            try
               ret = obj.(MemberName);
               catch
                    error(['EQuUs:', class(obj), ':Read'], ['No property to read with name: ',  MemberName]);
               end   
        end      
        
    end
%% Clear
%> @brief Clears the value of an attribute in the class.
%> @param MemberName The name of the attribute to be cleared.
    function Clear(obj, MemberName)
        
        try
          obj.(MemberName) = [];
          catch
               error(['EQuUs:', class(obj), ':Clear'], ['No property to clear with name: ',  MemberName]);
          end  
        
    end

end %methods public



methods (Access=private)
    
%% Determine_Degenerate_k_subspaces
%> @brief Determine the degenerate k subspaces for which the current operator needs to be diagonalized.
    function Determine_Degenerate_k_subspaces(obj) 
        
        degeneracy_tolerance = 1e-12;
        
        obj.degenerate_k_subspaces = obj.expk*transpose(1./obj.expk);
        obj.degenerate_k_subspaces = abs(obj.degenerate_k_subspaces - 1) < degeneracy_tolerance; 
 
    end


%% NormalizeModes
%> @brief Renormalizes the wave functions with the overlap integrals (Eq (9) in PRB 78 035407
    function NormalizeModes(obj) 
        usingDualModes = obj.Opt.usingDualModes;
            
        modusmtx_adj = obj.modusmtx';
        [S0, S1, S1adj] = obj.Get_Effective_Overlaps();
        if ~isempty(S0)
            for idx = 1:length(obj.expk)
                expk = obj.expk(idx);
                if length(obj.q) < 2 && ~isempty(obj.q) && ~obj.HamiltoniansDecimated
                    S0 = S0 + obj.S1_transverse*diag(exp(-1i*obj.q)) + obj.S1_transverse'*diag(exp(1i*obj.q)); 
                else
                    S0 = S0;
                end
                S = S0 + S1*expk + S1adj*(1/expk);
                
                l_n = modusmtx_adj(idx,:)*S*obj.modusmtx(:,idx);
                obj.modusmtx(:,idx) = obj.modusmtx(:,idx)/sqrt(l_n);
                    
                % normalize the left sided eigenvectors
                if ~usingDualModes
                    l_n = obj.modusmtx_left(idx,:)*S*obj.modusmtx(:,idx);
                    obj.modusmtx_left(idx,:) = obj.modusmtx_left(idx,:)/(l_n);
                end
            end
        else
            %no overlap integrals are present
            normfactors = diag( modusmtx_adj*obj.modusmtx );
            obj.modusmtx = obj.modusmtx*diag(1./sqrt(normfactors));
            if ~usingDualModes
                normfactors = diag( obj.modusmtx_left*obj.modusmtx );
                obj.modusmtx_left = diag(1./(normfactors))*obj.modusmtx_left;
            end
        end
    end            
    
    
%% Extend_Wavefncs --- DEVELOPMENT
%> @brief Extend the reduced wave functions to the original basis before the SVD regularization  (Eq (45) in PRB 78 035407) This function is in a development stage.
    function Extend_Wavefncs(obj) %NEW   
        wavefnc_extended = zeros(obj.M, obj.Get_Neff());
        for idx = 1:length(obj.expk)
            wavefnc_extended(:,idx) = obj.Extend_Wavefnc( obj.modusmtx(:,idx), obj.expk(idx) );
        end
        
        obj.modusmtx = wavefnc_extended;
    end
        



end % methods private


methods (Access=protected)

%% Initialize
%> @brief Initializes object properties.
    function Initialize(obj)
        Initialize@SVDregularizationLead(obj);
        obj.sort_tolerance = 1e-6;    
    end
    
end
    
    
    
end