InterfaceRegion.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - InterfaceRegion
%    Copyright (C) 2009-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 InterfaceRegion.m
%> @brief A class describing the interface region between the scattering region and a lead.
%> @} 
%> @brief A class describing the interface region between the scattering region and a lead.
%> @Available
%> EQuUs v4.8 or later
%%
classdef InterfaceRegion < handle & CommonFunctions & SVDregularizationLead
    
    
properties ( Access = protected )
    %> Effective Hamiltonian of the interface region
    K00 
    %> Hinterface -E*Sinterface
    Kinterface   
     %> Hamiltonian of the interface region
    Hinterface      
     %> Hcoupling-E*Scoupling
    Kcoupling   
     %> Hcouplingadj-E*Scoupling'
    Kcouplingadj    
     %> Coupling Hamiltonian from the interface region to the scattering region
    Hcoupling  
     %> Coupling Hamiltonian from the scattering region to the interface region
    Hcouplingadj        
     %> Overlap integrals corresponding to Hcoupling
    Scoupling  
     %> Overlap integrals corresponding to Hinterface
    Sinterface    
    %> Coordinates of the surface sites of the scattering region to be connected to. 
    coordinates2
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 = InterfaceRegion(Opt, param, varargin)    
        obj = obj@SVDregularizationLead( Opt, param, varargin{:} ); 

        if strcmpi(class(obj), 'InterfaceRegion')       
            obj.Reset(); 
        end       
    end
    
    
%% ApplyOverlapMatrices
%> @brief Applies the overlap matrices to the Hamiltonians: K = H-ES
%> @param E The energy value.
    function ApplyOverlapMatrices(obj, E)
        
        if obj.OverlapApplied
            obj.display('EQuUs:InterfaceRegion:ApplyOverlapMatrices: Overlap matrices were already applied.');
            return;
        end
        
        ApplyOverlapMatrices@CreateLeadHamiltonians( obj, E );
        
        if ~isempty( obj.Scoupling )
            obj.Kcoupling = obj.Hcoupling - E*obj.Scoupling; 
            obj.Kcouplingadj = obj.Hcouplingadj - E*obj.Scoupling'; 
        else
            obj.Kcoupling = obj.Hcoupling;
            obj.Kcouplingadj = obj.Hcouplingadj;
        end
        
        if ~isempty( obj.Sinterface )
            obj.Kinterface = obj.Hinterface - E*obj.Sinterface; 
        else
            obj.Kinterface = obj.Hinterface;
        end        
            
    end
    
    
%% Transform2BdG
%> @brief Transforms the Hamiltonians and the overlap matrices into the BdG model.
    function Transform2BdG( obj )
    
        if isempty(obj.Opt.BdG) || ~obj.Opt.BdG
            return
        end
        
        if ~isempty( obj.coordinates.BdG_u )
            % already transformed into BdG model
            return
        end
        
        Transform2BdG@CreateLeadHamiltonians( obj );
        
        if isempty( obj.coordinates2 )        
            fnames = fieldnames( obj.coordinates2 );
            for idx = 1:length(fnames)
                fname = fnames{idx};
                if strcmp(fname, 'a') || strcmp(fname, 'b')
                    continue
                end                    
                obj.coordinates2.(fname) = [ obj.coordinates2.(fname); obj.coordinates2.(fname) ];
            end
            obj.coordinates2.BdG_u = [ true(size(obj.H0,1),1); false(size(obj.H0,1),1)];
        end
                
                  
        pair_potential = obj.params.pair_potential;
        
        % transforming the coupling Hamiltonians and overlaps        
        if isempty(obj.Scoupling)
            Scoupling = sparse([],[],[], size(obj.Hcoupling,1), size(obj.Hcoupling,2));
        else
            Scoupling = obj.Scoupling;
        end                
        
        obj.Hcoupling = [obj.Hcoupling, Scoupling*pair_potential; Scoupling*conj(pair_potential), -conj(obj.Hcoupling)];
        obj.Hcouplingadj = [obj.Hcouplingadj, Scoupling'*pair_potential; Scoupling'*conj(pair_potential), -conj(obj.Hcouplingadj)];
        
        if ~isempty( obj.Scoupling )
            obj.Scoupling = [obj.Scoupling, sparse([], [], [], size(obj.Scoupling,1), size(obj.Scoupling,2)); sparse([], [], [], size(obj.Scoupling,1), size(obj.Scoupling,2)), obj.Scoupling];
        end
             
        
    end
    
