SVDregularizationLead.m

Go to the documentation of this file.

%%   Eotvos Quantum Transport Utilities - SVDregularizationLead
%    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 SVDregularizationLead.m
%> @brief Class to regulerize the H1 problem of the leads by SVD decompozition.
%> @} 
%> @brief Class to regulerize the H1 problem of the leads by SVD decompozition.
%> The notations and the structure of the Hamiltonian is defined accroding to the following image:
%> @image html Lead_Hamiltonian.jpg
%> @image latex Lead_Hamiltonian.jpg
%> @Available
%> EQuUs v4.8 or later
%%
classdef SVDregularizationLead < handle & CreateLeadHamiltonians
    
    
properties ( Access = protected )
    %> true if the Hamiltonians were SVD transformed, false otherwise
    is_SVD_transformed
    %> U matrix from the SVD decompozition, see Eq  (41) of PRB 78, 035407 
    U
    %> S matrix from the SVD decompozition, see Eq  (41) of PRB 78, 035407 
    S
    %> V matrix from the SVD decompozition, see Eq  (41) of PRB 78, 035407 
    V
    %> SVD tolerance to identify singular values
    tolerance
    %> Somethimes it is needed to perform another SVD cycle to regularize the H1 matrix
    next_SVD_cycle
    %> Effective number of sites after the elimination of the singular values
    Neff    


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 #CreateLeadHamiltonians.CreateLeadHamiltonians.
%> @return An instance of the class
    function obj = SVDregularizationLead(Opt, param, varargin)    
        obj = obj@CreateLeadHamiltonians( Opt, param, varargin{:} );  
    end

    
    
