Lead.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - Lead
%    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 Lead.m
%> @brief A class to calculate the Green functions and self energies of a translational invariant lead
%> @} 
%> @brief A class to calculate the Green functions and self energies of a translational invariant lead
%> 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 Lead < handle & CommonFunctions & EigenProblemLead %NEW
    
    
properties ( Access = protected )
     %>function handle to calculate the self energy
    SelfEnergyCalc
     %>Function handle to calculate the retarded infinite Greens function.
    GinfCalc    
     %>The normalization matrix in the Greens function.
    Normamtx
     %>The inverse of the normalization matrix.
    invNormamtx
     %>The matrix of the surface Greens function of the semi-infinite ribbon.
    gsurf
     %>The inverse of the retarded surface Greens function.
    gsurfinv
     %>The matrix of the retarded Greens function of the infinite ribbon.
    ginf
     %>The matrix of the retarded surface Greens function of a finite ribbon.
    gfin
     %>The inverse of the retarded surface Greens function of a finite ribbon.
    gfininv
     %>The retarded self-energy of the semi-infinite ribbon.
    Sigma
end





methods ( Access = public )
%% Lead
%> @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. For details see #EigenProblemLead.EigenProblemLead.
%> @return An instance of the class
    function obj = Lead(Opt, param, varargin)    
        obj = obj@EigenProblemLead( Opt, param, varargin{:} ); 
    end

%% InfGreenFunction
%> @brief Calculates the Green function of the infinite lead for eregy #E between slabs z1 and z2. The calculated Greens function is stored by attribute #ginf. 
%> @param z1 The index of the first slab.
%> @param z2 The index of the second slab.
%> @param varargin Optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'z1points' Array of sites in the slab z1 to include in the calculations. By default each site in the given slab is included.
%> @param 'z2points' Array of sites in the slab z2 to include in the calculations. By default each site in the given slab is included.
%> @return Return with the  Green function of the infinite lead between the given slabs.
    function ret = InfGreenFunction(obj, z1,z2, varargin)
        obj.szetvalaszto();
        
        p = inputParser;
        if z1 >= z2
            p.addParameter('z1points', 1:size(obj.modusmtx_p,1));
            p.addParameter('z2points', 1:size(obj.modusmtx_p,2));
        else
            p.addParameter('z1points', 1:size(obj.modusmtx_m,1));
            p.addParameter('z2points', 1:size(obj.modusmtx_m,2));            
        end
        p.parse(varargin{:});

        z1points = p.Results.z1points;
        z2points = p.Results.z2points;
        
        obj.ginf = obj.GinfCalc(z1, z2, z1points, z2points);
        ret = obj.ginf;
        
%        d_modusmtx_m = [];
%        d_modusmtx_p = [];
    end
    
end

methods ( Access = protected )
%%    GinfCalcWithDualModes     
%> @brief Calculates the Green function of the infinite ribbon with dual modes. The calculated Greens function is stored by attribute #ginf. 
%> @param z1 The index of the first slab.
%> @param z2 The index of the second slab.
%> @param z1points Site indexes in the slab at z1.
%> @param z2points Site indexes in the slab at z2.
%> @return Return with the  Greens function of the infinite lead between the given slabs.
    function ginf = GinfCalcWithDualModes(obj, z1, z2, z1points, z2points) 
        obj.NormamtxSzamolo();
        obj.ginf = 0;
        
        if z1 >= z2
            ginf = obj.modusmtx_p(z1points,:)*diag(obj.expk_p.^(obj.Lead_Orientation*(z1-z2)))*obj.d_modusmtx_p;
        else
%             for idx = 1:M          
%                 xi_m    = modusmtx_m(z1points,idx);
%                 d_xi_m  = d_modusmtx_m(idx,:);
%                 ginf = ginf + xi_m*d_xi_m*exp(1i*k_m(idx)*Lead_Orientation*(z1-z2));  
%             end
            ginf = obj.modusmtx_m(z1points,:)*diag(obj.expk_m.^(obj.Lead_Orientation*(z1-z2)))*obj.d_modusmtx_m;
        end
        
        ginf = ginf*obj.invNormamtx(:,z2points);
    end
    %--------------------------------------------------------------------------
    