%% Calc_Effective_Hamiltonians
%> @brief Calculates the effective Hamiltonians according to Eq (48) of of PRB 78, 035407 
%> @param E The energy value
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'Lead' An instance of class #SVDregularizationLead (or its subclass) providing a source for the SVD matrices.
    function Calc_Effective_Hamiltonians( obj, E, varargin )  
        
        p = inputParser;
        p.addParameter('Lead', []);
        p.parse(varargin{:});
        Lead  = p.Results.Lead;
                
        obj.ApplyOverlapMatrices(E);
        
        
        if ~isempty(obj.kulso_szabfokok) && length(obj.kulso_szabfokok) ~= size(obj.K0,1)
            obj.Decimate_Hamiltonians();
            return
        elseif Lead.Read('is_SVD_transformed')
               obj.SVD_transform( Lead );
               return
        else
            % in case no SVD regularization or simple decimation is needed 
            obj.Neff = size(obj.K0,1);
            obj.K00 = obj.K0;
            return
        end
        
       
        
                
    end    


%% Get_Effective_Hamiltonians
%> @brief Gets the effective Hamiltonians of the interface region
%> @return [1] The effective Hamiltonian of the unit cell 
%> @return [2] The effective coupling between the unit cells 
%> @return [3] The adjungate of the effective coupling between the unit cells 
%> @return [4] The effective coupling between the lead and the scattering region 
%> @return [5] The adjungate of the effective coupling between the lead and the scattering region 
    function [K0_eff, K1_eff, K1adj_eff, Kcoupling, Kcouplingadj] = Get_Effective_Hamiltonians(obj)
        if isempty(obj.next_SVD_cycle)
            K0_eff = obj.K00;
            if obj.Lead_Orientation == 1
                 K1_eff = obj.K1;
             K1adj_eff = obj.K1adj;
               elseif obj.Lead_Orientation == -1
                 K1_eff = obj.K1adj;
             K1adj_eff = obj.K1;
               end
               Kcoupling = obj.Kcoupling;
               Kcouplingadj = obj.Kcouplingadj;
        elseif strcmpi( class(obj.next_SVD_cycle), 'InterfaceRegion' )
            [K0_eff, K1_eff, K1adj_eff, Kcoupling, Kcouplingadj] = obj.next_SVD_cycle.Get_Effective_Hamiltonians();
        else
            error('EQuUs:InterfaceRegion:Get_Effective_Hamiltonians', 'Unrecognized type of attribute next_SVD_cycle');
        end
        
    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)
            obj.K00 = Umtx*obj.K00*Umtx';
            obj.K1 = Umtx*obj.K1*Umtx';
            obj.K1adj = Umtx*obj.K1adj*Umtx';
               %obj.Kcoupling = Umtx*obj.Kcoupling;
               %obj.Kcouplingadj = obj.Kcouplingadj*Umtx';  
            %5
        elseif strcmpi( class(obj.next_SVD_cycle), 'InterfaceRegion' )
            obj.next_SVD_cycle.Unitary_Transform(Umtx);
        else
            error('EQuUs:InterfaceRegion:Unitary_Transform', 'Unrecognized type of attribute next_SVD_cycle');
        end
        
    end        