%% SVDdecompozition
%> @brief Calculates the SVD decomposition of the matrix H1
%> @param The transverse momentum dependent coulping between the unit cell (K1 is transverse momentum dependent, only if skew coupling #H1_skew_right or #H1_skew_left and the transverse momentum #q are not empty).
    function SVDdecompozition( obj, K1 )   
        if issparse(K1)
            [obj.U, obj.S, obj.V] = svd(full(K1));
          else
               [obj.U, obj.S, obj.V] = svd(K1);
        end  
        
        obj.U = sparse(obj.U);
        obj.V = sparse(obj.V);
        
        obj.S = diag(obj.S);
    end
    
%% is_SVD_needed
%> @brief Decides whether SVD regularization is needed or not.
%> @return Returns with true if SVD regularization is needed, false otherwise
    function ret = is_SVD_needed(obj)        
        if 1/condest(obj.H1)<obj.tolerance
            ret = true;
        else
            ret = false;
        end
    end
    
    
%% Calc_Effective_Hamiltonians
%> @brief Calculates the effective Hamiltonians according to Eq (48) of of PRB 78, 035407 
%> @param E The energy value
    function Calc_Effective_Hamiltonians( obj, E ) 
        
        obj.ApplyOverlapMatrices(E);
        
        if ~isempty(obj.kulso_szabfokok) && length(obj.kulso_szabfokok) ~= size(obj.K0,1)
            obj.Decimate_Hamiltonians();
            return
        elseif 1/condest(obj.K1)<obj.tolerance && isempty(obj.V)
            obj.SVD_transform();
            return
        else
            % in case no SVD regularization or simple decimation is needed 
            obj.Neff = size(obj.K0,1);
            return
        end
        
               
    end
    

%% Get_Effective_Hamiltonians
%> @brief Gets the effective Hamiltonians K0_eff, K1_eff, K1adj_eff according to Eq (48) of of PRB 78, 035407 
%> @return [1] The effective Hamiltonian of one unit cell
%> @return [2] The effective coupling between unit cells K1_eff
%> @return [3] The adjungate of the effective coupling between unit cells K1adj_eff 
    function [K0_eff, K1_eff, K1adj_eff] = Get_Effective_Hamiltonians(obj)
        if isempty(obj.next_SVD_cycle)
            [K0_eff, K1_eff, K1adj_eff] = obj.qDependentHamiltonians();
        elseif strcmpi( class(obj.next_SVD_cycle), 'SVDregularizationLead' )
            [K0_eff, K1_eff, K1adj_eff] = obj.next_SVD_cycle.Get_Effective_Hamiltonians();
        else
            error('EQuUs:SVDregularizationLead:Get_Effective_Hamiltonians', 'Unrecognized type of attribute next_SVD_cycle');
        end
        
    end
    
%% Get_Effective_Overlaps
%> @brief Gets the effective Hamiltonians S0_eff, S1_eff, S1adj_eff according to Eq (48) of of PRB 78, 035407 
%> @return [1] The effective overlap matrix of one unit cell S0_eff
%> @return [2] The effective overlap matrix between unit cells S1_eff
%> @return [3] The adjungate of the effective overlap matrix between unit cells S1adj_eff 
    function [S0_eff, S1_eff, S1adj_eff] = Get_Effective_Overlaps(obj)
        if isempty(obj.next_SVD_cycle)
            [S0_eff, S1_eff, S1adj_eff] = obj.qDependentOverlaps();
        elseif strcmpi( class(obj.next_SVD_cycle), 'SVDregularizationLead' )
            [S0_eff, S1_eff, S1adj_eff] = obj.next_SVD_cycle.Get_Effective_Overlaps();
        else
            error('EQuUs:SVDregularizationLead:Get_Effective_Overlaps', 'Unrecognized type of attribute next_SVD_cycle');
        end
        
    end    
    
%% Get_Effective_Coordinates
%> @brief Gets the coordinates of the sites of the effective Hamiltonians. (Has sense if the singular sites were given directly)
%> @return [1] A class Coordinates containing the coordinates of the sites of the effective Hamiltonian
    function coordinates = Get_Effective_Coordinates(obj)
        if isempty(obj.next_SVD_cycle)
            coordinates = obj.coordinates;
        elseif strcmpi( class(obj.next_SVD_cycle), 'SVDregularizationLead' )
            coordinates = obj.next_SVD_cycle.Get_Effective_Coordinates();
        else
            error('EQuUs:SVDregularizationLead:Get_Effective_Coordinates', 'Unrecognized type of attribute next_SVD_cycle');
        end
        
    end
    
    
%% SVD_transform
%> @brief Regularize the Hamiltonians of the lead by SVD regularization
    function SVD_transform( obj  ) 
        
        % determine the Hamiltonians K0 and K1, K1adj used in transport calculations
        [K0, K1, K1adj] = obj.qDependentHamiltonians();        
        
        obj.SVDdecompozition( K1 );
        
        % Eqs (48) in PRB 78, 035407 
        % here we use diiferent transformation for the leads with diffferent orientation
        if obj.Lead_Orientation == 1
            U = obj.V;
        elseif obj.Lead_Orientation == -1
            U = obj.U;
        end
        K0 = U'*K0*U;
        K1 = U'*K1*U;
        K1adj = U'*K1adj*U;
        
        obj.Neff = size(obj.K0, 1);    
        
        % determine singular values
        S = abs(obj.S);
        regular_sites_slab = transpose(find(S > obj.tolerance*max(S)));  
        
        K = [K0, K1; K1adj, K0];       
        
        
        regular_sites = [regular_sites_slab, size(K0,1)+regular_sites_slab];
                     
        CreateH = CreateHamiltonians(obj.Opt, obj.param);
        CreateH.Write('Hscatter', K);
        CreateH.Write('kulso_szabfokok', regular_sites);
        CreateH.Write('HamiltoniansCreated', true);
        CreateH.Write('HamiltoniansDecimated', false);

        Decimation_Procedure = Decimation(obj.Opt); 
        Decimation_Procedure.DecimationFunc(0, CreateH, 'Hscatter', 'kulso_szabfokok');
        
        K = CreateH.Read('Hscatter');
        CreateH.Clear('Hscatter'); 
        
        Neff_new = length(regular_sites_slab);     
        
        if obj.Lead_Orientation == 1
            K0 = K(Neff_new+1:end, Neff_new+1:end);
        elseif obj.Lead_Orientation == -1            
            K0 = K(1:Neff_new, 1:Neff_new);        
        end
        K1 = K(1:Neff_new, Neff_new+1:end);
        K1adj = K(Neff_new+1:end, 1:Neff_new);     
        
        
        obj.next_SVD_cycle = SVDregularizationLead(obj.Opt, obj.param, obj.varargin{:});
        obj.next_SVD_cycle.Write('K0', K0);
        obj.next_SVD_cycle.Write('K1', K1);
        obj.next_SVD_cycle.Write('K1adj', K1adj);
        obj.next_SVD_cycle.Write('OverlapApplied', true); 
        obj.next_SVD_cycle.Write('HamiltoniansDecimated', true);     
        obj.next_SVD_cycle.Write('Neff', Neff_new); 
        
        
        obj.is_SVD_transformed = true;
        obj.HamiltoniansDecimated = true;
        
        if issparse(K1)
            condition_num = 1/condest(K1);
        else
            condition_num = rcond(K1);
        end
        
        if condition_num<obj.tolerance
            obj.next_SVD_cycle.Calc_Effective_Hamiltonians( 0 );
        end             
    end
    
%% Decimate_Hamiltonians
%> @brief Decimates the Hamiltonians (if the singular sites are predefined).
    function Decimate_Hamiltonians( obj ) 
            % determine the Hamiltonians K0 and K1, K1adj used in transport calculations
            [K0loc, K1loc, K1adjloc] = obj.qDependentHamiltonians();   
            
            K = [K0loc, K1loc;
                    K1adjloc, K0loc];       
     
            regular_sites_slab = obj.kulso_szabfokok;
            regular_sites = [regular_sites_slab, size(obj.K0,1)+regular_sites_slab];            

            
            CreateH = CreateHamiltonians(obj.Opt, obj.param);
            CreateH.Write('Hscatter', K);
            CreateH.Write('kulso_szabfokok', regular_sites);
            CreateH.Write('HamiltoniansCreated', true);
            CreateH.Write('HamiltoniansDecimated', false);

            Decimation_Procedure = Decimation(obj.Opt); 
            Decimation_Procedure.DecimationFunc(0, CreateH, 'Hscatter', 'kulso_szabfokok');
        
            K = CreateH.Read('Hscatter');
            CreateH.Clear('Hscatter'); 
            coordinates = CreateH.Read('coordinates');
        
            Neff_new = length(regular_sites_slab);     
        
            K0 = K(Neff_new+1:end, Neff_new+1:end);
            K1 = K(1:Neff_new, Neff_new+1:end);
            K1adj = K(Neff_new+1:end, 1:Neff_new); 
                        
            obj.next_SVD_cycle = SVDregularizationLead(obj.Opt, obj.param, obj.varargin{:});
            obj.next_SVD_cycle.Write('K0', K0);
            obj.next_SVD_cycle.Write('K1', K1);
            obj.next_SVD_cycle.Write('K1adj', K1adj);
            obj.next_SVD_cycle.Write('OverlapApplied', true); 
            obj.next_SVD_cycle.Write('HamiltoniansDecimated', true);     
            obj.next_SVD_cycle.Write('Neff', Neff_new); 
            
         
            obj.HamiltoniansDecimated = true;
            
            
            % determine the truncated coordinates
            indexes2remove = true( size( obj.coordinates.x ) );
            indexes2remove( regular_sites_slab ) = false;
            coordinates = obj.coordinates.RemoveSites( indexes2remove );
            obj.next_SVD_cycle.Write('coordinates', coordinates); 
            
            
            
%             Decimation_Procedure_Lead = Decimation(obj.Opt, 'ForcedDecimation', obj.Opt.Decimation);
%             kulso_szabfokok_tmp = obj.kulso_szabfokok;
%             obj.kulso_szabfokok = [obj.kulso_szabfokok, obj.kulso_szabfokok+size(obj.H0,1)];       
%             Decimation_Procedure_Lead.DecimationFunc(0, obj, 'Hamiltonian2Dec', 'kulso_szabfokok'); 
%             obj.kulso_szabfokok = kulso_szabfokok_tmp;
%             obj.Neff = size(obj.K0,1);        
        
    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)        
        
        if isempty(obj.next_SVD_cycle)
            if ~isempty(obj.K0)
                obj.K0 = Umtx*obj.K0*Umtx';
            end
            
            if ~isempty(obj.K1)
                obj.K1 = Umtx*obj.K1*Umtx';
            end
            
            if ~isempty(obj.K1adj)
                obj.K1adj = Umtx*obj.K1adj*Umtx';
            end            
            
        elseif strcmpi( class(obj.next_SVD_cycle), 'SVDregularizationLead' )
            obj.next_SVD_cycle.Unitary_Transform(Umtx);
        else
            error('EQuUs:SVDregularizationLead:Unitary_Transform', 'Unrecognized type of attribute next_SVD_cycle');
        end
        
    end    
    
