NTerminal.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - NTerminal
%    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 utilities Utilities
%> @{
%> @file NTerminal.m
%> @brief   A class describing an N-terminal geometry for equilibrium calculations mostly in the zero temperature limit.
%> @image html two_terminal_structure.jpg
%> @image latex two_terminal_structure.jpg
%> @} 
%> @brief  A class describing an N-terminal geometry for equilibrium calculations mostly in the zero temperature limit.
%> @Available
%> EQuUs v4.9 or later
%> <tr class="heading"><td colspan="2"><h2 class="groupheader"><a name="avail"></a>Structure described by the class</h2></td></tr>
%> @image html two_terminal_structure.jpg
%> @image latex two_terminal_structure.jpg
%> The drawing represents a two-terminal structure made of two leads, of a scattering region and of two interface regions between the leads and scattering center.
%> However, additional lead can be added to the structure.
%> Each rectangle describes a unit cell including singular and non-singular sites. 
%> The scattering center is, on the other hand, represented by a larger rectangle corresponding to the Hamiltonian of the scattering region.
%> Arrows indicate the hopping direction stored by the attributes in the corresponding classes (see attributes #CreateLeadHamiltonians.H1 and #InterfaceRegion.Hcoupling for details).
%> The orientation of the lead is +1 if the lead is terminated by the interface region in the positive direction, and -1 if the lead is terminated by the interface region in the negative direction.
%> (see attribute #CreateLeadHamiltonians.Lead_Orientation for details)
%%
classdef NTerminal < handle & Messages & CommonFunctions

    
    properties (Access = protected )   
        %> An instance of strucutre #param
        param
        %> Green operator of the scattering region
        G 
        %> The inverse of the Green operator G
        Ginv
        %> The energy used in the calculations
        E
        %> The Fermi energy. Attribute #E is measured from this value. (Use for equilibrium calculations in the zero temperature limit.) 
        EF
        %> The transverse momentum quantum number
        q    
        %> Function handle to create custom Hamiltonians.  Has the same inputs ans outputs as #Custom_Hamiltonians.LoadHamiltonians
        CustomHamiltoniansHandle        
    end
    
    
    properties (Access = public)
        %> An instance of class #CreateHamiltonians or its subclass to manipulate the Hamiltonian of the scattering region
        CreateH
        %> An instance of class #Transport_Interface (or its subclass) for transport calculations.
        FL_handles
        %> A cell array of classes #InterfaceRegion to describe the interface region between the leads and the scattering region
        Interface_Regions
        %> An instance of class #Peierls object to describe the peirls substitution in the scattering region
        PeierlsTransform_Scatter
        %> An instance of class #Peierls object to describe the peirls substitution in the leads
        PeierlsTransform_Leads
        %> Function handle S = f( x,y) of the gauge transformation in the scattering center. (S is a N x 1 vector, where N is the number of the points given by the x and y coordinates.)
        gauge_field
        %> function handle for individual physical model of the leads
        leadmodel
        %> function handle for individual physical model of the interface regions
        interfacemodel        
        %> Input filename containing the computational parameters. (Obsolete)
        filenameIn
        %> Output filename to export the computational parameters
        filenameOut
        %> A string of the working directory
        WorkingDir
        %> An instance of class #Custom_Hamiltonians to load the Hamiltonians from external source
        cCustom_Hamiltonians         
        %> if true, no output messages are print
        silent
    end


methods ( Access = public )
    
    
%% Contructor of the class
%> @brief Constructor of the class.
%> @param varargin Cell array of optional parameters. For details see #InputParsing.
%> @return An instance of the class
function obj = NTerminal( varargin )
    
    obj.G = []; % Greens function of the scattering region
    obj.Ginv = []; % inverse of G
    obj.param = [];
    obj.Opt  = [];
    
    obj.CreateH                  = [];
    obj.FL_handles               = [];
    obj.Interface_Regions        = [];
    obj.PeierlsTransform_Scatter = [];
    obj.PeierlsTransform_Leads   = [];
    obj.gauge_field              = [];  
    obj.leadmodel                = []; 
    obj.filenameIn = [pwd, filesep, 'Basic_Input_zigzag_leads_orig.xml'];
    obj.filenameOut = [pwd, filesep, 'Basic_Input_running_parameters.xml'];
    obj.WorkingDir = pwd;
    obj.silent     = false;    
    obj.E = 0;
    obj.EF = [];
    obj.q  = [];
    
    
    if strcmpi( class(obj), 'NTerminal') 
        % processing the optional parameters
        obj.InputParsing( varargin{:} );
    
        obj.cCustom_Hamiltonians = [];

        obj.display(['EQuUs:Utils:', class(obj), ':Creating NTerminal class.'])

    
        % exporting calculation parameters into an XML format
        createOutput( obj.filenameOut, obj.Opt, obj.param );            
            
        % create class intances and initializing class attributes
        obj.CreateHandles();
    
        % create Hamiltonian of the scattering region (or load from external source)
        obj.CreateNTerminalHamiltonians();
    end
    


end




%% Transport
%> @brief Calculates the transport through the two terminal setup. Use for development pupose only.
%> @param Energy The energy value.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @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 'gfininvfromHamiltonian' Logical value. Set true calculate the surface Greens function of the scattering region from the Hamiltonaian of the scattering region, or false (default) to calculate it by the fast way (see Phys. Rev. B 90, 125428 (2014) for details).
%> @param 'decimateDyson' Logical value. Set true (default) to decimate the sites of the scattering region in the Dyson equation.
%> @param 'PotInScatter' Obsolete parameter. Use 'scatterPotential' instead.
%> @param 'scatterPotential' A function handle pot=f( #Coordinates ) or pot=f( #CreateHamiltonians, Energy) for the potential to be applied in the Hamiltonian (used when FiniteGreensFunctionFromHamiltonian=true).
%> @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.
%> @return [1] Conductance_matrix The conductance tensor
%> @return [2] ny Array of the open channel in the leads.
%> @return [3] S The scattering matrix.
    function [Conductance_matrix,ny,S] = Transport(obj, Energy, varargin)
        
        p = inputParser;
        p.addParameter('constant_channels', false);
        p.addParameter('gfininvfromHamiltonian', false); 
        p.addParameter('decimateDyson', true);
        p.addParameter('PotInScatter', []) %OBSOLETE use scatterPotential instead
        p.addParameter('scatterPotential', []) %NEW overrides optional argument 'PotInScatter'
        p.addParameter('selfEnergy', false)
        p.parse(varargin{:});
        constant_channels     = p.Results.constant_channels;
        gfininvfromHamiltonian = p.Results.gfininvfromHamiltonian;
        
        scatterPotential               = p.Results.PotInScatter;        
        if ~isempty( p.Results.scatterPotential )
            scatterPotential               = p.Results.scatterPotential;
        end
        
        decimateDyson          = p.Results.decimateDyson;
        selfEnergy             = p.Results.selfEnergy; 
        
        obj.CalcSpectralFunction(Energy, 'constant_channels', constant_channels, 'gfininvfromHamiltonian', gfininvfromHamiltonian, ...
            'decimateDyson', decimateDyson, 'scatterPotential', scatterPotential, 'SelfEnergy', selfEnergy );
        
        
        [S,ny] = obj.FL_handles.SmatrixCalc();
        
        norma = norm(S*S'-eye(sum(ny)));
            
        if norma >= 1e-3
            obj.display( ['error of the unitarity of S-matrix: ',num2str(norma)] )
        end
        
        
        Conductance_matrix = obj.FL_handles.Conduktance();
                
    end


%% CalcSpectralFunction
%> @brief Calculates the spectral density and the Green operator.
%> @param Energy The energy value.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @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 'gfininvfromHamiltonian' Logical value. Set true calculate the surface Greens function of the scattering region from the Hamiltonaian of the scattering region, or false (default) to calculate it by the fast way (see Phys. Rev. B 90, 125428 (2014) for details).
%> @param 'decimateDyson' Logical value. Set true (default) to decimate the sites of the scattering region in the Dyson equation.
%> @param 'PotInScatter' Obsolete parameter. Use 'scatterPotential' instead.
%> @param 'scatterPotential' A function handle pot=f( #Coordinates ) or pot=f( #CreateHamiltonians, Energy) for the potential to be applied in the Hamiltonian (used when FiniteGreensFunctionFromHamiltonian=true).
%> @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.
%> @return [1] The spectral density.
%> @return [2] The Green operator.
    function [A,G] = CalcSpectralFunction(obj, Energy, varargin)
        obj.setEnergy( Energy )
        
        p = inputParser;
        p.addParameter('constant_channels', false);
        p.addParameter('gfininvfromHamiltonian', false);
        p.addParameter('decimateDyson', true);
        p.addParameter('PotInScatter', []) %OBSOLETE use scatterPotential instead
        p.addParameter('scatterPotential', []) %NEW overrides optional argument 'PotInScatter'
        p.addParameter('selfEnergy', false)
        p.parse(varargin{:});
        constant_channels     = p.Results.constant_channels;
        gfininvfromHamiltonian = p.Results.gfininvfromHamiltonian;
        
        scatterPotential               = p.Results.PotInScatter;        
        if ~isempty( p.Results.scatterPotential )
            scatterPotential               = p.Results.scatterPotential;
        end
        
        decimateDyson          = p.Results.decimateDyson;
        selfEnergy             = p.Results.selfEnergy;         
        if ~gfininvfromHamiltonian
            obj.CalcFiniteGreensFunction('gauge_trans', true, 'onlyGinv', true); 
        else
            obj.CalcFiniteGreensFunctionFromHamiltonian('gauge_trans', true, 'scatterPotential', scatterPotential, 'onlyGinv', true); 
        end
        
        Dysonfunc = @()(obj.CustomDysonFunc( 'constant_channels', constant_channels, 'decimate', decimateDyson, 'SelfEnergy', selfEnergy ));
        try
            G = obj.FL_handles.DysonEq( 'CustomDyson', Dysonfunc );
        catch errCause
            err = MException(['EQuUs:Utils:', class(obj),':CalcSpectralFunction'], 'Error occured in calculating the Greens function');
            err = addCause(err, errCause);
            save(['Error_', class(obj), '_CalcSpectralFunction.mat']);
            obj.display(errCause.identifier, 1 );
              obj.display( errCause.stack(1), 1 );            
            throw( errCause );
        end
        
        A = -imag( trace(G))/pi;
        
        
    end

%% getCoordinates
%> @brief Gets the coordinates of the central region
%> @return [1] Coordinates of the central region.
%> @return [2] Coordinates of the interface region.
function [coordinates, coordinates_interface] = getCoordinates( obj )
    try
        coordinates = obj.CreateH.Read('coordinates');
        
    catch errCause
        err = MException(['EQuUs:Utils:', class(obj), ':getCoordinates'], 'Error occured when retriving the coordinates');
        err = addCause(err, errCause);
        save('Error_Ribbon_getCoordinatesFromRibbon.mat');
        throw( err );
    end
    
    
    try
        if isempty( obj.Interface_Regions )
            coordinates_interface = [];
            return
        end
        
        coordinates_interface = cell( length(obj.Interface_Regions), 1);
        for idx = 1:length( obj.Interface_Regions )
            coordinates_interface{idx} = obj.Interface_Regions{idx}.Read('coordinates');
        end
        
    catch errCause
        err = MException(['EQuUs:Utils:', class(obj), ':getCoordinates', 'Error occured when retriving the coordinates']);
        err = addCause(err, errCause);
        save('Error_Ribbon_getCoordinatesFromRibbon.mat');
        throw( err );
    end    
        
end


%% CreateScatter
%> @brief Creates an instance of class #CreateHamiltonians for storing and manipulate the Hamiltonian of the the scattering region. The created object is stored in attribute #CreateH.
    function CreateH = CreateScatter( obj ) 
        
        %> Create an instance of class CreateHamiltonians to create the Hamiltonian of the scattering region
        if ~strcmpi( class(obj.CreateH), 'CreateHamiltonians')
            obj.CreateH = CreateHamiltonians(obj.Opt, obj.param, 'q', obj.q);
        else
            CreateH = obj.CreateH;
        end
        
        %y Create the Hamiltonian of the scattering region if not created
        if ~(obj.CreateH.Read('HamiltoniansCreated') && (~obj.CreateH.Read('HamiltoniansDecimated')))            
            obj.CreateH.CreateScatterH( 'CustomHamiltonian', obj.cCustom_Hamiltonians );    
        end
        
        %> apply magnetic field in scatter
        if ~isempty( obj.PeierlsTransform_Scatter ) && isempty(obj.q) %for finite q vector potential must be parallel to q, and perpendicular to the unit cell vector
               obj.display('Applying magnetic field in the scattering region')
               obj.PeierlsTransform_Scatter.PeierlsTransfor( obj.CreateH );
        end     
        
    end

%% setEnergy
%> @brief Sets the energy for the calculations
%> @param Energy The value of the energy in the same units as the Hamiltonian.
    function setEnergy( obj, Energy )
        obj.E = Energy;

        % resetting the class describing the N-terminal geometry
        if ~isempty( obj.FL_handles ) && strcmpi(class(obj.FL_handles), 'Transport_Interface')        
            obj.FL_handles.setEnergy( obj.E ); 
        end
        
        % reset the class describing the scattering region
        if strcmpi( class(obj.CreateH), 'CreateHamiltonians')
            obj.CreateH.Reset();
        end        
         
        % create interface regions
        for idx = 1:length(obj.Interface_Regions)
          obj.Interface_Regions{idx}.Reset(); 
        end  
        
        % recreate the Hamiltonian of the scattering region
        if ~isempty( obj.CreateH ) && strcmpi( class(obj), 'NTerminal' )      
            obj.CreateNTerminalHamiltonians();
        end
        
    end
    
    
%% getEnergy
%> @brief Retrives the energy value (attribute #E) used in the calculations
%> @return The energy value
    function ret = getEnergy(obj)
      ret = obj.E;    
    end      


%% getFermiEnergy
%> @brief Retrives the Fermi level (attribute #EF) used in the calculations
%> @return The Femi level
    function ret = getFermiEnergy(obj)
      ret = obj.EF;    
    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').
%> @param 'UseHamiltonian' Set true if the interface region is matched to the whole Hamiltonian of the scattering center, or false (default) if the surface Green operator of the scattering center is used in the calculations.
%> @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.addParameter('UseHamiltonian', false); %true if the interface region is matched to the whole Hamiltonian of the scattering center, false if the surface Green operator of the scattering center is used in the calculations.
    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;
    UseHamiltonian    = p.Results.UseHamiltonian;
    
    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), 'UseHamiltonian', UseHamiltonian );
            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       
           
            %> getting dimensions
            rows_interface= zeros(size(Neffs));
            cols_interface= zeros(size(Neffs));
            for idx = 1:length(Neffs)
                rows_interface(idx) = size(K_interface{idx},1);
                cols_interface(idx) = size(K_interface{idx},2);
            end
            cols_scatter = size(gfininv,2);     
            
            % check whether gfininv is sparse
            if issparse(gfininv)
                Ginverz = sparse( [], [], [], sum(Neffs)+sum(rows_interface)+size(gfininv,1), sum(Neffs)+sum(cols_interface)+size(gfininv,2) );
            else
                Ginverz = zeros( sum(Neffs)+sum(rows_interface)+size(gfininv,1), sum(Neffs)+sum(cols_interface)+size(gfininv,2) );
            end
            
            for idx = 1:length(Neffs)
                % surface Green function of leads
                row_indexes_lead = sum(Neffs(1:idx-1))+1:sum(Neffs(1:idx));
                col_indexes_lead = sum(Neffs(1:idx-1))+1:sum(Neffs(1:idx));
                Ginverz( row_indexes_lead, col_indexes_lead ) = Gsurfinverz{idx};
                
                %interface region and coupling to the leads
                row_indexes_interface = sum(Neffs) + (sum(rows_interface(1:idx-1))+1:sum(rows_interface(1:idx)));  
                col_indexes_interface = sum(Neffs) + (sum(cols_interface(1:idx-1))+1:sum(cols_interface(1:idx)));                
                Ginverz( row_indexes_lead, col_indexes_interface ) = -K1{idx};
                Ginverz( row_indexes_interface, col_indexes_lead ) = -K1adj{idx};
                Ginverz( row_indexes_interface, col_indexes_interface ) = -K_interface{idx};
                 
                %scattering region and coupling to the interface regions
                col_indexes_scatter = sum(Neffs) + sum(cols_interface) + (1:cols_scatter);                
                Ginverz( row_indexes_interface, col_indexes_scatter ) = -K1_sc{idx};
                Ginverz( col_indexes_scatter, col_indexes_interface ) = -K1adj_sc{idx};
            end
            
            Ginverz( end-size(gfininv,1)+1:end,  end-size(gfininv,2)+1:end) = gfininv;
            
            %[K0_lead, K1_lead, K1adj_lead] = Leads{1}.Get_Effective_Hamiltonians();
            %Ginverz = [Gsurfinverz{1}, -K1_lead; -K1adj_lead, Gsurfinverz{2}];
        end
        
        
        %-------------------------------
        % calculate the inverse Greens function with the self energy
        function Ginverz = CalcGinverzwithSelfEnergy()
            SelfEnergy = cell(length(Neffs),1);
            for iidx = 1:length(Neffs)
                SelfEnergy{iidx} = Leads{iidx}.Read('Sigma');
            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), 'UseHamiltonian', UseHamiltonian );
            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    
            
            % getting dimensions
            rows_interface= zeros(size(Neffs));
            cols_interface= zeros(size(Neffs));
            for idx = 1:length(Neffs)
                rows_interface(idx) = size(K_interface{idx},1);
                cols_interface(idx) = size(K_interface{idx},2);
            end
            cols_scatter = size(gfininv,2);    
            
            % check whether gfininv is sparse
            if issparse(gfininv)
                Ginverz = sparse( [], [], [], sum(rows_interface)+size(gfininv,1), sum(cols_interface)+size(gfininv,2) );
            else
                Ginverz = zeros( sum(rows_interface)+size(gfininv,1), sum(cols_interface)+size(gfininv,2) );
            end            
            
            for idx = 1:length(Neffs)                
                %interface region with the self energies of the leads
                K_interface{idx}(1:Neffs(idx), 1:Neffs(idx)) = K_interface{idx}(1:Neffs(idx), 1:Neffs(idx)) + SelfEnergy{idx};
                row_indexes_interface = (sum(rows_interface(1:idx-1))+1:sum(rows_interface(1:idx)));  
                col_indexes_interface = (sum(cols_interface(1:idx-1))+1:sum(cols_interface(1:idx)));       
                Ginverz( row_indexes_interface, col_indexes_interface ) = -K_interface{idx};
                 
                %scattering region and coupling to the interface regions
                col_indexes_scatter = sum(cols_interface) + (1:cols_scatter);                
                Ginverz( row_indexes_interface, col_indexes_scatter ) = -K1_sc{idx};
                Ginverz( col_indexes_scatter, col_indexes_interface ) = -K1adj_sc{idx};
            end
            
            Ginverz( end-size(gfininv,1)+1:end,  end-size(gfininv,2)+1:end) = gfininv; 
            
% [K0, K1, K1adj] = Leads{1}.Get_Effective_Hamiltonians(); 
% %[K_interface{1}, K1_interface{1}, K1adj_interface{1}, K1_sc{1}, K1adj_sc{1}] = obj.Interface_Regions{1}.Get_Effective_Hamiltonians();
% Hscatter = obj.CreateH.Read('Hscatter');
% Sscatter = obj.CreateH.Read('Sscatter');
% Kscatter = Hscatter - Sscatter*obj.E;
% 
% Ginverz2 = [-K_interface{1}, -K1, zeros(size(K_interface{2})); 
%           -K1', gfininv, -K1;
%           zeros(size(K_interface{2})), -K1', -K_interface{2}];
%  
% 5

            
        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


%% DecimationFunction
%> @brief Performs the decimation procedure on the inverse Green operator.
%> @param kulso_szabfokok Array of sites to be kept after the decimation.
%> @param ginv The inverse of the Green operator.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'coordinates' An instance of the structure #Coordinates containing the coordinates of the sites.
%> @return [1] The matrix of the decimated inverse Greens operator. 
%> @return [2] An instance of structure #Coordinates containing the sites remained after the decimation.
    function [ret, coordinates_ret] = DecimationFunction( obj, kulso_szabfokok, ginv, varargin)
        p = inputParser;
        p.addParameter('coordinates', [] );        
        p.parse(varargin{:});       
        coordinates    = p.Results.coordinates;
        
        CreateH = CreateHamiltonians(obj.Opt, obj.param);
        
        Hamiltoni = eye(size(ginv))*obj.E - ginv;
        ginv = [];
    
        CreateH.Write('Hscatter', Hamiltoni);
        clear Hamiltoni;
     
        CreateH.Write('kulso_szabfokok', kulso_szabfokok);
        CreateH.Write('coordinates', coordinates);
        CreateH.Write('HamiltoniansCreated', 1);
        CreateH.Write('HamiltoniansDecimated', 0);
        
        Decimation_Procedure = Decimation(obj.Opt); 
        Decimation_Procedure.DecimationFunc(obj.E, CreateH, 'Hscatter', 'kulso_szabfokok');
    
        Hamiltoni = CreateH.Read('Hscatter');
        CreateH.Clear('Hscatter');
    
        ret = eye(size(Hamiltoni))*obj.E - Hamiltoni;
        Hamiltoni = [];
        
        
        coordinates_ret = CreateH.Read('coordinates');
        
        
    end


%% GetFiniteGreensFunction
%> @brief Query fo the calculated Green operator of the scattering center.
%> @return [1] The Green operator.
%> @return [2] The inverse Green operator.
    function [Gret, Ginvret] = GetFiniteGreensFunction( obj )
        Gret = obj.G;
        Ginvret = obj.Ginv;
    end

%% CalcFiniteGreensFunction
%> @brief Calculates the Green operator of the scattering region.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'gauge_trans' Logical value. Set true to perform gauge transformation on the Green operator and on the Hamiltonians.
%> @param 'onlyGinv' Logical value. Set true to calculate only the inverse of the surface Greens function #Ginv, or false (default) to calculate #G as well. In the latter case the attribute #Ginv is set to empty at the end.
    function CalcFiniteGreensFunction( obj, varargin )
    
     p = inputParser;
        p.addParameter('gauge_trans', false); % logical: true if want to perform gauge transformation on the Green's function and Hamiltonians 
        p.addParameter('onlyGinv', false);
        p.parse(varargin{:});
        gauge_trans     = p.Results.gauge_trans;
        onlyGinv        = p.Results.onlyGinv;     
        
          obj.CalcFiniteGreensFunctionFromHamiltonian('gauge_trans', gauge_trans, 'onlyGinv', onlyGinv);
        
    end


%% CalcFiniteGreensFunctionFromHamiltonian
%> @brief Calculates the Green operator of the scattering region from the whole Hamiltonian.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'gauge_trans' Logical value. Set true to perform gauge transformation on the Green operator and on the Hamiltonians.
%> @param 'onlyGinv' Logical value. Set true to calculate only the inverse of the surface Greens function #Ginv, or false (default) to calculate #G as well. In the latter case the attribute #Ginv is set to empty at the end.
%> @param 'PotInScatter' Obsolete parameter. Use 'scatterPotential' instead.
%> @param 'scatterPotential' A function handle pot=f( #Coordinates ) or pot=f( #CreateHamiltonians, Energy) for the potential to be applied in the Hamiltonian (used when FiniteGreensFunctionFromHamiltonian=true). 
    function CalcFiniteGreensFunctionFromHamiltonian( obj, varargin )
        
        
        
     p = inputParser;
        p.addParameter('gauge_trans', false); % logical: true if want to perform gauge transformation on the Green's function and Hamiltonians     
        p.addParameter('onlyGinv', false);
        p.addParameter('PotInScatter', []) %OBSOLETE use scatterPotential instead
        p.addParameter('scatterPotential', []) %NEW overrides optional argument 'PotInScatter'
        p.parse(varargin{:});
        gauge_trans              = p.Results.gauge_trans;
        onlyGinv                 = p.Results.onlyGinv;
        
        
        scatterPotential               = p.Results.PotInScatter;        
        if ~isempty( p.Results.scatterPotential )
            scatterPotential               = p.Results.scatterPotential;
        end
        
        
        obj.CreateScatter();
%**************************************** 
        CreateH = obj.CreateH;
        Hscatter = CreateH.Read('Hscatter');
        Sscatter = CreateH.Read('Sscatter');
        coordinates = CreateH.Read('coordinates');
        kulso_szabfokok = CreateH.Read( 'kulso_szabfokok' );
        
        
        if ~isempty(scatterPotential)          
            obj.ApplyPotentialInScatter( CreateH, scatterPotential)  
            Hscatter = CreateH.Read('Hscatter');
        end     
                        
              
        q = obj.CreateH.Read('q');
        if ~isempty( q ) && ~obj.CreateH.Read('HamiltoniansDecimated')
            Hscatter_transverse = obj.CreateH.Read('Hscatter_transverse');
            Hscatter = Hscatter + Hscatter_transverse*diag(exp(-1i*q)) + Hscatter_transverse'*diag(exp(1i*q)); 
            obj.CreateH.Clear('Hscatter_transverse');
        end
        
        % Applying overlap matrices
        if isempty( Sscatter )
            Hscatter = (sparse(1:size(Hscatter,1), 1:size(Hscatter,1), obj.E, size(Hscatter,1),size(Hscatter,1))   - Hscatter);
        else
            Hscatter = obj.E*Sscatter - Hscatter;
        end      

        indexes = false( size(Hscatter,1), 1);
        indexes(kulso_szabfokok) = true;
        Hscatter = [ Hscatter( ~indexes, ~indexes) , Hscatter(~indexes, indexes);
                     Hscatter(indexes, ~indexes), Hscatter(indexes, indexes)];

        obj.G = obj.partialInv( Hscatter, length(kulso_szabfokok) );        

          
%         Hscatter = obj.CreateH.Read('Hscatter');
%         Hscatter = (sparse(1:size(Hscatter,1), 1:size(Hscatter,1), obj.E, size(Hscatter,1),size(Hscatter,1))   - Hscatter);
%         obj.G = inv(Hscatter);
%**************************************************        
        
        % gauge transformation of the vector potential in the effective Hamiltonians
        if gauge_trans
            try
                surface_coordinates = obj.getCoordinates();
               % gauge transformation on Green's function
                if ~isempty(obj.G) && ~isempty(obj.PeierlsTransform_Scatter) && isempty(obj.q)
                    % gauge transformation on the inverse Green's function
                    obj.G = obj.PeierlsTransform_Scatter.gaugeTransformation( obj.G, surface_coordinates, obj.gauge_field );
                end
            catch  errCause  
                err = MException(['EQuUs:Utils:', class(obj), ':CalcFiniteGreensFunctionFromHamiltonian'], 'Unable to perform gauge transformation');
                err = addCause(err, errCause);
                save('Error_Ribbon_CalcFiniteGreensFunctionFromHamiltonian.mat')
                throw(err);
            end
        end
              
        if onlyGinv
            rcond_G = rcond(obj.G);
            if isnan(rcond_G ) || abs( rcond_G ) < 1e-15
                obj.display( ['EQuUs:Utils:', class(obj), ':CalcFiniteGreensFunctionFromHamiltonian: Regularizing Ginv by SVD'],1);
                obj.Ginv = obj.inv_SVD( obj.G );
            else
                obj.Ginv = inv(obj.G);
            end       
            obj.G = [];
        else
            obj.Ginv = [];
        end
        
    end




%% setInterfaceRegions
%> @brief Replaces the attribute Interface_Region with the given value.
%> @param Interface_Regions A two component array of classes InterfaceRegion 
    function setInterfaceRegions( obj, Interface_Regions )
       if length( Interface_Regions ) ~= length(obj.param.Leads)
           err = MException(['EQuUs:Utils:', class(obj), ':setInterfaceRegions'], ['The length of Interface_Regions must be ', num2str(length(obj.param.Leads))]);
           save('Error_Ribbon_setInterfaceRegions.mat');
           throw(err);          
       end
       
       if ~iscell( Interface_Regions )
           err = MException(['EQuUs:Utils:', class(obj), ':setInterfaceRegions'], 'The Interface_Regions must be a cell containing of instances of class InterfaceRegion.');
           save('Error_Ribbon_setInterfaceRegions.mat');
           throw(err);          
       end 
       
       obj.Interface_Regions = Interface_Regions;
        
    end

%% CreateInterface
%> @brief Creates the Hamiltonians for the interface regions between the leads and scattering center.
%> @param idx Identification number of the interface region. 
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'UseHamiltonian' Set true if the Dyson equation is constructed from the whole Hamiltonian of the scattering region and the interface region is not needed to be truncated
     function CreateInterface( obj, idx, varargin )
        
        p = inputParser;
        p.addParameter('UseHamiltonian', false); %Set true if the Dyson equation is constructed from the whole Hamiltonian of the scattering region and the interface region is not needed to be truncated
        p.parse(varargin{:});
        UseHamiltonian    = p.Results.UseHamiltonian;
        
        Interface_Region = obj.Interface_Regions{idx} ;

        if ~obj.Interface_Regions{idx}.Read( 'HamiltoniansCreated' )    
          obj.display(['EQuUs:Utils:', class(obj), ':CreateInterface: Creating interface region ', num2str(idx)]) 
            H0 = obj.cCustom_Hamiltonians.Read( 'H0' );
            S0 = obj.cCustom_Hamiltonians.Read( 'S0' );
            H1 = obj.cCustom_Hamiltonians.Read( 'H1' );
            S1 = obj.cCustom_Hamiltonians.Read( 'S1' );
            H_coupling = obj.cCustom_Hamiltonians.Read( 'Hcoupling' );
            S_coupling = obj.cCustom_Hamiltonians.Read( 'Scoupling' );
            coordinates = obj.cCustom_Hamiltonians.Read( 'coordinates' );
            Interface_Region.Write( 'H0', H0{idx} );
            if ~isempty( S0 )
                Interface_Region.Write( 'S0', S0{idx} );
            end

            Interface_Region.Write( 'H1', H1{idx} );
            Interface_Region.Write( 'H1adj', H1{idx}' ); 
            if ~isempty( S1 )
                Interface_Region.Write( 'S1', S1{idx} );  
            end  
            

            
            Interface_Region.Write('Hcoupling', H_coupling{idx});
            Interface_Region.Write('Hcouplingadj', H_coupling{idx}');
            if ~isempty( S_coupling )
                Interface_Region.Write('Scoupling', S_coupling{idx});     
            end
            
            Interface_Region.Write( 'M', size(H0{idx},1) );  
            Interface_Region.Write( 'coordinates', coordinates{idx} ); 
            Interface_Region.Write( 'coordinates2', obj.CreateH.Read('coordinates') );
            
            % Transforming the Hamiltonians according to the BdG model
            Interface_Region.Transform2BdG();
            
            % truncating the coupling Hamiltonians to fit the scattering region
            if ~UseHamiltonian
                surface_points = obj.CreateH.Read( 'kulso_szabfokok' );
                H_coupling = Interface_Region.Read('Hcoupling');   
                indexes = false( size(H_coupling, 2), 1);
                indexes( surface_points ) = true;
                H_coupling = H_coupling(:,indexes);            
                Interface_Region.Write('Hcoupling', H_coupling);
                Interface_Region.Write('Hcouplingadj', H_coupling');
            
            
                S_coupling = Interface_Region.Read('Scoupling');
                if ~isempty( S_coupling )
                    S_coupling = S_coupling(:,surface_points);
                    Interface_Region.Write('Scoupling', S_coupling);    
                end
            
                % keeping data on the surface sites of the scattering region
                coordinates2 = obj.Interface_Regions{idx}.Read('coordinates2');             
                coordinates2 = coordinates2.KeepSites( indexes );  
                obj.Interface_Regions{idx}.Write('coordinates2', coordinates2);
            
            end
                               
        end
        
        
        % apply magnetic field in the interface region
        if ~isempty( obj.PeierlsTransform_Leads )      
            % In superconducting lead one must not include nonzero magnetic
            % field.
            % Hamiltonians in transverse computations must remain
            % traslational invariant. 
            if ~Interface_Region.isSuperconducting() %&& isempty( obj.Interface_Regions{idx}.Read('q') )
                obj.display(['EQuUs:Utils:', class(obj), 'CreateInterface: Applying magnetic field in the Hamiltonians of interface region ', num2str(idx)])
                    obj.PeierlsTransform_Leads.PeierlsTransformLeads( Interface_Region );  
            else %&& isempty( obj.Interface_Regions{idx}.Read('q') )
                obj.display(['EQuUs:Utils:', class(obj), 'CreateInterface: Applying gauge transformation in the Hamiltonians of interface region ', num2str(idx)])
                    obj.PeierlsTransform_Leads.gaugeTransformationOnLead( Interface_Region, obj.gauge_field );
            end
        end  
        
        Interface_Region.ApplyOverlapMatrices( obj.E );    
        
        % method to adjust the interface region and coupling to the scattering region by an external function.
        if ~isempty( obj.interfacemodel )
            obj.interfacemodel( Interface_Region ); 
        end
        
        % The regularization of the interface is performed according to the Leads
        Leads = obj.FL_handles.Read( 'Leads' );
        Lead = Leads{idx};
        Interface_Region.Calc_Effective_Hamiltonians( obj.E, 'Lead', Lead );
    end
    





%% setHandlesForMagneticField
%> @brief Sets the function handles of the vector potential and gauge transformation.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'scatter' Function handle A = f( x,y) of the vector potential to be used in the scattering region (A is a N x 2 vector, where N is the number of the points given by the x and y coordinates.)
%> @param 'lead' Function handle A = f( x,y) of the vector potential to be used in the leads. (A is a N x 2 vector, where N is the number of the points given by the x and y coordinates.)
%> @param 'gauge_field' Function handle S = f( x,y) of the gauge transformation. (S is a N x 1 vector, where N is the number of the points given by the x and y coordinates.)
    function setHandlesForMagneticField( obj, varargin )     
    
        p = inputParser;
        p.addParameter('scatter', [] );
        p.addParameter('lead', [] );
        p.addParameter('gauge_field', [] ); 
        
        p.parse(varargin{:});

       
        hscatter    = p.Results.scatter;
        hlead       = p.Results.lead;
        obj.gauge_field = p.Results.gauge_field;
        
        
        obj.display(['EQuUs:Utils:', class(obj), ':setHandlesForMagneticField: Adding handles of the magnetic field to the scattering region']) 
        if obj.Opt.magnetic_field && isempty( obj.PeierlsTransform_Scatter )
            % setting the vector potential for the scattering region
          obj.PeierlsTransform_Scatter = Peierls(obj.Opt, 'Vectorpotential', hscatter);
        elseif ~isempty( obj.PeierlsTransform_Scatter ) && ~isempty(hscatter)
            % setting the vector potential for the scattering region if the Peierls class was already constructed
          obj.PeierlsTransform_Scatter.setVectorPotential( hscatter );
        end
        
        if obj.Opt.magnetic_field && isempty( obj.PeierlsTransform_Leads )
            % setting the vector potential for the leads
          obj.PeierlsTransform_Leads   = Peierls(obj.Opt, 'Vectorpotential', hlead);
        elseif ~isempty( obj.PeierlsTransform_Leads ) && ~isempty(hlead)
            % setting the vector potential for the leads if the Peierls class was already constructed
          obj.PeierlsTransform_Leads.setVectorPotential( hlead );
          end
          
    end



%% CreateClone
%> @brief Creates a clone of the present object.
%> @return Returns with the cloned object.
    function ret = CreateClone( obj )
        
        % cerating new instance of class NTerminal
        ret = NTerminal( ...
            'filenameIn', obj.filenameIn, ...
            'filenameOut', obj.filenameOut, ...
            'E', obj.E, ...
            'EF', obj.EF, ...
            'phi', obj.phi, ...
            'silent', obj.silent, ...
               'Opt', obj.Opt, ...
               'param', obj.param, ...
            'q', obj.q, ...
            'leadmodel', obj.leadmodel, ...
            'interfacemodel', obj.interfacemodel );
        
        %> cloning the individual attributes
        ret.CreateH = obj.CreateH.CreateClone();
        ret.FL_handles = obj.FL_handles.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
        
        %> setting the gauge field
        ret.gauge_field              = obj.gauge_field;  % function handle for the scalar field to transform the vector potential from Landauy to Landaux
        
    end   
    
   
    
%% setParam
%> @brief Sets the structure #param in the attributes.
%> @param param An instance of structure #param.
    function setParam(obj, param)
        obj.param = param;    
      
        if ~isempty( obj.CreateH )
            obj.CreateH.Write('param', param);
        end
        
        if ~isempty( obj.FL_handles )
            obj.FL_handles.Write('param', param);
        end
        
        if ~isempty( obj.cCustom_Hamiltonians )
            obj.cCustom_Hamiltonians.Write('param', param);
        end
        
        for idx = 1:length( obj.Interface_Regions )
            obj.Interface_Regions{idx}.Write('param', param);
        end

    end  
    
%% getParam
%> @brief Retrives a copy of the structure #param used in the current calculations.
%> @return param An instance of structure #param.
    function ret = getParam(obj)
      ret = obj.param;    
    end    
    
%% getTransverseMomentum
%> @brief Retrives the value of the transverse momentum quantum number.
%> @return q The transverse momentum quantum number.
    function ret = getTransverseMomentum(obj)
      ret = obj.q;    
    end          
    

%% ApplyPotentialInScatter
%> @brief Applies the potential in the scattering region
%> @param CreateH An instance of class #CreateHamiltonians containing the Hamiltonian of the scattering center.
%> @param scatterPotential A function handle pot=f( #Coordinates ) or pot=f( #CreateHamiltonians, Energy) for the potential to be applied in the Hamiltonian (used when FiniteGreensFunctionFromHamiltonian=true). 
    function ApplyPotentialInScatter( obj, CreateH, scatterPotential )
        

        if nargin( scatterPotential ) == 1
            %> calculating the potential if the scattering potential has one input arguments
            %> obtaining coordinates of the scattering region
            coordinates = CreateH.Read('coordinates');        
            %> calculating the potential from the coordinates
            potential = scatterPotential( coordinates );
        elseif nargin( scatterPotential ) == 2
            %> calculating the potential if the scattering potential has two input arguments
            potential = scatterPotential( CreateH, obj.E );
            %> obtaining coordinates of the scattering region that might have been changed inside the scatterPotential function
            coordinates = CreateH.Read('coordinates');   % coordinates might be changed in function scatterPotential
        else
            error('EQuUs:Ribbon:ApplyPotentialInScatter', 'To many input arguments in function handle scatterPotential');
        end
        
        % obtaining the scattering Hamiltonian
        Hscatter = CreateH.Read('Hscatter'); 
        
        if isvector(potential) && length(potential) == size(Hscatter,1)
            % applying potential in case of a diagonal on-site potential
            
            % transforming the potential for the case of the BdG model.
            if isprop(coordinates, 'BdG_u')
                potential(~coordinates.BdG_u) = -potential(~coordinates.BdG_u);
            end
            
            %creating diagonal sparse matrix from the potential
               potential = sparse(1:size(Hscatter,1), 1:size(Hscatter,1), potential, size(Hscatter,1), size(Hscatter,1));
            
        elseif length(size(potential)) == 2 && norm( size(potential)-size(Hscatter) ) == 0
            % applying potential in case, when the potential is given in a matrix form            
            
            % transforming the potential for the case of the BdG model.
            if isprop(coordinates, 'BdG_u')
                potential(~coordinates.BdG_u, ~coordinates.BdG_u) = -potential(~coordinates.BdG_u, ~coordinates.BdG_u);
            end
        else
            error(['EQuUs:utils:', class(obj), ':ApplyPotentialInScatter'], 'Wrong output format of the function handle scatterPotential');
        end
        
        
        % Ensure not to overload the memory
        if issparse(Hscatter) && ~issparse(potential)
            error(['EQuUs:utils:', class(obj), ':ApplyPotentialInScatter'], 'If the Hamiltonian is sparse, potential should also be sparse');
        end
        
        % adding potential to the scattering region
        Hscatter = Hscatter + potential;

        % save the Hamiltonain 
        CreateH.Write('Hscatter', Hscatter);

    end        
      
    
end %methods public    

methods ( Access = protected )
    
%% CreateNTerminalHamiltonians
%> @brief Extracts the Hamiltonians from the external source.
    function CreateNTerminalHamiltonians( obj )
        
        % Creating attribute Custom_Hamiltonian providing Hamiltonians from external source
        if ~strcmpi( class(obj.cCustom_Hamiltonians), 'Custom_Hamiltonian' )
            obj.cCustom_Hamiltonians =  Custom_Hamiltonians( obj.Opt, obj.param, 'EF', obj.EF, 'CustomHamiltoniansHandle', obj.CustomHamiltoniansHandle );
        end
        
        % load Hamiltonians from the external source
        if ~obj.cCustom_Hamiltonians.Read( 'Hamiltonians_loaded' )
            obj.cCustom_Hamiltonians.LoadHamiltonians( 'WorkingDir', obj.WorkingDir );
        end
        
        % cerating the Hamiltonain for the scattering region
        obj.CreateScatter();
    end
   
     
%% CreateHandles
%> @brief Initializes the attributes of the class.
    function CreateHandles( obj )  
        
        % creating classes for managing the magnetic field
        obj.setHandlesForMagneticField();
        
        % create class storing the Hamiltonian of the scattering region
        obj.CreateH = CreateHamiltonians(obj.Opt, obj.param, 'q', obj.q);
        
        % create class for transport calculations
        obj.FL_handles = Transport_Interface(obj.E, obj.Opt, obj.param, 'CreateH', obj.CreateH, 'PeierlsTransform', obj.PeierlsTransform_Leads );
        
        % read out structure Opt containing the computational parameters
        obj.Opt = obj.FL_handles.Read( 'Opt' );
        
        % creating classes of the interface regions
        obj.Interface_Regions = cell(length(obj.param.Leads),1);
        for idx = 1:length(obj.Interface_Regions)
            Lead_Orientation = obj.param.Leads{idx}.Lead_Orientation;
            obj.Interface_Regions{idx} = InterfaceRegion(obj.Opt, obj.param, 'Hanyadik_Lead', idx, 'q', obj.q, 'Lead_Orientation', Lead_Orientation);
        end
        
    end

    
end % protected methods
       
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 '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 'CustomHamiltoniansHandle' function handle for the custom Hamiltonians. Has the same inputs as #Custom_Hamiltonians.LoadHamiltonians and output values defined by the example #Hamiltonians.
%> @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 '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)
    
        p = inputParser;
        p.addParameter('filenameIn', obj.filenameIn, @ischar);
        p.addParameter('filenameOut', obj.filenameOut, @ischar);
        p.addParameter('WorkingDir', obj.WorkingDir, @ischar);
        p.addParameter('CustomHamiltoniansHandle', obj.CustomHamiltoniansHandle); %function handle for the custom Hamiltonians
        p.addParameter('E', obj.E, @isscalar);
        p.addParameter('EF', obj.EF); 
        p.addParameter('silent', obj.silent);   
        p.addParameter('leadmodel', obj.leadmodel); %individual physical model for the contacts
        p.addParameter('interfacemodel', obj.interfacemodel); %individual physical model for the interface regions
        p.addParameter('Opt', obj.Opt);
        p.addParameter('param', obj.param);
        p.addParameter('q', obj.q);
        
        p.parse(varargin{:});

       
        obj.filenameIn   = p.Results.filenameIn;        
        obj.filenameOut  = p.Results.filenameOut;
        obj.WorkingDir   = p.Results.WorkingDir;
        obj.CustomHamiltoniansHandle = p.Results.CustomHamiltoniansHandle;
        obj.E            = p.Results.E;
        obj.EF           = p.Results.EF;
        obj.silent       = p.Results.silent;
        obj.leadmodel    = p.Results.leadmodel;
        obj.interfacemodel  = p.Results.interfacemodel;
        obj.q            = p.Results.q;
          obj.Opt          = p.Results.Opt;
          obj.param        = p.Results.param;       
        
                
          if isempty(obj.Opt) || isempty(obj.param)   
               [ obj.Opt, obj.param ] = parseInput( obj.filenameIn );
          end

    end

end % methods private

end