%% SVD_transform
%> @brief Transforms the Hamiltonians of the interface region according to the SVD regularization in the lead.
%> @param Lead an instance of class #SVDregularizationLead or its subclass representing the attached lead.
    function SVD_transform( obj, Lead ) 

        
          obj.Neff = Lead.Read('Neff');

        if ~Lead.Read('is_SVD_transformed')
            return
        end
        
        if obj.Lead_Orientation == 1
            obj.V = Lead.Read('V'); 
            U = obj.V;            
        elseif obj.Lead_Orientation == -1   
            obj.U = Lead.Read('U'); 
            U = obj.U;             
        end        
               

        %> Performing SVD regularization     
        
        if length(obj.q) < 2 && ~isempty(obj.q) && ~obj.HamiltoniansDecimated
            K0 = obj.K0;
            K0 = K0 + obj.K1_transverse*diag(exp(1i*obj.q)) + obj.K1_transverse'*diag(exp(-1i*obj.q));
        else
            K0 = Lead.Read('K0');
        end
        
        % Eqs (48) in PRB 78, 035407 
        % here we use diiferent transformation for the leads with diffferent orientation   
        Neff = obj.Neff;
        K0 = U'*K0*U;
        K1 = obj.K1;
        K1adj = obj.K1adj;
        K1(1:Neff,1:Neff) = U'*obj.K1(1:Neff,1:Neff)*U;
        K1adj(1:Neff,1:Neff) = U'*obj.K1adj(1:Neff,1:Neff)*U;
        
        if isempty(obj.K00)
            K00 = K0;
        else
            K00 = obj.K00;
            K00(1:Neff, 1:Neff) =  U'*obj.K00(1:Neff, 1:Neff)*U;
            K00(1:Neff, end-Neff+1:end) = U'*obj.K00(1:Neff, end-Neff+1:end);
            K00(end-Neff+1:end, 1:Neff) = obj.K00(end-Neff+1:end, 1:Neff)*U;
        end
        
        S = abs(Lead.Read('S'));
        non_singular_sites_slab = transpose(find(S > obj.tolerance*max(S)));        
        
        if obj.Lead_Orientation == 1
            K = [K0, K1; K1adj, K00];
            non_singular_sites = [non_singular_sites_slab, size(K0,1)+1:size(K0,1)+size(K00,1)];
        elseif obj.Lead_Orientation == -1
            K = [K00, K1; K1adj, K0];
            non_singular_sites = [1:size(K00,1), size(K00,1) + non_singular_sites_slab];
        end        
          

          
                     
        CreateH = CreateHamiltonians(obj.Opt, obj.param);
        CreateH.Write('Hscatter', K);
        CreateH.Write('kulso_szabfokok', non_singular_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(non_singular_sites_slab);

        if obj.Lead_Orientation == 1            
            K00 = 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);            
        elseif obj.Lead_Orientation == -1
            K00 = K(1:end-Neff_new, 1:end-Neff_new);
            K1 = K(1:end-Neff_new, end-Neff_new+1:end);
            K1adj = K(end-Neff_new+1:end, 1:end-Neff_new);                     
        end
        
            
        obj.HamiltoniansDecimated = true;
        
        %> SVD transform the coupling Hamiltonians
        Kcoupling = obj.Kcoupling;
        Kcouplingadj = obj.Kcouplingadj;
        if obj.Lead_Orientation == 1
            Kcoupling(1:Neff,:) = U'*Kcoupling(1:Neff,:);
            Kcouplingadj(:,1:Neff) = Kcouplingadj(:,1:Neff)*U;
        else
            Kcoupling(1:Neff,:) = U'*Kcoupling(1:Neff,:);
            Kcouplingadj(:,1:Neff) = Kcouplingadj(:,1:Neff)*U;        
        end        
        
        
        obj.next_SVD_cycle = InterfaceRegion(obj.Opt, obj.param, obj.varargin{:});
        obj.next_SVD_cycle.Write('K00', K00);
        obj.next_SVD_cycle.Write('K1', K1);
        obj.next_SVD_cycle.Write('K1adj', K1adj);  
        obj.next_SVD_cycle.Write('Kcoupling', Kcoupling);
        obj.next_SVD_cycle.Write('Kcouplingadj', Kcouplingadj);   
        obj.next_SVD_cycle.Write('HamiltoniansDecimated', true); 
        obj.next_SVD_cycle.Write('Neff', Neff_new); 
        
        obj.is_SVD_transformed = true;
        obj.HamiltoniansDecimated = true;
        
        if ~isempty(obj.next_SVD_cycle)  && strcmpi( class(obj.next_SVD_cycle), 'InterfaceRegion' )
            obj.next_SVD_cycle.SVD_transform( Lead.Read('next_SVD_cycle') );
        elseif ~isempty(obj.next_SVD_cycle)
            error('EQuUs:InterfaceRegion:Calc_Effective_Hamiltonians', 'Unrecognized type of attribute next_SVD_cycle');
        end

     end   