%%    GinfCalcWithLeftModes
%> @brief Calculates the Green function of the infinite ribbon with with the left-sided modes. The calculated Greens function is stored by attribute #ginf. 
%> @param z1 The index of the first slab.
%> @param z2 The index of the second slab.
%> @param z1points Site indexes in the slab at z1.
%> @param z2points Site indexes in the slab at z2.
%> @return Return with the  Greens function of the infinite lead between the given slabs.
    function ginf = GinfCalcWithLeftModes(obj, z1, z2, z1points, z2points) 
      
        if z1 >= z2
            ginf = obj.modusmtx_p(z1points,:)*diag(obj.expk_p.^(obj.Lead_Orientation*(z1-z2))./(1i*obj.vcsop_p*obj.Lead_Orientation))*obj.modusmtx_p_left(:,z2points);
        else
            ginf = obj.modusmtx_m(z1points,:)*diag(obj.expk_m.^(obj.Lead_Orientation*(z1-z2))./(-1i*obj.vcsop_m*obj.Lead_Orientation))*obj.modusmtx_m_left(:,z2points);
        end    
        
    end
    %--------------------------------------------------------------------------
    
end    


methods ( Access = public )
    
%% Gamma
%> @brief Calculates the effective coupling of the lead to the scattering region according to Eq (3) in Eur. Phys. J. B 53, 537-549 (2006)
%> @return Return with the  matrix of the effective coupling.
    function ret = Gamma( obj )
       if isempty( obj.Sigma )
           obj.SelfEnergy();
       end
       
       ret = 1i*( obj.Sigma - (obj.Sigma)' );
       
       
    end
%% SelfEnergy
%> @brief Calculates the retarded self energy of the semi-infinite lead according to Eq (36) of Ref PRB 78, 035407 (2008). The calculated self energy is stored by attribure #Sigma. 
%> @return Returns with the calculated self energy.
    function ret = SelfEnergy( obj )
        obj.szetvalaszto();        
        
        obj.Sigma = obj.SelfEnergyCalc();
        ret = obj.Sigma;
        
    end
    
    
end


methods ( Access = protected )    

%% SelfEnergyCalcWithDualModes
%> @brief Calculates the self energy with dual modes according to Eqs (17) and (36) in PRB 78 035407
%> @return Return with the self energy.
    function Sigma = SelfEnergyCalcWithDualModes( obj ) 
    %>Eq (17) and (36) in PRB 78 035407
    %>g00 = inv(Normamtx)
        obj.NormamtxSzamolo();      
        [~, K1_eff, K1adj_eff] = obj.Get_Effective_Hamiltonians();
        Sigma = obj.modusmtx_m*diag(obj.expk_m.^(-obj.Lead_Orientation))*obj.d_modusmtx_m;
        if obj.Lead_Orientation == 1
            Sigma = K1adj_eff*Sigma;
        else
            Sigma =  K1_eff*Sigma;
        end
    end
    %--------------------------------------------------------------------------

%% SelfEnergyCalcWithLeftModes
%> @brief Calculates the self energy with the left-sided modes according to Eqs (20), (17) and (36) in PRB 78 035407
%> @return Return with the self energy.
    function Sigma = SelfEnergyCalcWithLeftModes( obj ) 
    %>Eq (20), (17) and (36) in PRB 78 035407
    %>g00 = inv(Normamtx)
        obj.NormamtxSzamolo(); 
        [~, K1_eff, K1adj_eff] = obj.Get_Effective_Hamiltonians();
        Sigma = obj.modusmtx_m*diag(obj.expk_m.^(-obj.Lead_Orientation)./(-1i*obj.vcsop_m*obj.Lead_Orientation))*obj.modusmtx_m_left;
        if obj.Lead_Orientation == 1
            Sigma = K1adj_eff*Sigma*obj.Normamtx;
        else
            Sigma = K1_eff*Sigma*obj.Normamtx;
        end                
    end
    %--------------------------------------------------------------------------  
    
end

methods ( Access = public )

%% FiniteGreenFunction
%> @brief Calculates the Green function of a finite piece of the lead for eregy #E between slabs z1 and z2. The calculated Greens function is stored by attribute #gfin.
%> @param z1 The index of the first slab.
%> @param z2 The index of the second slab.
%> @param varargin Optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'onlygfininv' If true, only the inverse of the Green function is calculated. (Default value is false)
%> @return Return with the  Green function of the infinite lead between the given slabs.
    function ret = FiniteGreenFunction(obj, z1,z2, varargin)
        
        p = inputParser;
        p.addParameter('onlygfininv', false );
        
        p.parse(varargin{:});       
        onlygfininv    = p.Results.onlygfininv;
        
        obj.szetvalaszto();
        obj.NormamtxSzamolo();   
        
        Neff = obj.Get_Neff();        
        obj.gfin = zeros(4*Neff,4*Neff);
        
        if z1 >= z2
            points = 1:size(obj.modusmtx_p,1);
            dpoints = 1:size(obj.modusmtx_p,1);
        else
            points = 1:size(obj.modusmtx_m,1);
            dpoints = 1:size(obj.modusmtx_m,1);  
        end        
         
        gtmp = obj.GinfCalc(z1-1, z1-1, points, dpoints);
        obj.gfin(1:Neff,1:Neff)              = gtmp;
        obj.gfin(1+Neff:2*Neff, 1+Neff:2*Neff)     = gtmp;
        obj.gfin(1+2*Neff:3*Neff, 1+2*Neff:3*Neff) = gtmp;              
        obj.gfin(1+3*Neff:4*Neff, 1+3*Neff:4*Neff) = gtmp;
        
        gtmp = obj.GinfCalc(z1-1, z1, points, dpoints);
        obj.gfin(1:Neff,1+Neff:2*Neff) = gtmp;
        obj.gfin(1+2*Neff:3*Neff,1+3*Neff:4*Neff) = gtmp;
        
        gtmp = obj.GinfCalc(z1, z1-1, points, dpoints);
        obj.gfin(1+Neff:2*Neff,1:Neff) = gtmp;
        obj.gfin(1+3*Neff:4*Neff, 1+2*Neff:3*Neff) = gtmp;    
        
        gtmp = obj.GinfCalc(z1-1, z2, points, dpoints);
        obj.gfin(1:Neff,1+2*Neff:3*Neff) = gtmp; 
        
        gtmp = obj.GinfCalc(z2, z1-1, points, dpoints);
        obj.gfin(1+2*Neff:3*Neff,1:Neff) = gtmp; 
                
        gtmp = obj.GinfCalc(z1-1, z2+1, points, dpoints);
        obj.gfin(1:Neff,1+3*Neff:4*Neff) = gtmp; 
                
        gtmp = obj.GinfCalc(z2+1, z1-1, points, dpoints);
        obj.gfin(1+3*Neff:4*Neff,1:Neff) = gtmp;   
        
        gtmp = obj.GinfCalc(z1, z2, points, dpoints);
        obj.gfin(1+Neff:2*Neff,1+2*Neff:3*Neff) = gtmp; 
        
        gtmp = obj.GinfCalc(z2, z1, points, dpoints);
        obj.gfin(1+2*Neff:3*Neff,1+Neff:2*Neff) = gtmp; 
        
        gtmp = obj.GinfCalc(z1, z2+1, points, dpoints);
        obj.gfin(1+Neff:2*Neff,1+3*Neff:4*Neff) = gtmp; 
        
        gtmp = obj.GinfCalc(z2+1, z1, points, dpoints);
        obj.gfin(1+3*Neff:4*Neff,1+Neff:2*Neff) = gtmp;   
        
        rcond_gfin = rcond(obj.gfin);
        if isnan(rcond_gfin) || abs( rcond_gfin ) < 1e-15
            obj.display( 'EQuUs:Lead:FiniteGreenFunction: Regularizing gfin by SVD',1);
            obj.gfininv = obj.inv_SVD( obj.gfin );
        else
            obj.gfininv = inv(obj.gfin);
        end       
        obj.gfininv = obj.gfininv( 1+Neff:3*Neff, 1+Neff:3*Neff );
        
        if onlygfininv
            ret = [];
        else
            rcond_gfinin = rcond(obj.gfininv);
            if isnan(rcond_gfinin ) || abs( rcond_gfinin ) < 1e-15
                obj.display( 'EQuUs:Lead:FiniteGreenFunction: Regularizing gfininv by SVD',1);
                obj.gfin = obj.inv_SVD( obj.gfininv );
            else
                obj.gfin = inv(obj.gfininv);
            end       
            obj.gfininv = [];
            ret = obj.gfin;
        end
        
    end

%% SurfaceGreenFunction
%> @brief Calculates the surface Green function of a semi-infinite lead for eregy #E. The calculated Green function is stored in attribute #gsurf.
%> @param varargin Optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'OnlyInverse' If true (default), only the inverse of the Green function is calculated.
%> @param 'CalcInverse' If true (default), the inverse of the Green function is also calculated.
    function SurfaceGreenFunction( obj, varargin ) 

     p = inputParser;
     p.addParameter('OnlyInverse', true ); 
     p.addParameter('CalcInverse', true );
     p.parse(varargin{:});
     OnlyInverse = p.Results.OnlyInverse;     
          CalcInverse = p.Results.CalcInverse;     
            
        obj.szetvalaszto();
        obj.NormamtxSzamolo(); 
        if obj.Opt.usingDualModes
            obj.GsurfSzamoloWithDualModes();
        else
            obj.GsurfSzamoloWithLeftModes();
        end 

          if ~CalcInverse
               return
          end

        rcond_gsurf = rcond(obj.gsurf);
        if isnan(rcond_gsurf) || abs( rcond_gsurf ) < 1e-15
            obj.display( ['EQuUs:', class(obj), ':SurfaceGreenFunction: Regularizing gsurf by SVD'],1);
            obj.gsurfinv = obj.inv_SVD( obj.gsurf );
        else
            obj.gsurfinv = inv(obj.gsurf);
        end  

          if  OnlyInverse            
          obj.gsurf = [];
          end      
        
    end
    
end

methods ( Access = protected )
    
%> @brief Calculates the surface Green function with dual modes. The calculated Green function is stored by attribute #gsurf.
    function GsurfSzamoloWithDualModes( obj )  
            
            %test for the self energy
%              obj.SelfEnergy();
%              [K0_eff] = obj.Get_Effective_Hamiltonians();
%              gsurf2 = -inv(K0_eff + obj.Sigma);
%         
            
        gsurf_m = obj.modusmtx_m*diag(obj.expk_m.^(-obj.Lead_Orientation))*obj.d_modusmtx_m;
        gsurf_p = obj.modusmtx_p*diag(obj.expk_p.^(obj.Lead_Orientation))*obj.d_modusmtx_p;
             

        %gsurf = (eye(M) - gsurf_p*gsurf_m)*inv(Normamtx);
        obj.gsurf = (eye(obj.Get_Neff()) - gsurf_m*gsurf_p)*obj.invNormamtx;
        
    end
    %-------------------------------------------------------------------------- 
        
        
%> @brief Calculates the surface Green function with the left-sided modes according to Eq (33) in PRB 78, 035407 (2008). The calculated Green function is stored by attribute #gsurf.
    function GsurfSzamoloWithLeftModes( obj )  
          Neff = obj.Get_Neff();            
        points = 1:Neff;
        dpoints = 1:Neff;
            
    %>right SGF Eq (33) in PRB 78, 035407 (2008)
    %>g00 = inv(Normamtx)
%         g01 = GinfCalc(0, 1, points, dpoints);
%         g10 = GinfCalc(1, 0, points, dpoints);
%         gsurfR = invNormamtx - g10*Normamtx*g01;
            
    %>left SGF Eq (33) in PRB 78, 035407 (2008)
    %>g00 = inv(Normamtx)
        g0m1 = obj.GinfCalc(0, -1, points, dpoints);
        gm10 = obj.GinfCalc(-1, 0, points, dpoints);
        obj.gsurf = obj.invNormamtx - gm10*obj.Normamtx*g0m1;
                       
        
    end
    %--------------------------------------------------------------------------  
    
end

methods ( Access = public )
     
%%    NormamtxSzamolo
%> @brief The normalization matrix to calculate the Green function. The calculated normalization matrix is stored by attribute #Normamtx.
    function NormamtxSzamolo( obj )
            
            if ~isempty(obj.Normamtx) && ~isempty(obj.invNormamtx)
               return 
            elseif obj.Opt.usingDualModes
                obj.NormamtxSzamoloWithDualModes();
            else
                obj.NormamtxSzamoloWithLeftModes();
            end       
    end
    
end
                        

methods ( Access = protected )
%% NormamtxSzamoloWithDualModes
%> @brief Calculates the normalization matrix to calculate the Green function using the dual modes. The calculated normalization matrix is stored by attribute #Normamtx.
    function NormamtxSzamoloWithDualModes(obj)
            
        obj.calcDualModes();  
        obj.Normamtx = 0;
        obj.invNormamtx = 0;
            
            
        obj.Normamtx = (obj.modusmtx_p*diag(obj.expk_p.^(-1))*obj.d_modusmtx_p - ...
                    obj.modusmtx_m*diag(obj.expk_m.^(-1))*obj.d_modusmtx_m)*obj.Lead_Orientation;        
                
        [~, ~, K1adj_eff] = obj.Get_Effective_Hamiltonians();
        obj.Normamtx = K1adj_eff*obj.Normamtx;
        
        try
            rcond_Normamtx = rcond(obj.Normamtx);
            if isnan(rcond_Normamtx) || abs( rcond_Normamtx ) < 1e-15
                obj.display( 'EQuUs:Lead:NormamtxSzamoloWithDualModes: Regularizing Normamtx by SVD',1);
                obj.invNormamtx = obj.inv_SVD( obj.Normamtx );
            else
                obj.invNormamtx = obj.Normamtx\(eye(size(obj.Normamtx)));%inv(Normamtx);
            end                
        catch errCause
            err = MException('EQuUs:Lead:NormamtxSzamoloWithDualModes', 'The inversion of the Normamtx failed');
            err = addCause(err, errCause);
            save('Error_Lead_NormamtxSzamoloWithDualModes.mat');
            throw(err); 
        end
    end
    %---------------------------------------------------------------
%% NormamtxSzamoloWithLeftModes
%> @brief Calculates the normalization matrix to calculate the Green function using the left-sided modes. The calculated normalization matrix is stored by attribute #Normamtx.
    function NormamtxSzamoloWithLeftModes(obj)
           
        obj.invNormamtx = obj.InfGreenFunction(0,0);
        rcond_invNormamtx = rcond(obj.invNormamtx);
        if isnan(rcond_invNormamtx) || abs( rcond_invNormamtx ) < 1e-15
            obj.display( 'EQuUs:Lead:NormamtxSzamoloWithLeftModes: Regularizing invNormamtx by SVD',1);
            obj.Normamtx = obj.inv_SVD( obj.invNormamtx );
        else
            obj.Normamtx    = obj.invNormamtx\(eye(size(obj.invNormamtx)));
        end                
           
            
    end

%% calcDualModes
%> @brief Calculates the dual modes. The dual modes is stored by attributes #d_modusmtx_m and #d_modusmtx_p.
    function  calcDualModes( obj )
        if isempty(obj.d_modusmtx_m) 
            obj.display('EQuUs:Lead:calcDualModes: Calculating dual modes for m')
            if obj.Opt.usingDualModes
                    obj.d_modusmtx_m = obj.CalcDualModesEigenvec( obj.modusmtx_m ); %a dualis modusok
            else
                    obj.d_modusmtx_m = obj.CalcDualModesEigenvecWithLeftModes( obj.modusmtx_m_left, -obj.vcsop_m*obj.Lead_Orientation ); %a dualis modusok
            end
       end
       
       if isempty(obj.d_modusmtx_p) 
            obj.display('EQuUs:Lead:calcDualModes: Calculating dual modes for p')
               if obj.Opt.usingDualModes
                    obj.d_modusmtx_p = obj.CalcDualModesEigenvec( obj.modusmtx_p ); %a dualis modusok
               else
                    obj.d_modusmtx_p = obj.CalcDualModesEigenvecWithLeftModes( obj.modusmtx_p_left, obj.vcsop_p*obj.Lead_Orientation ); %a dualis modusok
               end
       end

%    obj.display('modusmtx inverse',1);     
%      obj.display(num2str(norm( obj.d_modusmtx_m*obj.modusmtx_m - eye(obj.Get_Neff()) )), 1);
%      obj.display(num2str(norm( obj.d_modusmtx_p*obj.modusmtx_p - eye(obj.Get_Neff()) )), 1);   
    end


%% CalcDualModesEigenvec
%> @brief Calculates the dual modes by inverting the right-sided eigenvectors. 
%> @param modusmtx The right-sided eigenvectors. 
%> @return Returns with the calculated dual modes.
    function d_modes = CalcDualModesEigenvec( obj, modusmtx )  
        try
            rcond_modusmtx = rcond(modusmtx);
            if isnan(rcond_modusmtx ) || abs( rcond_modusmtx ) < 1e-15
                obj.display( 'EQuUs:Lead:CalcDualModesEigenvec: Regularizing modusmtx by SVD',1);
                d_modes = obj.inv_SVD( modusmtx );
            else
                d_modes = (modusmtx)\eye(size(modusmtx)); % a dualis modusok
            end  
                        
        catch errCause
            err = MException('EQuUs::Lead::CalcDualModesEigenvec', 'The calculation of the dual wavefunctions failed.');
            err = addCause(err, errCause);
            save('Error_Lead_CalcDualModesEigenvec.mat');
            throw(err); 
        end
            
    end

%% CalcDualModesEigenvecWithLeftModes
%> @brief Calculates the dual modes using the left-sided eigenvectors. 
%> @param modusmtx_left The left-sided eigenvectors. 
%> @param vcsop The group velocities.
%> @return Returns with the calculated dual modes
    function d_modes = CalcDualModesEigenvecWithLeftModes( obj, modusmtx_left, vcsop ) 
    %>Eq (20) in  PRB 78, 035407 (2008)
          d_modes = diag(1./(1i*vcsop))*modusmtx_left*obj.Normamtx;
    end    
    
    
%% Initialize
%> @brief Initializes class attributes.
    function Initialize(obj)
        
        Initialize@EigenProblemLead(obj);
        
        if obj.Opt.usingDualModes
            obj.display( 'EQuUs:Lead:Initialize: Using dual modes in the calculations.');
            obj.GinfCalc = @(z1, z2, z1points, z2points)(obj.GinfCalcWithDualModes(z1, z2, z1points, z2points));
            obj.SelfEnergyCalc = @()(obj.SelfEnergyCalcWithDualModes);
        else
            obj.display( 'EQuUs:Lead:Initialize: Using left-sided modes in the calculations.');
            obj.GinfCalc = @(z1, z2, z1points, z2points)(obj.GinfCalcWithLeftModes(z1, z2, z1points, z2points));
            obj.SelfEnergyCalc = @()(obj.SelfEnergyCalcWithLeftModes);        
        end  
    end    
    
end
   
methods ( Access = public )

%% TrukkosSajatertekek
%> @brief Calculates the wave numbers corresponding to the propagating states at given energy. The calculated wave numbers are stored in attribute #expk. The normalization matrix is set to empty.
%> @param E The energy value used in the calculations.
    function TrukkosSajatertekek(obj, E)
        TrukkosSajatertekek@EigenProblemLead(obj, E);
        obj.Normamtx       = [];
        obj.invNormamtx    = [];        
    end
    
%% ShiftLead
%> @brief Shifts the on-site energies in the leads by a given energy. The normalization matrix is set to empty.
%> @param Energy The enrgy value.
    function ShiftLead( obj, Energy )        
        ShiftLead@EigenProblemLead( obj, Energy );
        obj.Normamtx = [];
    end

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

        AddPotential@EigenProblemLead( obj, V );
        obj.Normamtx = [];
    end
    
%% Unitary_Transform
%> @brief Transforms the effective Hamiltonians and the calculated surface Green operator and selfenergy by a unitary transformation
%> @param Umtx The matrix of the unitary transformation.
    function Unitary_Transform(obj, Umtx)
        Unitary_Transform@EigenProblemLead( obj, Umtx );
        
        if ~isempty(obj.gsurfinv)
            obj.gsurfinv = Umtx*obj.gsurfinv*Umtx';
        end
            
        if ~isempty(obj.Sigma)
            obj.Sigma = Umtx*obj.Sigma*Umtx';
        end  
        
        if ~isempty(obj.Normamtx)
            obj.Normamtx = Umtx*obj.Normamtx*Umtx';
        end  
        
    end        

%% CreateClone
%> @brief Creates a clone of the present Lead 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.
    function ret = CreateClone( obj, varargin )     
        
        p = inputParser;
        p.addParameter('empty', false );
        
        p.parse(varargin{:});       
        empty    = p.Results.empty;
        
        ret = Lead(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, 'SelfEnergyCalc' ) || strcmpi( prop_name, 'GinfCalc' )  
              continue
            elseif 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), 'Lead' )
            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  strcmpi(MemberName, 'd_modusmtx_p')
            obj.calcDualModes()
            ret = obj.d_modusmtx_p;
            return
        elseif strcmpi(MemberName, 'd_modusmtx_m')
            obj.calcDualModes()
            ret = obj.d_modusmtx_m;  
            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

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