%% Get_Neff
%> @brief Gets the effective number of sites after the elimination of the singular values.
%> @return Returns with the effective number of sites
    function Neff = Get_Neff(obj)
        if isempty(obj.next_SVD_cycle)
            Neff = obj.Neff;
        elseif strcmpi( class(obj.next_SVD_cycle), 'SVDregularizationLead' )
            Neff = obj.next_SVD_cycle.Get_Neff();
        else
            error('EQuUs:SVDregularizationLead:Get_Neff', 'Unrecognized type of attribute next_SVD_cycle');
        end    
    end



    
%% Get_U
%> @brief Gets the total transformation U related to the SVD transformation
%> @return Returns with the total transformation U
function Vtot = Get_V(obj)
        if isempty(obj.next_SVD_cycle)
            Vtot = obj.V;
        elseif strcmpi( class(obj.next_SVD_cycle), 'SVDregularizationLead' )
            Vtot_tmp = obj.next_SVD_cycle.Get_V();
            size_V = size(obj.V);
            size_Vtmp = size(Vtot_tmp);
            Vtot_tmp = [Vtot_tmp, zeros(size_Vtmp(1), size_V(2)-size_Vtmp(2));
                        zeros(size_V(1)-size_Vtmp(1), size_Vtmp(2)), eye(size_V(1)-size_Vtmp(1))];
            Vtot = obj.V*Vtot_tmp;
        else
            error('EQuUs:SVDregularizationLead:Get_V', 'Unrecognized type of attribute next_SVD_cycle');
        end
        