%% Decimate_Interface
%> @brief Decimates the Hamiltonians of the interface region.
    function Decimate_Hamiltonians( obj )  
        if ~isempty(obj.q)
            K0loc = obj.K0 + obj.K1_transverse*diag(exp(1i*obj.q)) + obj.K1_transverse'*diag(exp(-1i*obj.q));
        else
            K0loc = obj.K0;
        end            
        K = [K0loc,obj.K1;obj.K1adj,K0loc];    
        
        non_singular_sites_slab = obj.kulso_szabfokok; 
        Neff_new = length(non_singular_sites_slab);     
        
        
        non_singular_sites = [non_singular_sites_slab, size(K,1)/2+1:size(K,1)];
                     
        CreateH = CreateHamiltonians(obj.Opt, obj.param);
        CreateH.Write('Hscatter', K);
        CreateH.Write('kulso_szabfokok', non_singular_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'); 
        
        cCoordinates = obj.coordinates;
        coord_fieldnames = fieldnames(cCoordinates);
        if obj.Lead_Orientation == 1          
            K00 = 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);   
            
            % reordering the sites to get the non-singular sites at the top corner
            indexes = false( size(K00,1), 1);
            indexes(non_singular_sites_slab) = true;
            K00 = [K00(indexes, indexes), K00(indexes, ~indexes); K00(~indexes, indexes), K00(~indexes, ~indexes)];
            K1 = [K1(:, indexes), K1(:, ~indexes)];
            K1adj = [K1adj(indexes, :); K1adj(~indexes, :)];
            Kcoupling = obj.Kcoupling;
            Kcoupling = [Kcoupling(indexes, :); Kcoupling(~indexes, :)];
            Kcouplingadj = obj.Kcouplingadj;
            Kcouplingadj = [Kcouplingadj(:,indexes), Kcouplingadj(:, ~indexes)];
            
            % reordering the sites to get the non-singular sites at the top corner
            for idx = 1:length(coord_fieldnames)
                fieldname = coord_fieldnames{idx};
                if strcmpi(fieldname, 'a') || strcmpi(fieldname, 'b') || strcmpi(fieldname, 'LatticeConstant') 
                    continue
                elseif ~isempty(cCoordinates.(fieldname))
                    tmp = cCoordinates.(fieldname);
                    cCoordinates.(fieldname) = [tmp(indexes); tmp(~indexes)];
                end
            end            
            
        elseif obj.Lead_Orientation == -1  
            % keeping only regular sites
            K00 = K(1:Neff_new, 1:Neff_new);
            K1 = K(1:Neff_new, Neff_new+non_singular_sites_slab);
            K1adj = K(Neff_new+non_singular_sites_slab, 1:Neff_new);
            Kcoupling = obj.Kcoupling;
            Kcouplingadj = obj.Kcouplingadj;
            
            % keeping only regular sites
            cCoordinates = cCoordinates.KeepSites(non_singular_sites_slab);
        else
            error('EQuUs:InterfaceRegion:Decimate_Interface', 'Unknown lead orientation');            
        end    
        
        obj.next_SVD_cycle = InterfaceRegion(obj.Opt, obj.param, obj.varargin{:});
        obj.next_SVD_cycle.Write('K00', K00);
        obj.next_SVD_cycle.Write('K1', K1);
        obj.next_SVD_cycle.Write('K1adj', K1adj);  
        obj.next_SVD_cycle.Write('Kcoupling', Kcoupling);
        obj.next_SVD_cycle.Write('Kcouplingadj', Kcouplingadj);   
        obj.next_SVD_cycle.Write('HamiltoniansDecimated', true); 
        obj.next_SVD_cycle.Write('Neff', Neff_new); 
        obj.next_SVD_cycle.Write('coordinates', cCoordinates); 
          
            
        obj.HamiltoniansDecimated = true;
        obj.Neff = size(obj.K0,1);

     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), 'InterfaceRegion' )
            Neff = obj.next_SVD_cycle.Get_Neff();
        else
            error('EQuUs:InterfaceRegion:Get_Neff', 'Unrecognized type of attribute next_SVD_cycle');
        end
        
    end
   

%% CreateClone
%> @brief Creates a clone of the present class.
%> @return Returns with the cloned object.
    function ret = CreateClone( obj )
        
        
        ret = InterfaceRegion(obj.Opt, obj.param,...
                                            'Hanyadik_Lead', obj.Hanyadik_Lead,...
                                            'Lead_Orientation', obj.Lead_Orientation, ...
                                            'q', obj.q);  
                                        
        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), 'InterfaceRegion' )
            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              
        
        Reset@SVDregularizationLead( obj );   
             
    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 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 interface.
%> @param MemberName The name of the attribute to be set.
%> @return Returns with the value of the attribute.
    function ret = Read(obj, MemberName)
        
        try
            ret = obj.(MemberName);
          catch
               error(['EQuUs:', class(obj), ':Read'], ['No property to read with name: ',  MemberName]);
          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

%------------------------------------------------------------------
end % Global End