end     
    
%% CreateClone
%> @brief Creates a clone of the present object.
%> @return Returns with the cloned object.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'empty' Set true to create an empty clone, or false (default) to clone all atributes.
%> @return Returns with the cloned object.
    function ret = CreateClone( obj, varargin )        
        
        p = inputParser;
        p.addParameter('empty', false );
        
        p.parse(varargin{:});       
        empty    = p.Results.empty;
        
        ret = SVDregularizationLead(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), 'SVDregularizationLead' )
            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 interface.
%> @param MemberName The name of the attribute to be set.
%> @param input The value to be set.
    function Write(obj, MemberName, input)       
        
        
        if strcmpi(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 interface.
%> @param MemberName The name of the attribute to be set.
%> @return Returns with the value of the attribute.
    function ret = Read(obj, MemberName)
        
        if strcmpi(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 interface.
%> @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=protected)
          
    
%% Extend_Wavefnc ---- DEVELOPMENT
%> @brief Extend a reduced wave function to the original basis before the SVD regularization  (Eq (45) in PRB 78 035407
%> @param wavefnc_reduced The reduced wavefunction of the effective system
%> @param expk e^(i*k)
%> @return A matrix conatining the extended wave functions.
    function wavefnc_extended = Extend_Wavefnc(obj, wavefnc_reduced, expk) 
        
        if ~isempty(obj.next_SVD_cycle)  && strcmpi( class(obj.next_SVD_cycle), 'SVDregularizationLead' )
            wavefnc_reduced = obj.next_SVD_cycle.Extend_Wavefnc(wavefnc_reduced, expk);
        elseif ~isempty(obj.next_SVD_cycle)
            error('EQuUs:SVDregularizationLead:Extend_Wavefnc', 'Unrecognized type of attribute next_SVD_cycle');
        end
        
        Fn = -obj.invD*(obj.K1adju/expk + obj.C);
        wavefnc_extended = [wavefnc_reduced; Fn*wavefnc_reduced]; 
        wavefnc_extended = obj.U*wavefnc_extended;
        
    end
    
    
    
%% Initialize
%> @brief Initializes class properties.
    function Initialize(obj)
        Initialize@CreateLeadHamiltonians(obj);
        obj.tolerance = 1e-12;
        obj.is_SVD_transformed = false;        
    end    
        



end % methods protected


    
    
    
end